include/Pythia8/SimpleTimeShower.h

0001 // SimpleTimeShower.h is a part of the PYTHIA event generator.
0002 // Copyright (C) 2024 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 original simple timelike final-state showers.
0007 // TimeDipoleEnd: data on a radiating dipole end in FSR.
0008 // SimpleTimeShower: handles the showering description.
0009
0010 #ifndef Pythia8_SimpleTimeShower_H
0011 #define Pythia8_SimpleTimeShower_H
0012
0013 #include "Pythia8/TimeShower.h"
0014 #include "Pythia8/SimpleWeakShowerMEs.h"
0015
0016 namespace Pythia8 {
0017
0018 // Forward declaration of SplitOnia;
0019 class SplitOnia;
0020
0021 //==========================================================================
0022
0023 // Data on radiating dipole ends; only used inside SimpleTimeShower class.
0024
0025 class TimeDipoleEnd {
0026
0027 public:
0028
0029   // Constructors.
0030   TimeDipoleEnd() : iRadiator(-1), iRecoiler(-1), pTmax(0.), colType(0),
0031     chgType(0), gamType(0), weakType(0), isrType(0), system(0), systemRec(0),
0032     MEtype(0), iMEpartner(-1), weakPol(0), oniumType(0), isHiddenValley(false),
0033     colvType(0), MEmix(0.), MEorder(true), MEsplit(true), MEgluinoRec(false),
0034     isFlexible(false), hasJunction(false), flavour(), iAunt(),
0035     mRad(), m2Rad(), mRec(), m2Rec(), mDip(), m2Dip(), m2DipCorr(), pT2(),
0036     m2(), z(), mFlavour(), asymPol(), flexFactor(), pAccept() {}
0037   TimeDipoleEnd(int iRadiatorIn, int iRecoilerIn, double pTmaxIn = 0.,
0038     int colIn = 0, int chgIn = 0, int gamIn = 0, int weakTypeIn = 0,
0039     int isrIn = 0, int systemIn = 0, int MEtypeIn = 0, int iMEpartnerIn = -1,
0040     int weakPolIn = 0, int oniumTypeIn = 0,
0041     bool isHiddenValleyIn = false, int colvTypeIn = 0, double MEmixIn = 0.,
0042     bool MEorderIn = true, bool MEsplitIn = true, bool MEgluinoRecIn = false,
0043     bool isFlexibleIn = false, bool hasJunctionIn = false) :
0044     iRadiator(iRadiatorIn), iRecoiler(iRecoilerIn), pTmax(pTmaxIn),
0045     colType(colIn), chgType(chgIn), gamType(gamIn), weakType(weakTypeIn),
0046     isrType(isrIn), system(systemIn), systemRec(systemIn), MEtype(MEtypeIn),
0047     iMEpartner(iMEpartnerIn), weakPol(weakPolIn), oniumType(oniumTypeIn),
0048     isHiddenValley(isHiddenValleyIn), colvType(colvTypeIn), MEmix(MEmixIn),
0049     MEorder (MEorderIn), MEsplit(MEsplitIn), MEgluinoRec(MEgluinoRecIn),
0050     isFlexible(isFlexibleIn), hasJunction(hasJunctionIn),
0051     flavour(), iAunt(), mRad(), m2Rad(), mRec(), m2Rec(), mDip(), m2Dip(),
0052     m2DipCorr(), pT2(), m2(), z(), mFlavour(), asymPol(), flexFactor(),
0053     pAccept() {}
0054
0055   // Basic properties related to dipole and matrix element corrections.
0056   int    iRadiator, iRecoiler;
0057   double pTmax;
0058   int    colType, chgType, gamType, weakType, isrType, system, systemRec,
0059          MEtype, iMEpartner, weakPol, oniumType;
0060   bool   isHiddenValley;
0061   int    colvType;
0062   double MEmix;
0063   bool   MEorder, MEsplit, MEgluinoRec, isFlexible;
0064   bool   hasJunction;
0065
0066
0067   // Properties specific to current trial emission.
0068   int    flavour, iAunt;
0069   double mRad, m2Rad, mRec, m2Rec, mDip, m2Dip, m2DipCorr,
0070          pT2, m2, z, mFlavour, asymPol, flexFactor, pAccept;
0071   double m2A{0}, m2B{0}, m2C{0}, m2gg{-1};
0072
0073   // Pointer to an onium emission object, if present.
0074   shared_ptr<SplitOnia> emissionPtr{};
0075
0076 };
0077
0078 //==========================================================================
0079
0080 // The SimpleTimeShower class does timelike showers.
0081
0082 class SimpleTimeShower : public TimeShower {
0083
0084 public:
0085
0086   // Constructor.
0087   SimpleTimeShower() : hasWeaklyRadiated(), iSysSel(), pTmaxFudge(),
0088     pTLastBranch(), doQCDshower(), doQEDshowerByQ(), doQEDshowerByL(),
0089     doQEDshowerByOther(), doQEDshowerByGamma(), doWeakShower(),
0090     doMEcorrections(), doMEextended(), doMEafterFirst(), doPhiPolAsym(),
0091     doPhiPolAsymHard(), doInterleave(), doInterleaveResDec(),
0092     allowBeamRecoil(), dampenBeamRecoil(), useFixedFacScale(),
0093     allowRescatter(), canVetoEmission(), doHVshower(), brokenHVsym(),
0094     setLambdaHV(), globalRecoil(), useLocalRecoilNow(), doSecondHard(),
0095     hasUserHooks(), singleWeakEmission(), alphaSuseCMW(), vetoWeakJets(),
0096     allowMPIdipole(), weakExternal(), recoilDeadCone(), doDipoleRecoil(),
0097     doPartonVertex(), pTmaxMatch(), pTdampMatch(), alphaSorder(),
0098     alphaSnfmax(), nGluonToQuark(), weightGluonToQuark(), recoilStrategyRF(),
0099     alphaEMorder(), nGammaToQuark(), nGammaToLepton(), nCHV(), nFlavHV(),
0100     idHV(), alphaHVorder(), nMaxGlobalRecoil(), weakMode(), pTdampFudge(),
0101     mc(), mb(), m2c(), m2b(), renormMultFac(), factorMultFac(),
0102     fixedFacScale2(), alphaSvalue(), alphaS2pi(), Lambda3flav(),
0103     Lambda4flav(), Lambda5flav(), Lambda3flav2(), Lambda4flav2(),
0104     Lambda5flav2(), scaleGluonToQuark(), extraGluonToQuark(), pTcolCutMin(),
0105     pTcolCut(), pT2colCut(), pTchgQCut(), pT2chgQCut(), pTchgLCut(),
0106     pT2chgLCut(), pTweakCut(), pT2weakCut(), mMaxGamma(), m2MaxGamma(),
0107     mZ(), gammaZ(), thetaWRat(),
0108     mW(), gammaW(), CFHV(), alphaHVfix(), alphaHVref(), LambdaHV(),
0109     pThvCut(), pT2hvCut(), mHV(), pTmaxFudgeMPI(), weakEnhancement(),
0110     vetoWeakDeltaR2(), twoHard(), dopTlimit1(), dopTlimit2(), dopTdamp(),
0111     pT2damp(), kRad(), kEmt(), pdfScale2(), doTrialNow(), canEnhanceEmission(),
0112     canEnhanceTrial(), canEnhanceET(), doUncertaintiesNow(), dipSel(),
0113     iDipSel(), nHard(), nFinalBorn(), nMaxGlobalBranch(), nGlobal(),
0114     globalRecoilMode(), limitMUQ(), weakHardSize() { beamOffset = 0;
0115     pdfMode = 0; useSystems = true; }
0116
0117   // Destructor.
0118   virtual ~SimpleTimeShower() override {}
0119
0120   // Initialize alphaStrong and related pTmin parameters.
0121   virtual void init( BeamParticle* beamAPtrIn = 0,
0122     BeamParticle* beamBPtrIn = 0) override;
0123
0124   // Find whether to limit maximum scale of emissions, and whether to dampen.
0125   virtual bool limitPTmax( Event& event, double Q2Fac = 0.,
0126     double Q2Ren = 0.) override;
0127
0128   // Top-level routine to do a full time-like shower in resonance decay.
0129   virtual int shower( int iBeg, int iEnd, Event& event, double pTmax,
0130     int nBranchMax = 0) override;
0131
0132   // Top-level routine for QED radiation in hadronic decay to two leptons.
0133   virtual int showerQED( int i1, int i2, Event& event, double pTmax) override;
0134
0135   // Prepare process-level event for shower + interleaved resonance decays.
0136   // Usage: prepareProcess( process, event, iPos).
0137   // iPos provides mapping from process to event entries (before showering).
0138   virtual void prepareProcess( Event& process, Event&, vector<int>&) override;
0139
0140   // Global recoil: reset counters and store locations of outgoing partons.
0141   virtual void prepareGlobal( Event& event) override;
0142
0143   // Prepare system for evolution after each new interaction; identify ME.
0144   virtual void prepare( int iSys, Event& event, bool limitPTmaxIn = true)
0145     override;
0146
0147   // Update dipole list after a multiparton interactions rescattering.
0148   virtual void rescatterUpdate( int iSys, Event& event) override;
0149
0150   // Update dipole list after each ISR emission.
0151   virtual void update( int iSys, Event& event, bool hasWeakRad = false)
0152     override;
0153
0154   // Select next pT for FSR in downwards evolution.
0155   virtual double pTnext( Event& event, double pTbegAll, double pTendAll,
0156     bool isFirstTrial = false, bool doTrialIn = false) override;
0157
0158   // Select next pT for interleaved resonance decays.
0159   virtual double pTnextResDec() override;
0160
0161   // ME corrections and kinematics that may give failure.
0162   virtual bool branch( Event& event, bool isInterleaved = false) override;
0163
0164   // Do a resonance decay + resonance shower (including any nested decays).
0165   // May be called recursively for nested decays.
0166   // Usage: resonanceShower( process, event, iPos, pTmerge), where iPos
0167   // maps process to event entries, and pTmerge is the scale at which this
0168   // system should be merged into its parent system.
0169   virtual bool resonanceShower(Event& process, Event& event,
0170     vector<int>& iPos, double qRestart) override;
0171
0172   // Print dipole list; for debug mainly.
0173   virtual void list() const override;
0174
0175   // Initialize data members for calculation of uncertainty bands.
0176   virtual bool initUncertainties() override;
0177
0178   // Initialize data members for application of enhancements.
0179   virtual bool initEnhancements() override;
0180
0181   // Tell whether FSR has done a weak emission.
0182   virtual bool getHasWeaklyRadiated() override {return hasWeaklyRadiated;}
0183
0184   // Tell which system was the last processed one.
0185   virtual int system() const override {return iSysSel;}
0186
0187   // Potential enhancement factor of pTmax scale for hardest emission.
0188   virtual double enhancePTmax() override {return pTmaxFudge;}
0189
0190   // Provide the pT scale of the last branching in the above shower.
0191   virtual double pTLastInShower() override {return pTLastBranch;}
0192
0193   // Functions to directly extract the probability of no emission between two
0194   // scales. These functions are not used in the Pythia core code, but can be
0195   // used by external programs to interface with the shower directly.
0196   double noEmissionProbability( double pTbegAll, double pTendAll, double m2dip,
0197     int id, int type, double s = -1., double x = -1.) override;
0198   double pTnext( vector<TimeDipoleEnd> dipEnds, Event event, double pTbegAll,
0199     double pTendAll, double m2dip, int id, int type, double s = -1.,
0200     double x = -1.);
0201   int pdfMode;
0202   bool useSystems;
0203
0204 private:
0205
0206   // Constants: could only be changed in the code itself.
0207   static const double MCMIN, MBMIN, SIMPLIFYROOT, XMARGIN, XMARGINCOMB,
0208          WTRATIOMAX,TINYPDF, LARGEM2, THRESHM2, LAMBDA3MARGIN1ORD,
0209          LAMBDA3MARGIN2ORD,WEAKPSWEIGHT, WG2QEXTRA, REJECTFACTOR,
0210          PROBLIMIT;
0211   static const int NLOOPMAX;
0212   // Rescatter: try to fix up recoil between systems
0213   static const bool   FIXRESCATTER, VETONEGENERGY;
0214   static const double MAXVIRTUALITYFRACTION, MAXNEGENERGYFRACTION;
0215
0216   // Store properties to be returned by methods.
0217   bool   hasWeaklyRadiated;
0218   int    iSysSel;
0219   double pTmaxFudge, pTLastBranch;
0220
0221   // Initialization data, normally only set once.
0222   bool   doQCDshower, doQEDshowerByQ, doQEDshowerByL, doQEDshowerByOther,
0223          doQEDshowerByGamma, doWeakShower, doMEcorrections, doMEextended,
0224          doMEafterFirst, doPhiPolAsym, doPhiPolAsymHard, doInterleave,
0225          doInterleaveResDec, allowBeamRecoil, dampenBeamRecoil,
0226          useFixedFacScale, allowRescatter, canVetoEmission, doHVshower,
0227          brokenHVsym, setLambdaHV, globalRecoil, useLocalRecoilNow,
0228          doSecondHard, hasUserHooks, singleWeakEmission, alphaSuseCMW,
0229          vetoWeakJets, allowMPIdipole, weakExternal, recoilDeadCone,
0230          doDipoleRecoil, doPartonVertex;
0231   int    pdfModeSave;
0232   int    pTmaxMatch, pTdampMatch, alphaSorder, alphaSnfmax, nGluonToQuark,
0233          weightGluonToQuark, recoilStrategyRF, alphaEMorder, nGammaToQuark,
0234          nGammaToLepton, nCHV, nFlavHV, idHV, alphaHVorder, nMaxGlobalRecoil,
0235          weakMode;
0236   double pTdampFudge, mc, mb, m2c, m2b, renormMultFac, factorMultFac,
0237          fixedFacScale2, alphaSvalue, alphaS2pi, Lambda3flav, Lambda4flav,
0238          Lambda5flav, Lambda3flav2, Lambda4flav2, Lambda5flav2,
0239          scaleGluonToQuark, extraGluonToQuark, pTcolCutMin, pTcolCut,
0240          pT2colCut, pTchgQCut, pT2chgQCut, pTchgLCut, pT2chgLCut,
0241          pTweakCut, pT2weakCut, mMaxGamma, m2MaxGamma,
0242          mZ, gammaZ, thetaWRat, mW, gammaW, CFHV,
0243          alphaHVfix, alphaHVref, LambdaHV, pThvCut, pT2hvCut, mHV,
0244          pTmaxFudgeMPI, weakEnhancement, vetoWeakDeltaR2;
0245
0246   // alphaStrong, alphaEM and alpha_HV calculations.
0247   AlphaStrong alphaS;
0248   AlphaEM     alphaEM;
0249   AlphaSUN    alphaHV;
0250
0251   // Weak matrix elements used for corrections both of ISR and FSR.
0252   SimpleWeakShowerMEs  simpleWeakShowerMEs;
0253
0254   // Some current values.
0255   bool   twoHard, dopTlimit1, dopTlimit2, dopTdamp;
0256   double pT2damp, kRad, kEmt, pdfScale2;
0257
0258   // Bookkeeping of enhanced  actual or trial emissions (see EPJC (2013) 73).
0259   bool doTrialNow, canEnhanceEmission, canEnhanceTrial, canEnhanceET,
0260        doUncertaintiesNow;
0261   string splittingNameNow, splittingNameSel;
0262   map< double, pair<string,double> > enhanceFactors;
0263   void storeEnhanceFactor(double pT2, string name, double enhanceFactorIn)
0264     { enhanceFactors.insert(make_pair(pT2,make_pair(name,enhanceFactorIn)));}
0265
0266   // All dipole ends and a pointer to the selected hardest dipole end.
0267   vector<TimeDipoleEnd> dipEnd;
0268   TimeDipoleEnd* dipSel;
0269   int iDipSel;
0270
0271   // Setup a dipole end, either QCD, QED/photon, weak or Hidden Valley one.
0272   void setupQCDdip( int iSys, int i, int colTag,  int colSign, Event& event,
0273     bool isOctetOnium = false, bool limitPTmaxIn = true);
0274   void setupQEDdip( int iSys, int i, int chgType, int gamType, Event& event,
0275     bool limitPTmaxIn = true);
0276   void setupWeakdip( int iSys, int i,int weakType, Event& event,
0277     bool limitPTmaxIn = true);
0278   // Special setup for weak dipoles if already specified in info ptr.
0279   void setupWeakdipExternal(Event& event, bool limitPTmaxIn = true);
0280   void setupHVdip( int iSys, int i, int colvType, Event& event,
0281     bool limitPTmaxIn = true);
0282
0283   // Apply ME corrections for a specific subsystem.
0284   bool applyMECorrections(const Event& event, TimeDipoleEnd* dipBranch,
0285     int iSysBranch);
0286
0287   // Special setup for onium.
0288   void regenerateOniumDipoles(Event & event);
0289
0290   // Evolve a QCD dipole end.
0291   void pT2nextQCD( double pT2begDip, double pT2sel, TimeDipoleEnd& dip,
0292     Event& event);
0293
0294   // Evolve a QED dipole end, either charged or photon.
0295   void pT2nextQED( double pT2begDip, double pT2sel, TimeDipoleEnd& dip,
0296     Event& event);
0297
0298   // Evolve a weak dipole end.
0299   void pT2nextWeak( double pT2begDip, double pT2sel, TimeDipoleEnd& dip,
0300     Event& event);
0301
0302   // Evolve a Hidden Valley dipole end.
0303   void pT2nextHV( double pT2begDip, double pT2sel, TimeDipoleEnd& dip,
0304     Event& );
0305
0306   // Evolve an onium dipole end.
0307   void pT2nextOnium( double pT2begDip, double pT2sel, TimeDipoleEnd& dip,
0308     Event& );
0309
0310   // Find kind of QCD ME correction.
0311   void findMEtype( Event& event, TimeDipoleEnd& dip);
0312
0313   // Find type of particle; used by findMEtype.
0314   int findMEparticle( int id, bool isHiddenColour = false);
0315
0316   // Find mixture of V and A in gamma/Z: energy- and flavour-dependent.
0317   double gammaZmix( Event& event, int iRes, int iDau1, int iDau2);
0318
0319   // Set up to calculate QCD ME correction with calcMEcorr.
0320   double findMEcorr(TimeDipoleEnd* dip, Particle& rad, Particle& partner,
0321    Particle& emt, bool cutEdge = true);
0322
0323   // Set up to calculate weak ME corrections for t-channel processes.
0324   double findMEcorrWeak(TimeDipoleEnd* dip, Vec4 rad, Vec4 rec,
0325    Vec4 emt, Vec4 p3, Vec4 p4, Vec4 radBef, Vec4 recBef);
0326
0327   // Calculate value of QCD ME correction.
0328   double calcMEcorr( int kind, int combiIn, double mixIn, double x1,
0329     double x2, double r1, double r2, double r3 = 0., bool cutEdge = true);
0330
0331   // Find coefficient of azimuthal asymmetry from gluon polarization.
0332   void findAsymPol( Event& event, TimeDipoleEnd* dip);
0333
0334   // Compute scale for interleaved resonance decays
0335   double calcPTresDec(Particle& res);
0336
0337   // Rescatter: propagate dipole recoil to internal lines connecting systems.
0338   bool rescatterPropagateRecoil( Event& event, Vec4& pNew);
0339
0340   // Properties stored for (some) global recoil schemes.
0341   // Vectors of event indices defining the hard process.
0342   vector<int> hardPartons;
0343   // Number of partons in current hard event, number of partons in Born-type
0344   // hard event (to distinguish between S and H), maximally allowed number of
0345   // global recoil branchings.
0346   int nHard, nFinalBorn, nMaxGlobalBranch;
0347   // Number of proposed splittings in hard scattering systems.
0348   map<int,int> nProposed;
0349   // Number of splittings with global recoil (currently only 1).
0350   int nGlobal, globalRecoilMode;
0351   // Switch to constrain recoiling system.
0352   bool limitMUQ;
0353
0354   // Calculate uncertainty-band weights for accepted/rejected trial branching.
0355   void calcUncertainties( bool accept, double pAccept,
0356     double enhance, TimeDipoleEnd* dip, Particle* radPtr,
0357     Particle* emtPtr, Particle* recPtr );
0358
0359   // 2 -> 2 information needed for the external weak setup.
0360   vector<Vec4> weakMomenta;
0361   vector<int> weak2to2lines;
0362   int weakHardSize;
0363
0364   // Settings and member variables for interleaved resonance decays.
0365   bool doFSRinResonances{};
0366   int resDecScaleChoice{-1}, iHardResDecSav{}, nRecurseResDec{};
0367   vector<int> idResDecSav;
0368   vector<double> pTresDecSav;
0369
0370   // Settings and containers for control of MECs.
0371   bool skipFirstMECinHardProc;
0372   vector<int> skipFirstMECinResDecIDs{};
0373
0374   // Onium emissions.
0375   bool doOniumShower{false};
0376   vector<SplitOniaPtr> oniumEmissions;
0377   set<double> oniumThresholds;
0378
0379 };
0380
0381 //==========================================================================
0382
0383 } // end namespace Pythia8
0384
0385 #endif // Pythia8_SimpleTimeShower_H