include/Pythia8/SimpleTimeShower.h

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