include/Pythia8/DireSpace.h

0001 // DireSpace.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 spacelike initial-state showers.
0007 // DireSpaceEnd: radiating dipole end in ISR.
0008 // DireSpace: handles the showering description.
0009
0010 #ifndef Pythia8_DireSpace_H
0011 #define Pythia8_DireSpace_H
0012
0013 #define DIRE_SPACE_VERSION "2.002"
0014
0015 #include "Pythia8/Basics.h"
0016 #include "Pythia8/SpaceShower.h"
0017 #include "Pythia8/BeamParticle.h"
0018 #include "Pythia8/Event.h"
0019 #include "Pythia8/Info.h"
0020 #include "Pythia8/ParticleData.h"
0021 #include "Pythia8/PartonSystems.h"
0022 #include "Pythia8/PythiaStdlib.h"
0023 #include "Pythia8/Settings.h"
0024 #include "Pythia8/StandardModel.h"
0025 #include "Pythia8/UserHooks.h"
0026 #include "Pythia8/MergingHooks.h"
0027 #include "Pythia8/SimpleWeakShowerMEs.h"
0028 #include "Pythia8/DireBasics.h"
0029 #include "Pythia8/DireSplittingLibrary.h"
0030 #include "Pythia8/DireWeightContainer.h"
0031
0032 namespace Pythia8 {
0033
0034 //==========================================================================
0035
0036 // Data on radiating dipole ends, only used inside DireSpace.
0037
0038 class DireSpaceEnd {
0039
0040 public:
0041
0042   // Constructor.
0043   DireSpaceEnd( int systemIn = 0, int sideIn = 0, int iRadiatorIn = 0,
0044     int iRecoilerIn = 0, double pTmaxIn = 0., int colTypeIn = 0,
0045     int chgTypeIn = 0, int weakTypeIn = 0,  int MEtypeIn = 0,
0046     bool normalRecoilIn = true, int weakPolIn = 0,
0047     DireSingleColChain iSiblingsIn = DireSingleColChain(),
0048     vector<int> iSpectatorIn = vector<int>(),
0049     vector<double> massIn = vector<double>(),
0050     vector<int> allowedIn = vector<int>() ) :
0051     system(systemIn), side(sideIn), iRadiator(iRadiatorIn),
0052     iRecoiler(iRecoilerIn), pTmax(pTmaxIn), colType(colTypeIn),
0053     chgType(chgTypeIn), weakType(weakTypeIn), MEtype(MEtypeIn),
0054     normalRecoil(normalRecoilIn), weakPol(weakPolIn), nBranch(0),
0055     pT2Old(0.), zOld(0.5), mass(massIn), iSpectator(iSpectatorIn),
0056     allowedEmissions(allowedIn), iSiblings(iSiblingsIn) {
0057     idDaughter = idMother = idSister = iFinPol = 0;
0058     x1  = x2 = m2Dip = pT2 = z = xMo = Q2 = mSister = m2Sister = pT2corr
0059         = pT2Old = zOld = asymPol = sa1 = xa = pT2start = pT2stop = 0.;
0060     mRad = m2Rad = mRec = m2Rec = mDip = 0.;
0061     phi = phia1 = -1.;
0062   }
0063
0064   // Explicit copy constructor.
0065   DireSpaceEnd( const DireSpaceEnd& dip )
0066     : system(dip.system), side(dip.side), iRadiator(dip.iRadiator),
0067       iRecoiler(dip.iRecoiler), pTmax(dip.pTmax), colType(dip.colType),
0068       chgType(dip.chgType), weakType(dip.weakType), MEtype(dip.MEtype),
0069       normalRecoil(dip.normalRecoil), weakPol(dip.weakPol),
0070       nBranch(dip.nBranch), idDaughter(dip.idDaughter), idMother(dip.idMother),
0071       idSister(dip.idSister), iFinPol(dip.iFinPol), x1(dip.x1), x2(dip.x2),
0072       m2Dip(dip.m2Dip), pT2(dip.pT2), z(dip.z), xMo(dip.xMo), Q2(dip.Q2),
0073       mSister(dip.mSister), m2Sister(dip.m2Sister), pT2corr(dip.pT2corr),
0074       pT2Old(dip.pT2Old), zOld(dip.zOld), asymPol(dip.asymPol), phi(dip.phi),
0075       pT2start(dip.pT2start), pT2stop(dip.pT2stop),
0076       mRad(dip.mRad), m2Rad(dip.m2Rad), mRec(dip.mRec), m2Rec(dip.m2Rec),
0077       mDip(dip.mDip), sa1(dip.sa1), xa(dip.xa),
0078       phia1(dip.phia1), mass(dip.mass), iSpectator(dip.iSpectator),
0079       allowedEmissions(dip.allowedEmissions), iSiblings(dip.iSiblings) {}
0080
0081   // Assignment operator.
0082   DireSpaceEnd & operator=(const DireSpaceEnd &s) { if (this != &s)
0083     { system = s.system; side = s.side; iRadiator = s.iRadiator;
0084       iRecoiler = s.iRecoiler; pTmax = s.pTmax; colType = s.colType;
0085       chgType = s.chgType; weakType = s.weakType; MEtype = s.MEtype;
0086       normalRecoil = s.normalRecoil; weakPol = s.weakPol;
0087       nBranch = s.nBranch; idDaughter = s.idDaughter; idMother = s.idMother;
0088       idSister = s.idSister; iFinPol = s.iFinPol; x1 = s.x1; x2 = s.x2;
0089       m2Dip = s.m2Dip; pT2 = s.pT2; z = s.z; xMo = s.xMo; Q2 = s.Q2;
0090       mSister = s.mSister; m2Sister = s.m2Sister; pT2corr = s.pT2corr;
0091       pT2Old = s.pT2Old; zOld = s.zOld; asymPol = s.asymPol; phi = s.phi;
0092       pT2start = s.pT2start; pT2stop = s.pT2stop;
0093       mRad = s.mRad; m2Rad = s.m2Rad; mRec = s.mRec; m2Rec = s.m2Rec;
0094       mDip = s.mDip; sa1 = s.sa1; xa = s.xa; phia1 = s.phia1; mass = s.mass;
0095       iSpectator = s.iSpectator; allowedEmissions = s.allowedEmissions;
0096       iSiblings = s.iSiblings;} return *this; }
0097
0098   // Store values for trial emission.
0099   void store( int idDaughterIn, int idMotherIn, int idSisterIn,
0100     double x1In, double x2In, double m2DipIn, double pT2In, double zIn,
0101     double sa1In, double xaIn, double xMoIn, double Q2In, double mSisterIn,
0102     double m2SisterIn, double pT2corrIn, double phiIn = -1.,
0103     double phia1In = 1.) {
0104     idDaughter = idDaughterIn; idMother = idMotherIn;
0105     idSister = idSisterIn; x1 = x1In; x2 = x2In; m2Dip = m2DipIn;
0106     pT2 = pT2In; z = zIn; sa1 = sa1In; xa = xaIn; xMo = xMoIn; Q2 = Q2In;
0107     mSister = mSisterIn; m2Sister = m2SisterIn; pT2corr = pT2corrIn;
0108     mRad = m2Rad = mRec = m2Rec = mDip = 0.;
0109     phi = phiIn; phia1 = phia1In; }
0110
0111   // Basic properties related to evolution and matrix element corrections.
0112   int    system, side, iRadiator, iRecoiler;
0113   double pTmax;
0114   int    colType, chgType, weakType, MEtype;
0115   bool   normalRecoil;
0116   int    weakPol;
0117
0118   // Properties specific to current trial emission.
0119   int    nBranch, idDaughter, idMother, idSister, iFinPol;
0120   double x1, x2, m2Dip, pT2, z, xMo, Q2, mSister, m2Sister, pT2corr,
0121          pT2Old, zOld, asymPol, phi, pT2start, pT2stop,
0122          mRad, m2Rad, mRec, m2Rec, mDip;
0123
0124   // Properties of 1->3 splitting.
0125   double sa1, xa, phia1;
0126
0127   // Stored masses.
0128   vector<double> mass;
0129
0130   // Extended list of recoilers.
0131   vector<int> iSpectator;
0132
0133   vector<int> allowedEmissions;
0134
0135   // List of allowed emissions (to avoid double-counting, since one
0136   // particle can be part of many different dipoles.
0137   void appendAllowedEmt( int id) {
0138     if ( find(allowedEmissions.begin(), allowedEmissions.end(), id)
0139         == allowedEmissions.end() ) { allowedEmissions.push_back(id);}
0140   }
0141   void clearAllowedEmt() { allowedEmissions.resize(0); }
0142   bool canEmit() { return int(allowedEmissions.size() > 0); }
0143
0144   void init(const Event& state) {
0145     mRad   = state[iRadiator].m();
0146     mRec   = state[iRecoiler].m();
0147     mDip   = sqrt( abs(2. * state[iRadiator].p() * state[iRecoiler].p()));
0148     m2Rad  = pow2(mRad);
0149     m2Rec  = pow2(mRec);
0150     m2Dip  = pow2(mDip);
0151   }
0152
0153   void list() const {
0154     // Header.
0155     cout << "\n --------  DireSpaceEnd Listing  -------------- \n"
0156          << "\n    syst  side   rad   rec       pTmax  col  chg   ME rec \n"
0157          << fixed << setprecision(3);
0158     cout << setw(6) << system
0159          << setw(6) << side      << setw(6)  << iRadiator
0160          << setw(6) << iRecoiler << setw(12) << pTmax
0161          << setw(5) << colType   << setw(5)  << chgType
0162          << setw(5) << MEtype    << setw(4)
0163          << normalRecoil
0164          << setw(12) << m2Dip;
0165     for (int j = 0; j < int(allowedEmissions.size()); ++j)
0166       cout << setw(5) << allowedEmissions[j] << " ";
0167     cout << endl;
0168    // Done.
0169     cout << "\n --------  End DireSpaceEnd Listing  ------------"
0170          << "-------------------------------------------------------" << endl;
0171   }
0172
0173   DireSingleColChain iSiblings;
0174   void setSiblings(DireSingleColChain s) { clearSiblings(); iSiblings = s; }
0175   void clearSiblings() { iSiblings.clear(); }
0176
0177 };
0178
0179 //==========================================================================
0180
0181 // The DireSpace class does spacelike showers.
0182
0183 class DireSpace : public SpaceShower {
0184
0185 public:
0186
0187   // Constructor.
0188   DireSpace() {
0189     beamOffset          = 0;
0190     pTdampFudge         = 0.;
0191     mergingHooksPtr     = nullptr;
0192     splittingsPtr       = nullptr;
0193     weights             = 0;
0194     direInfoPtr         = nullptr;
0195     beamAPtr = beamBPtr = nullptr;
0196     printBanner         = true;
0197     nWeightsSave        = 0;
0198     isInitSave          = false;
0199     nMPI                = 0;
0200     usePDFalphas        = false;
0201     usePDF              = true;
0202     useSystems          = true;
0203     suppressLargeMECs   = false;
0204   }
0205
0206   DireSpace( MergingHooksPtr mergingHooksPtrIn, PartonVertexPtr ) :
0207       pTdampFudge(0.), mc(0.), mb(0.), m2c(0.), m2b(0.), m2cPhys(0.),
0208       m2bPhys(0.), renormMultFac(0.), factorMultFac(0.), fixedFacScale2(0.),
0209       alphaSvalue(0.), alphaS2pi(0.), Lambda3flav(0.), Lambda4flav(0.),
0210       Lambda5flav(0.), Lambda3flav2(0.), Lambda4flav2(0.), Lambda5flav2(0.),
0211       pT0Ref(0.), ecmRef(0.), ecmPow(0.), pTmin(0.), sCM(0.), eCM(0.), pT0(0.),
0212       pT20(0.), pT2min(0.), m2min(0.), mTolErr(0.), pTmaxFudgeMPI(0.),
0213       strengthIntAsym(0.), pT2minVariations(0.), pT2minEnhance(0.),
0214       pT2minMECs(0.), Q2minMECs(0.),
0215       alphaS2piOverestimate(0.), usePDFalphas(false), usePDFmasses(false),
0216       useSummedPDF(false), usePDF(true), useSystems(true),
0217       useGlobalMapIF(false), forceMassiveMap(false), useMassiveBeams(false),
0218       suppressLargeMECs(false) {
0219       mergingHooksPtr   = mergingHooksPtrIn;
0220       beamOffset        = 0;
0221       pTdampFudge       = 0.;
0222       splittingsPtr     = nullptr;
0223       weights           = 0;
0224       direInfoPtr       = nullptr;
0225       printBanner       = true;
0226       nWeightsSave      = 0;
0227       isInitSave        = false;
0228       nMPI = 0;
0229       beamAPtr = 0;
0230       beamBPtr = 0;
0231   }
0232
0233   // Destructor.
0234   virtual ~DireSpace() {}
0235
0236   // Initialize generation. Possibility to force re-initialization by hand.
0237   virtual void init(BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn);
0238
0239   bool initSplits() {
0240     if (splittingsPtr) splits = splittingsPtr->getSplittings();
0241     return (splits.size() > 0);
0242   }
0243
0244   // Initialize various pointers.
0245   // (Separated from rest of init since not virtual.)
0246   void reinitPtr(Info* infoPtrIn,  MergingHooksPtr mergingHooksPtrIn,
0247     DireSplittingLibrary* splittingsPtrIn, DireInfo* direInfoPtrIn) {
0248        infoPtr          = infoPtrIn;
0249        settingsPtr      = infoPtr->settingsPtr;
0250        particleDataPtr  = infoPtr->particleDataPtr;
0251        rndmPtr          = infoPtr->rndmPtr;
0252        partonSystemsPtr = infoPtr->partonSystemsPtr;
0253        userHooksPtr     = infoPtr->userHooksPtr;
0254        mergingHooksPtr  = mergingHooksPtrIn;
0255        splittingsPtr    = splittingsPtrIn;
0256        direInfoPtr      = direInfoPtrIn;
0257   }
0258
0259   void initVariations();
0260
0261   // Reset parton shower.
0262   void clear();
0263
0264   void setWeightContainerPtr(DireWeightContainer* weightsIn) {
0265     weights = weightsIn;}
0266
0267   // Find whether to limit maximum scale of emissions, and whether to dampen.
0268   virtual bool limitPTmax( Event& event, double Q2Fac = 0.,
0269     double Q2Ren = 0.);
0270
0271   // Potential enhancement factor of pTmax scale for hardest emission.
0272   virtual double enhancePTmax() const {return pTmaxFudge;}
0273
0274   // Prepare system for evolution; identify ME.
0275   void resetWeights();
0276   virtual void prepare( int iSys, Event& event, bool limitPTmaxIn = true);
0277
0278   // Update dipole list after each FSR emission.
0279   // Usage: update( iSys, event).
0280   virtual void update( int , Event&, bool = false);
0281
0282   // Update dipole list after initial-initial splitting.
0283   void updateAfterII( int iSysSelNow, int sideNow, int iDipSelNow,
0284     int eventSizeOldNow, int systemSizeOldNow, Event& event, int iDaughter,
0285     int iMother, int iSister, int iNewRecoiler, double pT2, double xNew);
0286
0287   // Update dipole list after initial-initial splitting.
0288   void updateAfterIF( int iSysSelNow, int sideNow, int iDipSelNow,
0289     int eventSizeOldNow, int systemSizeOldNow, Event& event, int iDaughter,
0290     int iRecoiler, int iMother, int iSister, int iNewRecoiler, int iNewOther,
0291     double pT2, double xNew);
0292
0293   // Select next pT in downwards evolution.
0294   virtual double pTnext( Event& event, double pTbegAll, double pTendAll,
0295     int nRadIn = -1, bool = false);
0296   double newPoint( const Event& event);
0297
0298   // Select next pT in downwards evolution, based only on dipole mass and
0299   // incoming momentum fraction.
0300   double pTnext( vector<DireSpaceEnd> dipEnds, Event event, double pTbegAll,
0301     double pTendAll, double m2dip, int type, double s = -1., double x = -1.);
0302   double noEmissionProbability( double pTbegAll, double pTendAll, double m2dip,
0303     int id, int type, double s = -1., double x = -1.);
0304
0305   // Setup branching kinematics.
0306   virtual bool branch( Event& event);
0307
0308   bool branch_II( Event& event, bool = false,
0309     DireSplitInfo* split = nullptr);
0310   bool branch_IF( Event& event, bool = false,
0311     DireSplitInfo* split = nullptr);
0312
0313   // Setup clustering kinematics.
0314   virtual Event clustered( const Event& state, int iRad, int iEmt, int iRecAft,
0315     string name) {
0316     pair <Event, pair<int,int> > reclus
0317       = clustered_internal(state, iRad, iEmt, iRecAft, name);
0318     if (reclus.first.size() > 0)
0319       reclus.first[0].mothers(reclus.second.first,reclus.second.second);
0320     return reclus.first;
0321   }
0322   pair <Event, pair<int,int> > clustered_internal( const Event& state,
0323     int iRad, int iEmt, int iRecAft, string name);
0324   bool cluster_II( const Event& state, int iRad,
0325     int iEmt, int iRecAft, int idRadBef, Particle& radBef, Particle& recBef,
0326     Event& partialState);
0327   bool cluster_IF( const Event& state, int iRad,
0328     int iEmt, int iRecAft, int idRadBef, Particle& radBef, Particle& recBef,
0329     Event& partialState);
0330
0331   // Return ordering variable.
0332   // From Pythia version 8.215 onwards no longer virtual.
0333   double pT2Space ( const Particle& rad, const Particle& emt,
0334     const Particle& rec) {
0335     if (rec.isFinal()) return pT2_IF(rad,emt,rec);
0336     return pT2_II(rad,emt,rec);
0337   }
0338
0339   double pT2_II ( const Particle& rad, const Particle& emt,
0340     const Particle& rec);
0341   double pT2_IF ( const Particle& rad, const Particle& emt,
0342     const Particle& rec);
0343
0344   // Return auxiliary variable.
0345   // From Pythia version 8.215 onwards no longer virtual.
0346   double zSpace ( const Particle& rad, const Particle& emt,
0347     const Particle& rec) {
0348     if (rec.isFinal()) return z_IF(rad,emt,rec);
0349     return z_II(rad,emt,rec);
0350   }
0351
0352   double z_II ( const Particle& rad, const Particle& emt,
0353     const Particle& rec);
0354   double z_IF ( const Particle& rad, const Particle& emt,
0355     const Particle& rec);
0356
0357   double m2dipSpace ( const Particle& rad, const Particle& emt,
0358     const Particle& rec) {
0359     if (rec.isFinal()) return m2dip_IF(rad,emt,rec);
0360     return m2dip_II(rad,emt,rec);
0361   }
0362   double m2dip_II ( const Particle& rad, const Particle& emt,
0363     const Particle& rec);
0364   double m2dip_IF ( const Particle& rad, const Particle& emt,
0365     const Particle& rec);
0366
0367   // From Pythia version 8.218 onwards.
0368   // Return the evolution variable.
0369   // Usage: getStateVariables( const Event& event,  int iRad, int iEmt,
0370   //                   int iRec, string name)
0371   // Important note:
0372   // - This map must contain an entry for the shower evolution variable,
0373   //   specified with key "t".
0374   // - This map must contain an entry for the shower evolution variable from
0375   //   which the shower would be restarted after a branching. This entry
0376   //   must have key "tRS",
0377   // - This map must contain an entry for the argument of \alpha_s used
0378   //   for the branching. This entry must have key "scaleAS".
0379   // - This map must contain an entry for the argument of the PDFs used
0380   //   for the branching. This entry must have key "scalePDF".
0381   virtual map<string, double> getStateVariables (const Event& state,
0382     int rad, int emt, int rec, string name);
0383
0384   // From Pythia version 8.215 onwards.
0385   // Check if attempted clustering is handled by timelike shower
0386   // Usage: isSpacelike( const Event& event,  int iRad, int iEmt,
0387   //                   int iRec, string name)
0388   virtual bool isSpacelike(const Event& state, int iRad, int, int, string)
0389     { return !state[iRad].isFinal(); }
0390
0391   // From Pythia version 8.215 onwards.
0392   // Return a string identifier of a splitting.
0393   // Usage: getSplittingName( const Event& event, int iRad, int iEmt, int iRec)
0394   virtual vector<string> getSplittingName( const Event& state, int iRad,
0395     int iEmt,int) { return splittingsPtr->getSplittingName(state,iRad,iEmt); }
0396
0397   // From Pythia version 8.215 onwards.
0398   // Return the splitting probability.
0399   // Usage: getSplittingProb( const Event& event, int iRad, int iEmt, int iRec)
0400   virtual double getSplittingProb( const Event& state, int iRad,
0401     int iEmt, int iRecAft, string);
0402
0403   virtual bool allowedSplitting( const Event& state, int iRad, int iEmt);
0404
0405   virtual vector<int> getRecoilers( const Event& state, int iRad, int iEmt,
0406     string name);
0407
0408   virtual double getCoupling( double mu2Ren, string name) {
0409     if (splits.find(name) != splits.end())
0410       return splits[name]->coupling(-1.,mu2Ren, 0., 1.);
0411     return 1.;
0412   }
0413
0414   bool isSymmetric( string name, const Particle* rad, const Particle* emt) {
0415     if (splits.find(name) != splits.end())
0416       return splits[name]->isSymmetric(rad,emt);
0417     return false;
0418   }
0419
0420   // Auxiliary function to return the position of a particle.
0421   // Should go int Event class eventually!
0422   int FindParticle( const Particle& particle, const Event& event,
0423     bool checkStatus = true );
0424
0425   // Print dipole list; for debug mainly.
0426   virtual void list() const;
0427
0428   Event makeHardEvent( int iSys, const Event& state, bool isProcess = false );
0429
0430   // Check that particle has sensible momentum.
0431   bool validMomentum( const Vec4& p, int id, int status);
0432
0433   // Check colour/flavour correctness of state.
0434   bool validEvent( const Event& state, bool isProcess = false );
0435
0436   // Check that mother-daughter-relations are correctly set.
0437   bool validMotherDaughter( const Event& state );
0438
0439   // Find index colour partner for input colour.
0440   int FindCol(int col, vector<int> iExc, const Event& event, int type,
0441     int iSys = -1);
0442
0443   // Pointers to the two incoming beams.
0444   BeamParticle*  getBeamA () { return beamAPtr; }
0445   BeamParticle*  getBeamB () { return beamBPtr; }
0446
0447   // Pointer to Standard Model couplings.
0448   CoupSM* getCoupSM () { return coupSMPtr; }
0449
0450   // Function to calculate the correct alphaS/2*Pi value, including
0451   // renormalisation scale variations + threshold matching.
0452   double alphasNow( double pT2, double renormMultFacNow = 1., int iSys = 0 );
0453
0454   bool isInit() { return isInitSave; }
0455
0456   // Function to calculate the absolute phase-sace boundary for emissions.
0457   double m2Max (int iDip, const Event& state) {
0458     if ( state[dipEnd[iDip].iRecoiler].isFinal()
0459       && state[dipEnd[iDip].iRadiator].isFinal())
0460       return dipEnd[iDip].m2Dip;
0461     int iSys = dipEnd[iDip].system;
0462     int inA = partonSystemsPtr->getInA(iSys);
0463     int inB = partonSystemsPtr->getInB(iSys);
0464     double x = 1.;
0465     int iRad(dipEnd[iDip].iRadiator), iRecNow(dipEnd[iDip].iRecoiler);
0466     if (hasPDF(state[iRad].id()) && iRad == inA)
0467       x *= state[inA].pPos() / state[0].m();
0468     if (hasPDF(state[iRad].id()) && iRad == inB)
0469       x *= state[inB].pNeg() / state[0].m();
0470     if (hasPDF(state[iRecNow].id()) && iRecNow == inA)
0471       x *= state[inA].pPos() / state[0].m();
0472     if (hasPDF(state[iRecNow].id()) && iRecNow == inB)
0473       x *= state[inB].pNeg() / state[0].m();
0474     return dipEnd[iDip].m2Dip/x;
0475   }
0476
0477   bool dryrun;
0478
0479 private:
0480
0481   friend class DireTimes;
0482
0483   // Number of times the same error message is repeated, unless overridden.
0484   static const int TIMESTOPRINT;
0485
0486   // Allow conversion from mb to pb.
0487   static const double CONVERTMB2PB;
0488
0489   // Colour factors.
0490   //static const double CA, CF, TR, NC;
0491   double CA, CF, TR, NC;
0492
0493   // Store common beam quantities.
0494   int    idASave, idBSave;
0495
0496 protected:
0497
0498   // Store properties to be returned by methods.
0499   int    iSysSel;
0500   double pTmaxFudge;
0501
0502 private:
0503
0504   // Constants: could only be changed in the code itself.
0505   static const int    MAXLOOPTINYPDF;
0506   static const double MCMIN, MBMIN, CTHRESHOLD, BTHRESHOLD, EVALPDFSTEP,
0507          TINYPDF, TINYKERNELPDF, TINYPT2, HEAVYPT2EVOL, HEAVYXEVOL,
0508          EXTRASPACEQ, LAMBDA3MARGIN, PT2MINWARN, LEPTONXMIN, LEPTONXMAX,
0509          LEPTONPT2MIN, LEPTONFUDGE, HEADROOMQ2Q, HEADROOMQ2G,
0510          HEADROOMG2G, HEADROOMG2Q, TINYMASS,
0511          PT2_INCREASE_OVERESTIMATE, KERNEL_HEADROOM;
0512   static const double DPHI_II, DPHI_IF;
0513   static const double G2QQPDFPOW1, G2QQPDFPOW2;
0514
0515   // Initialization data, normally only set once.
0516   bool   isInitSave, doQCDshower, doQEDshowerByQ, doQEDshowerByL,
0517          useSamePTasMPI, doMEcorrections, doMEafterFirst, doPhiPolAsym,
0518          doPhiIntAsym, doRapidityOrder, useFixedFacScale, doSecondHard,
0519          canVetoEmission, hasUserHooks, alphaSuseCMW, printBanner, doTrialNow;
0520   int    pTmaxMatch, pTdampMatch, alphaSorder, alphaSnfmax,
0521          nQuarkIn, enhanceScreening, nFinalMax, nFinalMaxMECs, kernelOrder,
0522          kernelOrderMPI, nWeightsSave, nMPI, asScheme;
0523   double pTdampFudge, mc, mb, m2c, m2b, m2cPhys, m2bPhys, renormMultFac,
0524          factorMultFac, fixedFacScale2, alphaSvalue, alphaS2pi, Lambda3flav,
0525          Lambda4flav, Lambda5flav, Lambda3flav2, Lambda4flav2, Lambda5flav2,
0526          pT0Ref, ecmRef, ecmPow, pTmin, sCM, eCM, pT0, pT20,
0527          pT2min, m2min, mTolErr, pTmaxFudgeMPI, strengthIntAsym,
0528          pT2minVariations, pT2minEnhance, pT2minMECs, Q2minMECs;
0529   double alphaS2piOverestimate;
0530   bool  usePDFalphas, usePDFmasses, useSummedPDF,  usePDF, useSystems,
0531         useGlobalMapIF, forceMassiveMap, useMassiveBeams, suppressLargeMECs;
0532
0533   unordered_map<int,double> pT2cutSave;
0534   double pT2cut(int id) {
0535     if (pT2cutSave.find(id) != pT2cutSave.end()) return pT2cutSave[id];
0536     // Else return maximal value.
0537     double ret = 0.;
0538     for ( unordered_map<int,double>::iterator it = pT2cutSave.begin();
0539       it != pT2cutSave.end(); ++it ) ret = max(ret, it->second);
0540     return ret;
0541   }
0542   double pT2cutMax(DireSpaceEnd* dip) {
0543     double ret = 0.;
0544     for (int i=0; i < int(dip->allowedEmissions.size()); ++i)
0545       ret = max( ret, pT2cut(dip->allowedEmissions[i]));
0546     return ret;
0547   }
0548   double pT2cutMin(DireSpaceEnd* dip) {
0549     double ret = 1e15;
0550     for (int i=0; i < int(dip->allowedEmissions.size()); ++i)
0551       ret = min( ret, pT2cut(dip->allowedEmissions[i]));
0552     return ret;
0553   }
0554
0555   bool doDecaysAsShower;
0556
0557   // alphaStrong and alphaEM calculations.
0558   AlphaStrong alphaS;
0559
0560   // Some current values.
0561   bool   sideA, dopTlimit1, dopTlimit2, dopTdamp;
0562   int    iNow, iRec, idDaughter, nRad, idResFirst, idResSecond;
0563   double xDaughter, x1Now, x2Now, m2Dip, m2Rec, pT2damp, pTbegRef, pdfScale2;
0564
0565   // List of emissions in different sides in different systems:
0566   vector<int> nRadA,nRadB;
0567
0568   // All dipole ends
0569   vector<DireSpaceEnd> dipEnd;
0570
0571   // Pointers to the current and hardest (so far) dipole ends.
0572   int iDipNow, iSysNow;
0573   DireSpaceEnd* dipEndNow;
0574   DireSplitInfo splitInfoSel;
0575   DireSplitting* splittingSel;
0576   int iDipSel;
0577   DireSpaceEnd* dipEndSel;
0578   unordered_map<string,double> kernelSel, kernelNow;
0579   double auxSel, overSel, boostSel, auxNow, overNow, boostNow;
0580
0581   void setupQCDdip( int iSys, int side, int colTag, int colSign,
0582     const Event& event, int MEtype, bool limitPTmaxIn);
0583
0584   void getGenDip( int iSys, int side, const Event& event,
0585     bool limitPTmaxIn, vector<DireSpaceEnd>& dipEnds );
0586
0587   void getQCDdip( int iRad, int colTag, int colSign,
0588     const Event& event, vector<DireSpaceEnd>& dipEnds);
0589
0590   // Function to set up and append a new dipole.
0591   bool appendDipole( const Event& state, int sys, int side,
0592     int iRad, int iRecNow, double pTmax, int colType,
0593     int chgType, int weakType, int MEtype, bool normalRecoil,
0594     int weakPolIn, vector<int> iSpectatorIn, vector<double> massIn,
0595     vector<DireSpaceEnd>& dipEnds);
0596
0597   vector<int> sharedColor(const Particle& rad, const Particle& rec);
0598
0599   // Function to set up and append a new dipole.
0600   void saveSiblings(const Event& state, int iSys = -1);
0601   void updateDipoles(const Event& state, int iSys = -1);
0602   bool updateAllowedEmissions( const Event& state, DireSpaceEnd* dip);
0603   bool appendAllowedEmissions( const Event& state, DireSpaceEnd* dip);
0604
0605   // Flag for failure in branch(...) that will force a retry of parton level.
0606   bool doRestart() const {return false;}
0607   // Tell if latest scattering was a gamma->qqbar.
0608   bool wasGamma2qqbar() { return false; }
0609   // Tell whether ISR has done a weak emission.
0610   bool getHasWeaklyRadiated() {return false;}
0611   int system() const { return iSysSel;}
0612
0613   // Evolve a QCD dipole end.
0614   void pT2nextQCD( double pT2begDip, double pT2endDip,
0615     DireSpaceEnd& dip, Event& event, double pT2endForce = -1.,
0616     double pT2freeze = 0., bool forceBranching = false);
0617   bool pT2nextQCD_II( double pT2begDip, double pT2endDip,
0618     DireSpaceEnd& dip, Event& event, double pT2endForce = -1.,
0619     double pT2freeze = 0., bool forceBranching = false);
0620   bool pT2nextQCD_IF( double pT2begDip, double pT2endDip,
0621     DireSpaceEnd& dip, Event& event, double pT2endForce = -1.,
0622     double pT2freeze = 0., bool forceBranching = false);
0623
0624   double tNextQCD( DireSpaceEnd*, double overestimateInt,
0625     double tOld, double tMin, double tFreeze=0., int algoType = 0);
0626   bool zCollNextQCD( DireSpaceEnd* dip, double zMin, double zMax,
0627     double tMin = 0., double tMax = 0.);
0628   bool virtNextQCD( DireSpaceEnd* dip, double tMin, double tMax,
0629     double zMin =-1., double zMax =-1.);
0630
0631   // Function to determine how often the integrated overestimate should be
0632   // recalculated.
0633   double evalpdfstep(int idRad, double pT2, double m2cp = -1.,
0634     double m2bp = -1.) {
0635     double ret = 0.2;
0636     if (m2cp < 0.) m2cp = particleDataPtr->m0(4);
0637     if (m2bp < 0.) m2bp = particleDataPtr->m0(5);
0638     // More steps close to the thresholds.
0639     if ( abs(idRad) == 4 && pT2 < 1.2*m2cp && pT2 > m2cp) ret = 1.0;
0640     if ( abs(idRad) == 5 && pT2 < 1.2*m2bp && pT2 > m2bp) ret = 1.0;
0641     return ret;
0642   }
0643
0644   double tinypdf( double x) {
0645     double xref = 0.01;
0646     return TINYPDF*log(1-x)/log(1-xref);
0647   }
0648
0649   // Function to check if id is part of the incoming hadron state.
0650   bool hasPDF (int id) {
0651     if ( !usePDF )                          return false;
0652     if ( particleDataPtr->colType(id) != 0) return true;
0653     if ( particleDataPtr->isLepton(id)
0654       && settingsPtr->flag("PDF:lepton"))   return true;
0655     return false;
0656   }
0657
0658   // Wrapper around PDF calls.
0659   double getXPDF( int id, double x, double t, int iSys = 0,
0660     BeamParticle* beam = nullptr, bool finalRec = false, double z = 0.,
0661     double m2dip = 0.) {
0662     // Return one if no PDF should be used.
0663     if (!hasPDF(id)) return 1.0;
0664     // Else get PDF from beam particle.
0665     BeamParticle* b = beam;
0666     if (b == nullptr) {
0667       if (beamAPtr != nullptr || beamBPtr != nullptr) {
0668         b = (beamAPtr != nullptr && particleDataPtr->isHadron(beamAPtr->id()))
0669             ? beamAPtr
0670           : (beamBPtr != nullptr && particleDataPtr->isHadron(beamBPtr->id()))
0671             ? beamBPtr : nullptr;
0672       }
0673       if (b == nullptr && beamAPtr != nullptr) b = beamAPtr;
0674       if (b == nullptr && beamBPtr != nullptr) b = beamBPtr;
0675     }
0676
0677     double scale2 = t;
0678     if (asScheme == 2 && z != 0) {
0679       if (!finalRec) {
0680         double xcs = (z * (1.-z) - t/m2dip) / (1.-z);
0681         double vcs = t/m2dip / (1.-z);
0682         double sab = m2dip/xcs;
0683         double saj = vcs*sab;
0684         double sjb = sab-saj-m2dip;
0685         scale2= abs(saj*sjb/sab);
0686       } else {
0687         double xcs = z;
0688         double ucs = t/m2dip / (1.-z);
0689         scale2 = (1-xcs)/xcs*ucs/(1-ucs)*m2dip;
0690       }
0691     }
0692
0693     double ret =  (useSummedPDF) ? b->xf(id, x, scale2)
0694                                  : b->xfISR(iSys,id, x, scale2);
0695     // Done.
0696     return ret;
0697   }
0698
0699   // Functions to extract beams w/o requiring parton systems pointer.
0700   int getInA ( int sys, const Event& state = Event() ) {
0701     if (useSystems) return partonSystemsPtr->getInA(sys);
0702     int inA = 0;
0703     for (int i=0; i < state.size(); ++i)
0704       if (state[i].mother1() == 1) {inA = i; break; }
0705     return inA;
0706   }
0707   int getInB ( int sys, const Event& state = Event() ) {
0708     if (useSystems) return partonSystemsPtr->getInB(sys);
0709     int inB = 0;
0710     for (int i=0; i < state.size(); ++i)
0711       if (state[i].mother1() == 2) {inB = i; break; }
0712     return inB;
0713   }
0714
0715
0716   DireSplittingLibrary* splittingsPtr;
0717
0718   // Number of proposed splittings in hard scattering systems.
0719   unordered_map<int,int> nProposedPT;
0720
0721   // Return headroom factors for integrated/differential overestimates.
0722   double overheadFactors( string, int, bool, double, double);
0723   double enhanceOverestimateFurther( string, int, double );
0724
0725   // Function to fill map of integrated overestimates.
0726   void getNewOverestimates( int, DireSpaceEnd*, const Event&, double,
0727     double, double, double, multimap<double,string>& );
0728
0729   // Function to fill map of integrated overestimates.
0730   double getPDFOverestimates( int, double, double, string, bool, double, int&,
0731     int&);
0732
0733   // Function to sum all integrated overestimates.
0734   void addNewOverestimates( multimap<double,string>, double&);
0735
0736   // Function to attach the correct alphaS weights to the kernels.
0737   void alphasReweight(double t, double talpha, int iSys, bool forceFixedAs,
0738     double& weight, double& fullWeight, double& overWeight,
0739     double renormMultFacNow);
0740
0741   // Function to evaluate the accept-probability, including picking of z.
0742   void getNewSplitting( const Event&, DireSpaceEnd*, double, double, double,
0743     double, double, int, string, bool, int&, int&, double&, double&,
0744     unordered_map<string,double>&, double&);
0745
0746   pair<bool, pair<double,double> > getMEC ( const Event& state,
0747     DireSplitInfo* splitInfo);
0748   bool applyMEC ( const Event& state, DireSplitInfo* splitInfo,
0749     vector<Event> auxEvent = vector<Event>() );
0750
0751   // Get particle masses.
0752   double getMass(int id, int strategy, double mass = 0.) {
0753     BeamParticle& beam = ( particleDataPtr->isHadron(beamAPtr->id()) )
0754                        ? *beamAPtr : *beamBPtr;
0755     bool usePDFmass = usePDFmasses
0756       && (toLower(settingsPtr->word("PDF:pSet")).find("lhapdf")
0757          != string::npos);
0758     double mRet = 0.;
0759     // Parton masses.
0760     if ( particleDataPtr->colType(id) != 0) {
0761       if (strategy == 1) mRet = particleDataPtr->m0(id);
0762       if (strategy == 2 &&  usePDFmass) mRet = beam.mQuarkPDF(id);
0763       if (strategy == 2 && !usePDFmass) mRet = particleDataPtr->m0(id);
0764       if (strategy == 3) mRet = mass;
0765       if (mRet < TINYMASS) mRet = 0.;
0766     // Masses of other particles.
0767     } else {
0768       mRet = particleDataPtr->m0(id);
0769       if (strategy == 3) mRet = mass;
0770       if (mRet < TINYMASS) mRet = 0.;
0771     }
0772     return pow2(max(0.,mRet));
0773   }
0774
0775   // Check if variables are in allowed phase space.
0776   bool inAllowedPhasespace(int kinType, double z, double pT2, double m2dip,
0777     double xOld, int splitType = 0, double m2RadBef = 0.,
0778     double m2r = 0.,  double m2s = 0., double m2e = 0.,
0779     vector<double> aux = vector<double>());
0780
0781   // Function to attach the correct alphaS weights to the kernels.
0782   // Auxiliary function to get number of flavours.
0783   double getNF(double pT2);
0784
0785   // Auxiliary functions to get beta function coefficients.
0786   double beta0 (double NF)
0787     { return 11./6.*CA - 2./3.*NF*TR; }
0788   double beta1 (double NF)
0789     { return 17./6.*pow2(CA) - (5./3.*CA+CF)*NF*TR; }
0790   double beta2 (double NF)
0791     { return 2857./432.*pow(CA,3)
0792     + (-1415./216.*pow2(CA) - 205./72.*CA*CF + pow2(CF)/4.) *TR*NF
0793     + ( 79.*CA + 66.*CF)/108.*pow2(TR*NF); }
0794
0795   // Identifier of the splitting
0796   string splittingNowName, splittingSelName;
0797
0798   // Weighted shower book-keeping.
0799   unordered_map<string, map<double,double> > acceptProbability;
0800   unordered_map<string, multimap<double,double> > rejectProbability;
0801
0802 public:
0803
0804   DireWeightContainer* weights;
0805   DireInfo* direInfoPtr;
0806   // List of splitting kernels.
0807   unordered_map<string, DireSplitting* > splits;
0808
0809 private:
0810
0811   bool doVariations;
0812
0813   // Dynamically adjustable overestimate factors.
0814   unordered_map<string, double > overhead;
0815   void scaleOverheadFactor(string name, double scale) {
0816     overhead[name] *= scale;
0817     return;
0818   }
0819   void resetOverheadFactors() {
0820     for ( unordered_map<string,double>::iterator it = overhead.begin();
0821       it != overhead.end(); ++it )
0822       it->second = 1.0;
0823     return;
0824   }
0825
0826   // Map to store some settings, to be passes to splitting kernels.
0827   unordered_map<string,bool> bool_settings;
0828
0829 };
0830
0831 //==========================================================================
0832
0833 } // end namespace Pythia8
0834
0835 #endif // Pythia8_DireSpace_H