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