PYTHIA  8.312
SimpleTimeShower.h
1 // SimpleTimeShower.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2024 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 
6 // Header file for the original simple timelike final-state showers.
7 // TimeDipoleEnd: data on a radiating dipole end in FSR.
8 // SimpleTimeShower: handles the showering description.
9 
10 #ifndef Pythia8_SimpleTimeShower_H
11 #define Pythia8_SimpleTimeShower_H
12 
13 #include "Pythia8/TimeShower.h"
14 #include "Pythia8/SimpleWeakShowerMEs.h"
15 
16 namespace Pythia8 {
17 
18 // Forward declaration of SplitOnia;
19 class SplitOnia;
20 
21 //==========================================================================
22 
23 // Data on radiating dipole ends; only used inside SimpleTimeShower class.
24 
26 
27 public:
28 
29  // Constructors.
30  TimeDipoleEnd() : iRadiator(-1), iRecoiler(-1), pTmax(0.), colType(0),
31  chgType(0), gamType(0), weakType(0), isrType(0), system(0), systemRec(0),
32  MEtype(0), iMEpartner(-1), weakPol(0), oniumType(0), isHiddenValley(false),
33  colvType(0), MEmix(0.), MEorder(true), MEsplit(true), MEgluinoRec(false),
34  isFlexible(false), hasJunction(false), flavour(), iAunt(),
35  mRad(), m2Rad(), mRec(), m2Rec(), mDip(), m2Dip(), m2DipCorr(), pT2(),
36  m2(), z(), mFlavour(), asymPol(), flexFactor(), pAccept() {}
37  TimeDipoleEnd(int iRadiatorIn, int iRecoilerIn, double pTmaxIn = 0.,
38  int colIn = 0, int chgIn = 0, int gamIn = 0, int weakTypeIn = 0,
39  int isrIn = 0, int systemIn = 0, int MEtypeIn = 0, int iMEpartnerIn = -1,
40  int weakPolIn = 0, int oniumTypeIn = 0,
41  bool isHiddenValleyIn = false, int colvTypeIn = 0, double MEmixIn = 0.,
42  bool MEorderIn = true, bool MEsplitIn = true, bool MEgluinoRecIn = false,
43  bool isFlexibleIn = false, bool hasJunctionIn = false) :
44  iRadiator(iRadiatorIn), iRecoiler(iRecoilerIn), pTmax(pTmaxIn),
45  colType(colIn), chgType(chgIn), gamType(gamIn), weakType(weakTypeIn),
46  isrType(isrIn), system(systemIn), systemRec(systemIn), MEtype(MEtypeIn),
47  iMEpartner(iMEpartnerIn), weakPol(weakPolIn), oniumType(oniumTypeIn),
48  isHiddenValley(isHiddenValleyIn), colvType(colvTypeIn), MEmix(MEmixIn),
49  MEorder (MEorderIn), MEsplit(MEsplitIn), MEgluinoRec(MEgluinoRecIn),
50  isFlexible(isFlexibleIn), hasJunction(hasJunctionIn),
51  flavour(), iAunt(), mRad(), m2Rad(), mRec(), m2Rec(), mDip(), m2Dip(),
52  m2DipCorr(), pT2(), m2(), z(), mFlavour(), asymPol(), flexFactor(),
53  pAccept() {}
54 
55  // Basic properties related to dipole and matrix element corrections.
56  int iRadiator, iRecoiler;
57  double pTmax;
58  int colType, chgType, gamType, weakType, isrType, system, systemRec,
59  MEtype, iMEpartner, weakPol, oniumType;
60  bool isHiddenValley;
61  int colvType;
62  double MEmix;
63  bool MEorder, MEsplit, MEgluinoRec, isFlexible;
64  bool hasJunction;
65 
66 
67  // Properties specific to current trial emission.
68  int flavour, iAunt;
69  double mRad, m2Rad, mRec, m2Rec, mDip, m2Dip, m2DipCorr,
70  pT2, m2, z, mFlavour, asymPol, flexFactor, pAccept;
71  double m2A{0}, m2B{0}, m2C{0}, m2gg{-1};
72 
73  // Pointer to an onium emission object, if present.
74  shared_ptr<SplitOnia> emissionPtr{};
75 
76 };
77 
78 //==========================================================================
79 
80 // The SimpleTimeShower class does timelike showers.
81 
82 class SimpleTimeShower : public TimeShower {
83 
84 public:
85 
86  // Constructor.
87  SimpleTimeShower() : hasWeaklyRadiated(), iSysSel(), pTmaxFudge(),
88  pTLastBranch(), doQCDshower(), doQEDshowerByQ(), doQEDshowerByL(),
89  doQEDshowerByOther(), doQEDshowerByGamma(), doWeakShower(),
90  doMEcorrections(), doMEextended(), doMEafterFirst(), doPhiPolAsym(),
91  doPhiPolAsymHard(), doInterleave(), doInterleaveResDec(),
92  allowBeamRecoil(), dampenBeamRecoil(), useFixedFacScale(),
93  allowRescatter(), canVetoEmission(), doHVshower(), brokenHVsym(),
94  setLambdaHV(), globalRecoil(), useLocalRecoilNow(), doSecondHard(),
95  hasUserHooks(), singleWeakEmission(), alphaSuseCMW(), vetoWeakJets(),
96  allowMPIdipole(), weakExternal(), recoilDeadCone(), doDipoleRecoil(),
97  doPartonVertex(), pTmaxMatch(), pTdampMatch(), alphaSorder(),
98  alphaSnfmax(), nGluonToQuark(), weightGluonToQuark(), recoilStrategyRF(),
99  alphaEMorder(), nGammaToQuark(), nGammaToLepton(), nCHV(), nFlavHV(),
100  idHV(), alphaHVorder(), nMaxGlobalRecoil(), weakMode(), pTdampFudge(),
101  mc(), mb(), m2c(), m2b(), renormMultFac(), factorMultFac(),
102  fixedFacScale2(), alphaSvalue(), alphaS2pi(), Lambda3flav(),
103  Lambda4flav(), Lambda5flav(), Lambda3flav2(), Lambda4flav2(),
104  Lambda5flav2(), scaleGluonToQuark(), extraGluonToQuark(), pTcolCutMin(),
105  pTcolCut(), pT2colCut(), pTchgQCut(), pT2chgQCut(), pTchgLCut(),
106  pT2chgLCut(), pTweakCut(), pT2weakCut(), mMaxGamma(), m2MaxGamma(),
107  mZ(), gammaZ(), thetaWRat(),
108  mW(), gammaW(), CFHV(), alphaHVfix(), alphaHVref(), LambdaHV(),
109  pThvCut(), pT2hvCut(), mHV(), pTmaxFudgeMPI(), weakEnhancement(),
110  vetoWeakDeltaR2(), twoHard(), dopTlimit1(), dopTlimit2(), dopTdamp(),
111  pT2damp(), kRad(), kEmt(), pdfScale2(), doTrialNow(), canEnhanceEmission(),
112  canEnhanceTrial(), canEnhanceET(), doUncertaintiesNow(), dipSel(),
113  iDipSel(), nHard(), nFinalBorn(), nMaxGlobalBranch(), nGlobal(),
114  globalRecoilMode(), limitMUQ(), weakHardSize() { beamOffset = 0;
115  pdfMode = 0; useSystems = true; }
116 
117  // Destructor.
118  virtual ~SimpleTimeShower() override {}
119 
120  // Initialize alphaStrong and related pTmin parameters.
121  virtual void init( BeamParticle* beamAPtrIn = 0,
122  BeamParticle* beamBPtrIn = 0) override;
123 
124  // Find whether to limit maximum scale of emissions, and whether to dampen.
125  virtual bool limitPTmax( Event& event, double Q2Fac = 0.,
126  double Q2Ren = 0.) override;
127 
128  // Top-level routine to do a full time-like shower in resonance decay.
129  virtual int shower( int iBeg, int iEnd, Event& event, double pTmax,
130  int nBranchMax = 0) override;
131 
132  // Top-level routine for QED radiation in hadronic decay to two leptons.
133  virtual int showerQED( int i1, int i2, Event& event, double pTmax) override;
134 
135  // Prepare process-level event for shower + interleaved resonance decays.
136  // Usage: prepareProcess( process, event, iPos).
137  // iPos provides mapping from process to event entries (before showering).
138  virtual void prepareProcess( Event& process, Event&, vector<int>&) override;
139 
140  // Global recoil: reset counters and store locations of outgoing partons.
141  virtual void prepareGlobal( Event& event) override;
142 
143  // Prepare system for evolution after each new interaction; identify ME.
144  virtual void prepare( int iSys, Event& event, bool limitPTmaxIn = true)
145  override;
146 
147  // Update dipole list after a multiparton interactions rescattering.
148  virtual void rescatterUpdate( int iSys, Event& event) override;
149 
150  // Update dipole list after each ISR emission.
151  virtual void update( int iSys, Event& event, bool hasWeakRad = false)
152  override;
153 
154  // Select next pT for FSR in downwards evolution.
155  virtual double pTnext( Event& event, double pTbegAll, double pTendAll,
156  bool isFirstTrial = false, bool doTrialIn = false) override;
157 
158  // Select next pT for interleaved resonance decays.
159  virtual double pTnextResDec() override;
160 
161  // ME corrections and kinematics that may give failure.
162  virtual bool branch( Event& event, bool isInterleaved = false) override;
163 
164  // Do a resonance decay + resonance shower (including any nested decays).
165  // May be called recursively for nested decays.
166  // Usage: resonanceShower( process, event, iPos, pTmerge), where iPos
167  // maps process to event entries, and pTmerge is the scale at which this
168  // system should be merged into its parent system.
169  virtual bool resonanceShower(Event& process, Event& event,
170  vector<int>& iPos, double qRestart) override;
171 
172  // Print dipole list; for debug mainly.
173  virtual void list() const override;
174 
175  // Initialize data members for calculation of uncertainty bands.
176  virtual bool initUncertainties() override;
177 
178  // Initialize data members for application of enhancements.
179  virtual bool initEnhancements() override;
180 
181  // Tell whether FSR has done a weak emission.
182  virtual bool getHasWeaklyRadiated() override {return hasWeaklyRadiated;}
183 
184  // Tell which system was the last processed one.
185  virtual int system() const override {return iSysSel;}
186 
187  // Potential enhancement factor of pTmax scale for hardest emission.
188  virtual double enhancePTmax() override {return pTmaxFudge;}
189 
190  // Provide the pT scale of the last branching in the above shower.
191  virtual double pTLastInShower() override {return pTLastBranch;}
192 
193  // Functions to directly extract the probability of no emission between two
194  // scales. These functions are not used in the Pythia core code, but can be
195  // used by external programs to interface with the shower directly.
196  double noEmissionProbability( double pTbegAll, double pTendAll, double m2dip,
197  int id, int type, double s = -1., double x = -1.) override;
198  double pTnext( vector<TimeDipoleEnd> dipEnds, Event event, double pTbegAll,
199  double pTendAll, double m2dip, int id, int type, double s = -1.,
200  double x = -1.);
201  int pdfMode;
202  bool useSystems;
203 
204 private:
205 
206  // Constants: could only be changed in the code itself.
207  static const double MCMIN, MBMIN, SIMPLIFYROOT, XMARGIN, XMARGINCOMB,
208  WTRATIOMAX,TINYPDF, LARGEM2, THRESHM2, LAMBDA3MARGIN1ORD,
209  LAMBDA3MARGIN2ORD,WEAKPSWEIGHT, WG2QEXTRA, REJECTFACTOR,
210  PROBLIMIT;
211  static const int NLOOPMAX;
212  // Rescatter: try to fix up recoil between systems
213  static const bool FIXRESCATTER, VETONEGENERGY;
214  static const double MAXVIRTUALITYFRACTION, MAXNEGENERGYFRACTION;
215 
216  // Store properties to be returned by methods.
217  bool hasWeaklyRadiated;
218  int iSysSel;
219  double pTmaxFudge, pTLastBranch;
220 
221  // Initialization data, normally only set once.
222  bool doQCDshower, doQEDshowerByQ, doQEDshowerByL, doQEDshowerByOther,
223  doQEDshowerByGamma, doWeakShower, doMEcorrections, doMEextended,
224  doMEafterFirst, doPhiPolAsym, doPhiPolAsymHard, doInterleave,
225  doInterleaveResDec, allowBeamRecoil, dampenBeamRecoil,
226  useFixedFacScale, allowRescatter, canVetoEmission, doHVshower,
227  brokenHVsym, setLambdaHV, globalRecoil, useLocalRecoilNow,
228  doSecondHard, hasUserHooks, singleWeakEmission, alphaSuseCMW,
229  vetoWeakJets, allowMPIdipole, weakExternal, recoilDeadCone,
230  doDipoleRecoil, doPartonVertex;
231  int pdfModeSave;
232  int pTmaxMatch, pTdampMatch, alphaSorder, alphaSnfmax, nGluonToQuark,
233  weightGluonToQuark, recoilStrategyRF, alphaEMorder, nGammaToQuark,
234  nGammaToLepton, nCHV, nFlavHV, idHV, alphaHVorder, nMaxGlobalRecoil,
235  weakMode;
236  double pTdampFudge, mc, mb, m2c, m2b, renormMultFac, factorMultFac,
237  fixedFacScale2, alphaSvalue, alphaS2pi, Lambda3flav, Lambda4flav,
238  Lambda5flav, Lambda3flav2, Lambda4flav2, Lambda5flav2,
239  scaleGluonToQuark, extraGluonToQuark, pTcolCutMin, pTcolCut,
240  pT2colCut, pTchgQCut, pT2chgQCut, pTchgLCut, pT2chgLCut,
241  pTweakCut, pT2weakCut, mMaxGamma, m2MaxGamma,
242  mZ, gammaZ, thetaWRat, mW, gammaW, CFHV,
243  alphaHVfix, alphaHVref, LambdaHV, pThvCut, pT2hvCut, mHV,
244  pTmaxFudgeMPI, weakEnhancement, vetoWeakDeltaR2;
245 
246  // alphaStrong, alphaEM and alpha_HV calculations.
247  AlphaStrong alphaS;
248  AlphaEM alphaEM;
249  AlphaSUN alphaHV;
250 
251  // Weak matrix elements used for corrections both of ISR and FSR.
252  SimpleWeakShowerMEs simpleWeakShowerMEs;
253 
254  // Some current values.
255  bool twoHard, dopTlimit1, dopTlimit2, dopTdamp;
256  double pT2damp, kRad, kEmt, pdfScale2;
257 
258  // Bookkeeping of enhanced actual or trial emissions (see EPJC (2013) 73).
259  bool doTrialNow, canEnhanceEmission, canEnhanceTrial, canEnhanceET,
260  doUncertaintiesNow;
261  string splittingNameNow, splittingNameSel;
262  map< double, pair<string,double> > enhanceFactors;
263  void storeEnhanceFactor(double pT2, string name, double enhanceFactorIn)
264  { enhanceFactors.insert(make_pair(pT2,make_pair(name,enhanceFactorIn)));}
265 
266  // All dipole ends and a pointer to the selected hardest dipole end.
267  vector<TimeDipoleEnd> dipEnd;
268  TimeDipoleEnd* dipSel;
269  int iDipSel;
270 
271  // Setup a dipole end, either QCD, QED/photon, weak or Hidden Valley one.
272  void setupQCDdip( int iSys, int i, int colTag, int colSign, Event& event,
273  bool isOctetOnium = false, bool limitPTmaxIn = true);
274  void setupQEDdip( int iSys, int i, int chgType, int gamType, Event& event,
275  bool limitPTmaxIn = true);
276  void setupWeakdip( int iSys, int i,int weakType, Event& event,
277  bool limitPTmaxIn = true);
278  // Special setup for weak dipoles if already specified in info ptr.
279  void setupWeakdipExternal(Event& event, bool limitPTmaxIn = true);
280  void setupHVdip( int iSys, int i, int colvType, Event& event,
281  bool limitPTmaxIn = true);
282 
283  // Apply ME corrections for a specific subsystem.
284  bool applyMECorrections(const Event& event, TimeDipoleEnd* dipBranch,
285  int iSysBranch);
286 
287  // Special setup for onium.
288  void regenerateOniumDipoles(Event & event);
289 
290  // Evolve a QCD dipole end.
291  void pT2nextQCD( double pT2begDip, double pT2sel, TimeDipoleEnd& dip,
292  Event& event);
293 
294  // Evolve a QED dipole end, either charged or photon.
295  void pT2nextQED( double pT2begDip, double pT2sel, TimeDipoleEnd& dip,
296  Event& event);
297 
298  // Evolve a weak dipole end.
299  void pT2nextWeak( double pT2begDip, double pT2sel, TimeDipoleEnd& dip,
300  Event& event);
301 
302  // Evolve a Hidden Valley dipole end.
303  void pT2nextHV( double pT2begDip, double pT2sel, TimeDipoleEnd& dip,
304  Event& );
305 
306  // Evolve an onium dipole end.
307  void pT2nextOnium( double pT2begDip, double pT2sel, TimeDipoleEnd& dip,
308  Event& );
309 
310  // Find kind of QCD ME correction.
311  void findMEtype( Event& event, TimeDipoleEnd& dip);
312 
313  // Find type of particle; used by findMEtype.
314  int findMEparticle( int id, bool isHiddenColour = false);
315 
316  // Find mixture of V and A in gamma/Z: energy- and flavour-dependent.
317  double gammaZmix( Event& event, int iRes, int iDau1, int iDau2);
318 
319  // Set up to calculate QCD ME correction with calcMEcorr.
320  double findMEcorr(TimeDipoleEnd* dip, Particle& rad, Particle& partner,
321  Particle& emt, bool cutEdge = true);
322 
323  // Set up to calculate weak ME corrections for t-channel processes.
324  double findMEcorrWeak(TimeDipoleEnd* dip, Vec4 rad, Vec4 rec,
325  Vec4 emt, Vec4 p3, Vec4 p4, Vec4 radBef, Vec4 recBef);
326 
327  // Calculate value of QCD ME correction.
328  double calcMEcorr( int kind, int combiIn, double mixIn, double x1,
329  double x2, double r1, double r2, double r3 = 0., bool cutEdge = true);
330 
331  // Find coefficient of azimuthal asymmetry from gluon polarization.
332  void findAsymPol( Event& event, TimeDipoleEnd* dip);
333 
334  // Compute scale for interleaved resonance decays
335  double calcPTresDec(Particle& res);
336 
337  // Rescatter: propagate dipole recoil to internal lines connecting systems.
338  bool rescatterPropagateRecoil( Event& event, Vec4& pNew);
339 
340  // Properties stored for (some) global recoil schemes.
341  // Vectors of event indices defining the hard process.
342  vector<int> hardPartons;
343  // Number of partons in current hard event, number of partons in Born-type
344  // hard event (to distinguish between S and H), maximally allowed number of
345  // global recoil branchings.
346  int nHard, nFinalBorn, nMaxGlobalBranch;
347  // Number of proposed splittings in hard scattering systems.
348  map<int,int> nProposed;
349  // Number of splittings with global recoil (currently only 1).
350  int nGlobal, globalRecoilMode;
351  // Switch to constrain recoiling system.
352  bool limitMUQ;
353 
354  // Calculate uncertainty-band weights for accepted/rejected trial branching.
355  void calcUncertainties( bool accept, double pAccept,
356  double enhance, TimeDipoleEnd* dip, Particle* radPtr,
357  Particle* emtPtr, Particle* recPtr );
358 
359  // 2 -> 2 information needed for the external weak setup.
360  vector<Vec4> weakMomenta;
361  vector<int> weak2to2lines;
362  int weakHardSize;
363 
364  // Settings and member variables for interleaved resonance decays.
365  bool doFSRinResonances{};
366  int resDecScaleChoice{-1}, iHardResDecSav{}, nRecurseResDec{};
367  vector<int> idResDecSav;
368  vector<double> pTresDecSav;
369 
370  // Settings and containers for control of MECs.
371  bool skipFirstMECinHardProc;
372  vector<int> skipFirstMECinResDecIDs{};
373 
374  // Onium emissions.
375  bool doOniumShower{false};
376  vector<SplitOniaPtr> oniumEmissions;
377  set<double> oniumThresholds;
378 
379 };
380 
381 //==========================================================================
382 
383 } // end namespace Pythia8
384 
385 #endif // Pythia8_SimpleTimeShower_H
SimpleTimeShower()
Constructor.
Definition: SimpleTimeShower.h:87
Definition: StandardModel.h:23
virtual double enhancePTmax() override
Potential enhancement factor of pTmax scale for hardest emission.
Definition: SimpleTimeShower.h:188
The Event class holds all info on the generated event.
Definition: Event.h:453
The TimeShower class does timelike showers.
Definition: TimeShower.h:33
Definition: BeamParticle.h:133
TimeDipoleEnd()
Constructors.
Definition: SimpleTimeShower.h:30
int flavour
Properties specific to current trial emission.
Definition: SimpleTimeShower.h:68
Definition: StandardModel.h:106
shared_ptr< SplitOnia > emissionPtr
Pointer to an onium emission object, if present.
Definition: SimpleTimeShower.h:74
Definition: StandardModel.h:226
virtual bool getHasWeaklyRadiated() override
Tell whether FSR has done a weak emission.
Definition: SimpleTimeShower.h:182
Data on radiating dipole ends; only used inside SimpleTimeShower class.
Definition: SimpleTimeShower.h:25
virtual double pTLastInShower() override
Provide the pT scale of the last branching in the above shower.
Definition: SimpleTimeShower.h:191
Definition: Event.h:32
virtual ~SimpleTimeShower() override
Destructor.
Definition: SimpleTimeShower.h:118
Definition: SimpleWeakShowerMEs.h:26
virtual int system() const override
Tell which system was the last processed one.
Definition: SimpleTimeShower.h:185
int iRadiator
Basic properties related to dipole and matrix element corrections.
Definition: SimpleTimeShower.h:56
Header for classes to set beam momentum and interaction vertex spread.
Definition: Analysis.h:20
The SimpleTimeShower class does timelike showers.
Definition: SimpleTimeShower.h:82
Definition: Basics.h:32