include/Pythia8/DireSplittings.h

0001 // DireSplittings.h is a part of the PYTHIA event generator.
0002 // Copyright (C) 2024 Stefan Prestel, Torbjorn Sjostrand.
0003 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
0004 // Please respect the MCnet Guidelines, see GUIDELINES for details.
0005
0006 // Header file for the OverheadInfo and DireSplitting classes.
0007
0008 #ifndef Pythia8_DireSplittings_H
0009 #define Pythia8_DireSplittings_H
0010
0011 #define DIRE_SPLITTINGS_VERSION "2.002"
0012
0013 #include "Pythia8/Basics.h"
0014 #include "Pythia8/BeamParticle.h"
0015 #include "Pythia8/ParticleData.h"
0016 #include "Pythia8/PythiaStdlib.h"
0017 #include "Pythia8/Settings.h"
0018 #include "Pythia8/StandardModel.h"
0019 #include "Pythia8/Info.h"
0020 #include "Pythia8/DireSplitInfo.h"
0021 #include "Pythia8/DireBasics.h"
0022
0023 namespace Pythia8 {
0024
0025 class DireSpace;
0026 class DireTimes;
0027
0028 //==========================================================================
0029
0030 class OverheadInfo {
0031
0032 public:
0033
0034   OverheadInfo(int nFinalIn, int idIn, double valIn, double xIn, double pT2In)
0035     : nFinal(nFinalIn), id(idIn), val(valIn), x(xIn), pT2(pT2In) {}
0036
0037   int nFinal, id;
0038   double val, x, pT2;
0039
0040   bool match(int idIn, int nfIn) { return (idIn==id && nfIn==nFinal); }
0041
0042   string list () const {
0043     ostringstream os;
0044     os << scientific << setprecision(6)
0045     << "pT2 " << setw(10) << pT2 << " x " << setw(10) << x
0046     << " id " << setw(4) << id << " nf " << setw(4) << nFinal
0047     << " val=" << val;
0048     return os.str();
0049   }
0050
0051 };
0052
0053 class DireSplitting {
0054
0055 public:
0056
0057   // Constructor and destructor.
0058   DireSplitting() : renormMultFac(0),
0059       id("void"), correctionOrder(0), settingsPtr(0),
0060       particleDataPtr(0), rndmPtr(0), beamAPtr(0),
0061       beamBPtr(0),  coupSMPtr(0), infoPtr(0), direInfoPtr(0),
0062       is_qcd(false), is_qed(false), is_ewk(false), is_fsr(false),
0063       is_isr(false), is_dire(false), nameHash(0) {}
0064   DireSplitting(string idIn, int softRS, Settings* settings,
0065     ParticleData* particleData, Rndm* rndm, BeamParticle* beamA,
0066     BeamParticle* beamB, CoupSM* coupSMPtrIn, Info* infoPtrIn,
0067                 DireInfo* direInfo) :
0068       renormMultFac(0), id(idIn), correctionOrder(softRS),
0069       settingsPtr(settings), particleDataPtr(particleData), rndmPtr(rndm),
0070       beamAPtr(beamA), beamBPtr(beamB), coupSMPtr(coupSMPtrIn),
0071       infoPtr(infoPtrIn), direInfoPtr(direInfo), is_qcd(false), is_qed(false),
0072       is_ewk(false), is_fsr(false), is_isr(false), is_dire(false),
0073        nameHash(0) { init(); splitInfo.storeName(name()); }
0074   virtual ~DireSplitting() {}
0075
0076   void init();
0077
0078 public:
0079
0080   double renormMultFac;
0081
0082   string id;
0083   int correctionOrder;
0084   Settings* settingsPtr;
0085   ParticleData* particleDataPtr;
0086   Rndm* rndmPtr;
0087   BeamParticle* beamAPtr;
0088   BeamParticle* beamBPtr;
0089   CoupSM* coupSMPtr;
0090   Info* infoPtr;
0091   DireInfo* direInfoPtr;
0092
0093   // Some short-cuts and string hashes to help avoid string comparisons.
0094   bool is_qcd, is_qed, is_ewk, is_fsr, is_isr, is_dire;
0095   ulong nameHash;
0096   bool is( ulong pattern ) {
0097     if (pattern == nameHash) return true;
0098     return false;
0099   }
0100
0101   unordered_map<string,double> kernelVals;
0102
0103   string name () {return id;}
0104
0105   virtual bool canRadiate ( const Event&, pair<int,int>,
0106     unordered_map<string,bool> = unordered_map<string,bool>(),
0107     Settings* = NULL, PartonSystems* = NULL, BeamParticle* = NULL)
0108     {return false;}
0109
0110   // Discard below the cut-off for the splitting.
0111   virtual bool aboveCutoff( double, const Particle&, const Particle&, int,
0112     PartonSystems* = NULL) { return true; }
0113
0114   virtual bool useFastFunctions() { return false; }
0115   virtual bool canRadiate ( const Event&, int, int,
0116     Settings* = NULL, PartonSystems* = NULL, BeamParticle* = NULL)
0117     {return false;}
0118
0119   // Function to return an identifier for the phase space mapping
0120   // that is used for setting up this splitting.
0121   // return values: 1 --> Default Dire mapping.
0122   //                2 --> Dire 1->3 mapping.
0123   virtual int kinMap () {return 1;}
0124
0125   // Return id of mother after splitting.
0126   virtual int motherID(int) {return 0;}
0127
0128   // Return id of emission.
0129   virtual int sisterID(int) {return 0;}
0130
0131   // Return a pair of ids for the radiator and emission after
0132   // the splitting.
0133   virtual vector <int> radAndEmt(int, int) { return vector<int>(); }
0134   virtual vector < pair<int,int> > radAndEmtCols( int, int, Event)
0135     { return vector<pair<int,int> >(); }
0136   virtual bool canUseForBranching() { return false; }
0137   virtual bool isPartial()  { return false; }
0138   virtual int  nEmissions() { return 0; }
0139
0140   virtual bool swapRadEmt() { return false; }
0141   virtual bool isSymmetric( const Particle* = NULL, const Particle* = NULL)
0142     { return false; }
0143
0144   // Return a vector of all possible recoiler positions, given the
0145   // positions of the radiator and emission after the splitting.
0146   virtual vector <int> recPositions( const Event&, int, int)
0147     {return vector<int>();}
0148
0149   // Return id of recombined radiator (before splitting!)
0150   virtual int radBefID(int, int) {return 0;}
0151
0152   // Return colours of recombined radiator (before splitting!)
0153   virtual pair<int,int> radBefCols(int, int, int, int)
0154     {return make_pair(0,0);}
0155
0156   // Return color factor for splitting.
0157   virtual double gaugeFactor (int, int) {return 1.;}
0158
0159   // Return symmetry factor for splitting.
0160   virtual double symmetryFactor (int, int) {return 1.;}
0161
0162   // Return an identifier for the interaction that causes the
0163   // branching.
0164   // return values -1 --> Type not defined.
0165   //                1 --> QCD splitting (i.e. proportional to alphaS)
0166   //                2 --> QED splitting (i.e. proportional to alphaEM)
0167   //                3 --> EW splitting (i.e. proportional to sinThetaW)
0168   //                4 --> Yukawa splitting (i.e. proportional to y)
0169   virtual int couplingType (int, int) {return -1;}
0170
0171   // Return the value of the coupling that should be used for this branching.
0172   // Note that the last input allows easy access to the PS evolution variable.
0173   // return values -1         --> Coupling value not defined.
0174   //               double > 0 --> Value to be used for this branching.
0175   virtual double coupling (double = 0., double = 0., double = 0.,
0176     double = -1., pair<int,bool> = pair<int,bool>(),
0177     pair<int,bool> = pair<int,bool>()) {
0178     return -1.;
0179   }
0180   virtual double couplingScale2 (double = 0., double = 0., double = 0.,
0181     pair<int,bool> = pair<int,bool>(), pair<int,bool> = pair<int,bool>()) {
0182     return -1.;
0183   }
0184
0185   // Pick z for new splitting.
0186   virtual double zSplit(double, double, double) {return 0.5;}
0187
0188   // New overestimates, z-integrated versions.
0189   virtual double overestimateInt(double, double, double, double, int = -1)
0190     { return 0.;}
0191
0192   // Return kernel for new splitting.
0193   virtual double overestimateDiff(double, double, int = -1) {return 1.;}
0194
0195   // Functions to store and retrieve all the variants of the kernel.
0196   virtual double getKernel(string = "");
0197   virtual unordered_map<string,double> getKernelVals() { return kernelVals; }
0198   virtual void   clearKernels()         { kernelVals.clear(); }
0199
0200   DireSplitInfo splitInfo;
0201
0202   // Functions to calculate the kernel from SplitInfo information.
0203   virtual bool calc(const Event& = Event(), int = -1) { return false; }
0204
0205   shared_ptr<DireSpace> isr;
0206   shared_ptr<DireTimes> fsr;
0207   shared_ptr<DireTimes> fsrDec;
0208   void setTimesPtr(shared_ptr<DireTimes> fsrIn)    { fsr=fsrIn;}
0209   void setTimesDecPtr(shared_ptr<DireTimes> fsrIn) { fsrDec=fsrIn;}
0210   void setSpacePtr(shared_ptr<DireSpace> isrIn)    { isr=isrIn;}
0211
0212   // Functions that allow different ordering variables for emissions.
0213   // Note: Only works after splitInfo has been properly filled.
0214   virtual double getJacobian( const Event& = Event(),
0215     PartonSystems* = 0) { return 0.;}
0216   // Map filled identical to shower state variables map.
0217   virtual unordered_map<string, double> getPhasespaceVars(
0218     const Event& = Event(),
0219     PartonSystems* = 0) { return unordered_map<string,double>(); }
0220
0221   // Treatment of additional virtual corrections.
0222   virtual bool allow_z_endpoint_for_kinematics()   { return false; }
0223   virtual bool allow_pT2_endpoint_for_kinematics() { return false; }
0224   virtual bool allow_sai_endpoint_for_kinematics() { return false; }
0225   virtual bool allow_xa_endpoint_for_kinematics()  { return false; }
0226   // Functions to set if kernel should contribute to a kinematical endpoint.
0227   virtual void try_z_endpoint()                    { return; }
0228   virtual void try_pT2_endpoint()                  { return; }
0229   virtual void try_sai_endpoint()                  { return; }
0230   virtual void try_xa_endpoint()                   { return; }
0231   // Return endpoint information.
0232   virtual bool is_z_endpoint()                     { return false; }
0233   virtual bool is_pT2_endpoint()                   { return false; }
0234   virtual bool is_sai_endpoint()                   { return false; }
0235   virtual bool is_xa_endpoint()                    { return false; }
0236
0237   // Functions to calculate Dire variables from the evolution variables.
0238   virtual double tdire_ff(double, double t, double) { return t; }
0239   virtual double zdire_ff(double z, double, double) { return z; }
0240   virtual double tdire_fi(double, double t, double) { return t; }
0241   virtual double zdire_fi(double z, double, double) { return z; }
0242   virtual double tdire_if(double, double t, double) { return t; }
0243   virtual double zdire_if(double z, double, double) { return z; }
0244   virtual double tdire_ii(double, double t, double) { return t; }
0245   virtual double zdire_ii(double z, double, double) { return z; }
0246
0247   virtual bool hasMECBef(const Event&, double) { return false; }
0248   virtual bool hasMECAft(const Event&, double) { return false; }
0249
0250   virtual void storeOverhead(double pT2, double x, int radid, int nf,
0251     double val) { overhead_map.insert(make_pair(pT2, OverheadInfo(nf, radid,
0252     val, x, pT2))); }
0253   multimap<double,OverheadInfo> overhead_map;
0254
0255   virtual double overhead(double pT2, int idd, int nf) {
0256
0257     if (overhead_map.empty()) return 1.;
0258
0259     multimap<double,OverheadInfo>::iterator lo = overhead_map.lower_bound(pT2);
0260     if (lo != overhead_map.begin()) lo--;
0261     if (lo != overhead_map.begin()) lo--;
0262     multimap<double,OverheadInfo>::iterator hi = overhead_map.upper_bound(pT2);
0263     if (hi != overhead_map.end()) hi++;
0264     if (hi == overhead_map.end()) hi--;
0265
0266     int n(0);
0267     double sum = 0.;
0268     for ( multimap<double,OverheadInfo>::iterator it = lo;
0269       it != hi; ++it ){
0270       if (!it->second.match(idd,nf)) continue;
0271       sum += it->second.val;
0272       n++;
0273     }
0274
0275     if (hi->second.match(idd,nf)) {
0276       sum += hi->second.val;
0277       n++;
0278     }
0279
0280     return max(sum/max(1,n),1.);
0281
0282   }
0283
0284 };
0285
0286 //==========================================================================
0287
0288 } // end namespace Pythia8
0289
0290 #endif // Pythia8_DireSplittings_H