PYTHIA  8.312
DireSpace.h
1 // DireSpace.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2024 Stefan Prestel, 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 spacelike initial-state showers.
7 // DireSpaceEnd: radiating dipole end in ISR.
8 // DireSpace: handles the showering description.
9 
10 #ifndef Pythia8_DireSpace_H
11 #define Pythia8_DireSpace_H
12 
13 #define DIRE_SPACE_VERSION "2.002"
14 
15 #include "Pythia8/Basics.h"
16 #include "Pythia8/SpaceShower.h"
17 #include "Pythia8/BeamParticle.h"
18 #include "Pythia8/Event.h"
19 #include "Pythia8/Info.h"
20 #include "Pythia8/ParticleData.h"
21 #include "Pythia8/PartonSystems.h"
22 #include "Pythia8/PythiaStdlib.h"
23 #include "Pythia8/Settings.h"
24 #include "Pythia8/StandardModel.h"
25 #include "Pythia8/UserHooks.h"
26 #include "Pythia8/MergingHooks.h"
27 #include "Pythia8/SimpleWeakShowerMEs.h"
28 #include "Pythia8/DireBasics.h"
29 #include "Pythia8/DireSplittingLibrary.h"
30 #include "Pythia8/DireWeightContainer.h"
31 
32 namespace Pythia8 {
33 
34 //==========================================================================
35 
36 // Data on radiating dipole ends, only used inside DireSpace.
37 
38 class DireSpaceEnd {
39 
40 public:
41 
42  // Constructor.
43  DireSpaceEnd( int systemIn = 0, int sideIn = 0, int iRadiatorIn = 0,
44  int iRecoilerIn = 0, double pTmaxIn = 0., int colTypeIn = 0,
45  int chgTypeIn = 0, int weakTypeIn = 0, int MEtypeIn = 0,
46  bool normalRecoilIn = true, int weakPolIn = 0,
47  DireSingleColChain iSiblingsIn = DireSingleColChain(),
48  vector<int> iSpectatorIn = vector<int>(),
49  vector<double> massIn = vector<double>(),
50  vector<int> allowedIn = vector<int>() ) :
51  system(systemIn), side(sideIn), iRadiator(iRadiatorIn),
52  iRecoiler(iRecoilerIn), pTmax(pTmaxIn), colType(colTypeIn),
53  chgType(chgTypeIn), weakType(weakTypeIn), MEtype(MEtypeIn),
54  normalRecoil(normalRecoilIn), weakPol(weakPolIn), nBranch(0),
55  pT2Old(0.), zOld(0.5), mass(massIn), iSpectator(iSpectatorIn),
56  allowedEmissions(allowedIn), iSiblings(iSiblingsIn) {
57  idDaughter = idMother = idSister = iFinPol = 0;
58  x1 = x2 = m2Dip = pT2 = z = xMo = Q2 = mSister = m2Sister = pT2corr
59  = pT2Old = zOld = asymPol = sa1 = xa = pT2start = pT2stop = 0.;
60  mRad = m2Rad = mRec = m2Rec = mDip = 0.;
61  phi = phia1 = -1.;
62  }
63 
64  // Explicit copy constructor.
65  DireSpaceEnd( const DireSpaceEnd& dip )
66  : system(dip.system), side(dip.side), iRadiator(dip.iRadiator),
67  iRecoiler(dip.iRecoiler), pTmax(dip.pTmax), colType(dip.colType),
68  chgType(dip.chgType), weakType(dip.weakType), MEtype(dip.MEtype),
69  normalRecoil(dip.normalRecoil), weakPol(dip.weakPol),
70  nBranch(dip.nBranch), idDaughter(dip.idDaughter), idMother(dip.idMother),
71  idSister(dip.idSister), iFinPol(dip.iFinPol), x1(dip.x1), x2(dip.x2),
72  m2Dip(dip.m2Dip), pT2(dip.pT2), z(dip.z), xMo(dip.xMo), Q2(dip.Q2),
73  mSister(dip.mSister), m2Sister(dip.m2Sister), pT2corr(dip.pT2corr),
74  pT2Old(dip.pT2Old), zOld(dip.zOld), asymPol(dip.asymPol), phi(dip.phi),
75  pT2start(dip.pT2start), pT2stop(dip.pT2stop),
76  mRad(dip.mRad), m2Rad(dip.m2Rad), mRec(dip.mRec), m2Rec(dip.m2Rec),
77  mDip(dip.mDip), sa1(dip.sa1), xa(dip.xa),
78  phia1(dip.phia1), mass(dip.mass), iSpectator(dip.iSpectator),
79  allowedEmissions(dip.allowedEmissions), iSiblings(dip.iSiblings) {}
80 
81  // Assignment operator.
82  DireSpaceEnd & operator=(const DireSpaceEnd &s) { if (this != &s)
83  { system = s.system; side = s.side; iRadiator = s.iRadiator;
84  iRecoiler = s.iRecoiler; pTmax = s.pTmax; colType = s.colType;
85  chgType = s.chgType; weakType = s.weakType; MEtype = s.MEtype;
86  normalRecoil = s.normalRecoil; weakPol = s.weakPol;
87  nBranch = s.nBranch; idDaughter = s.idDaughter; idMother = s.idMother;
88  idSister = s.idSister; iFinPol = s.iFinPol; x1 = s.x1; x2 = s.x2;
89  m2Dip = s.m2Dip; pT2 = s.pT2; z = s.z; xMo = s.xMo; Q2 = s.Q2;
90  mSister = s.mSister; m2Sister = s.m2Sister; pT2corr = s.pT2corr;
91  pT2Old = s.pT2Old; zOld = s.zOld; asymPol = s.asymPol; phi = s.phi;
92  pT2start = s.pT2start; pT2stop = s.pT2stop;
93  mRad = s.mRad; m2Rad = s.m2Rad; mRec = s.mRec; m2Rec = s.m2Rec;
94  mDip = s.mDip; sa1 = s.sa1; xa = s.xa; phia1 = s.phia1; mass = s.mass;
95  iSpectator = s.iSpectator; allowedEmissions = s.allowedEmissions;
96  iSiblings = s.iSiblings;} return *this; }
97 
98  // Store values for trial emission.
99  void store( int idDaughterIn, int idMotherIn, int idSisterIn,
100  double x1In, double x2In, double m2DipIn, double pT2In, double zIn,
101  double sa1In, double xaIn, double xMoIn, double Q2In, double mSisterIn,
102  double m2SisterIn, double pT2corrIn, double phiIn = -1.,
103  double phia1In = 1.) {
104  idDaughter = idDaughterIn; idMother = idMotherIn;
105  idSister = idSisterIn; x1 = x1In; x2 = x2In; m2Dip = m2DipIn;
106  pT2 = pT2In; z = zIn; sa1 = sa1In; xa = xaIn; xMo = xMoIn; Q2 = Q2In;
107  mSister = mSisterIn; m2Sister = m2SisterIn; pT2corr = pT2corrIn;
108  mRad = m2Rad = mRec = m2Rec = mDip = 0.;
109  phi = phiIn; phia1 = phia1In; }
110 
111  // Basic properties related to evolution and matrix element corrections.
112  int system, side, iRadiator, iRecoiler;
113  double pTmax;
114  int colType, chgType, weakType, MEtype;
115  bool normalRecoil;
116  int weakPol;
117 
118  // Properties specific to current trial emission.
119  int nBranch, idDaughter, idMother, idSister, iFinPol;
120  double x1, x2, m2Dip, pT2, z, xMo, Q2, mSister, m2Sister, pT2corr,
121  pT2Old, zOld, asymPol, phi, pT2start, pT2stop,
122  mRad, m2Rad, mRec, m2Rec, mDip;
123 
124  // Properties of 1->3 splitting.
125  double sa1, xa, phia1;
126 
127  // Stored masses.
128  vector<double> mass;
129 
130  // Extended list of recoilers.
131  vector<int> iSpectator;
132 
133  vector<int> allowedEmissions;
134 
135  // List of allowed emissions (to avoid double-counting, since one
136  // particle can be part of many different dipoles.
137  void appendAllowedEmt( int id) {
138  if ( find(allowedEmissions.begin(), allowedEmissions.end(), id)
139  == allowedEmissions.end() ) { allowedEmissions.push_back(id);}
140  }
141  void clearAllowedEmt() { allowedEmissions.resize(0); }
142  bool canEmit() { return int(allowedEmissions.size() > 0); }
143 
144  void init(const Event& state) {
145  mRad = state[iRadiator].m();
146  mRec = state[iRecoiler].m();
147  mDip = sqrt( abs(2. * state[iRadiator].p() * state[iRecoiler].p()));
148  m2Rad = pow2(mRad);
149  m2Rec = pow2(mRec);
150  m2Dip = pow2(mDip);
151  }
152 
153  void list() const {
154  // Header.
155  cout << "\n -------- DireSpaceEnd Listing -------------- \n"
156  << "\n syst side rad rec pTmax col chg ME rec \n"
157  << fixed << setprecision(3);
158  cout << setw(6) << system
159  << setw(6) << side << setw(6) << iRadiator
160  << setw(6) << iRecoiler << setw(12) << pTmax
161  << setw(5) << colType << setw(5) << chgType
162  << setw(5) << MEtype << setw(4)
163  << normalRecoil
164  << setw(12) << m2Dip;
165  for (int j = 0; j < int(allowedEmissions.size()); ++j)
166  cout << setw(5) << allowedEmissions[j] << " ";
167  cout << endl;
168  // Done.
169  cout << "\n -------- End DireSpaceEnd Listing ------------"
170  << "-------------------------------------------------------" << endl;
171  }
172 
173  DireSingleColChain iSiblings;
174  void setSiblings(DireSingleColChain s) { clearSiblings(); iSiblings = s; }
175  void clearSiblings() { iSiblings.clear(); }
176 
177 };
178 
179 //==========================================================================
180 
181 // The DireSpace class does spacelike showers.
182 
183 class DireSpace : public SpaceShower {
184 
185 public:
186 
187  // Constructor.
189  beamOffset = 0;
190  pTdampFudge = 0.;
191  mergingHooksPtr = nullptr;
192  splittingsPtr = nullptr;
193  weights = 0;
194  direInfoPtr = nullptr;
195  beamAPtr = beamBPtr = nullptr;
196  printBanner = true;
197  nWeightsSave = 0;
198  isInitSave = false;
199  nMPI = 0;
200  usePDFalphas = false;
201  usePDF = true;
202  useSystems = true;
203  suppressLargeMECs = false;
204  }
205 
206  DireSpace( MergingHooksPtr mergingHooksPtrIn, PartonVertexPtr ) :
207  pTdampFudge(0.), mc(0.), mb(0.), m2c(0.), m2b(0.), m2cPhys(0.),
208  m2bPhys(0.), renormMultFac(0.), factorMultFac(0.), fixedFacScale2(0.),
209  alphaSvalue(0.), alphaS2pi(0.), Lambda3flav(0.), Lambda4flav(0.),
210  Lambda5flav(0.), Lambda3flav2(0.), Lambda4flav2(0.), Lambda5flav2(0.),
211  pT0Ref(0.), ecmRef(0.), ecmPow(0.), pTmin(0.), sCM(0.), eCM(0.), pT0(0.),
212  pT20(0.), pT2min(0.), m2min(0.), mTolErr(0.), pTmaxFudgeMPI(0.),
213  strengthIntAsym(0.), pT2minVariations(0.), pT2minEnhance(0.),
214  pT2minMECs(0.), Q2minMECs(0.),
215  alphaS2piOverestimate(0.), usePDFalphas(false), usePDFmasses(false),
216  useSummedPDF(false), usePDF(true), useSystems(true),
217  useGlobalMapIF(false), forceMassiveMap(false), useMassiveBeams(false),
218  suppressLargeMECs(false) {
219  mergingHooksPtr = mergingHooksPtrIn;
220  beamOffset = 0;
221  pTdampFudge = 0.;
222  splittingsPtr = nullptr;
223  weights = 0;
224  direInfoPtr = nullptr;
225  printBanner = true;
226  nWeightsSave = 0;
227  isInitSave = false;
228  nMPI = 0;
229  beamAPtr = 0;
230  beamBPtr = 0;
231  }
232 
233  // Destructor.
234  virtual ~DireSpace() {}
235 
236  // Initialize generation. Possibility to force re-initialization by hand.
237  virtual void init(BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn);
238 
239  bool initSplits() {
240  if (splittingsPtr) splits = splittingsPtr->getSplittings();
241  return (splits.size() > 0);
242  }
243 
244  // Initialize various pointers.
245  // (Separated from rest of init since not virtual.)
246  void reinitPtr(Info* infoPtrIn, MergingHooksPtr mergingHooksPtrIn,
247  DireSplittingLibrary* splittingsPtrIn, DireInfo* direInfoPtrIn) {
248  infoPtr = infoPtrIn;
249  settingsPtr = infoPtr->settingsPtr;
250  particleDataPtr = infoPtr->particleDataPtr;
251  rndmPtr = infoPtr->rndmPtr;
252  partonSystemsPtr = infoPtr->partonSystemsPtr;
253  userHooksPtr = infoPtr->userHooksPtr;
254  mergingHooksPtr = mergingHooksPtrIn;
255  splittingsPtr = splittingsPtrIn;
256  direInfoPtr = direInfoPtrIn;
257  }
258 
259  void initVariations();
260 
261  // Reset parton shower.
262  void clear();
263 
264  void setWeightContainerPtr(DireWeightContainer* weightsIn) {
265  weights = weightsIn;}
266 
267  // Find whether to limit maximum scale of emissions, and whether to dampen.
268  virtual bool limitPTmax( Event& event, double Q2Fac = 0.,
269  double Q2Ren = 0.);
270 
271  // Potential enhancement factor of pTmax scale for hardest emission.
272  virtual double enhancePTmax() const {return pTmaxFudge;}
273 
274  // Prepare system for evolution; identify ME.
275  void resetWeights();
276  virtual void prepare( int iSys, Event& event, bool limitPTmaxIn = true);
277 
278  // Update dipole list after each FSR emission.
279  // Usage: update( iSys, event).
280  virtual void update( int , Event&, bool = false);
281 
282  // Update dipole list after initial-initial splitting.
283  void updateAfterII( int iSysSelNow, int sideNow, int iDipSelNow,
284  int eventSizeOldNow, int systemSizeOldNow, Event& event, int iDaughter,
285  int iMother, int iSister, int iNewRecoiler, double pT2, double xNew);
286 
287  // Update dipole list after initial-initial splitting.
288  void updateAfterIF( int iSysSelNow, int sideNow, int iDipSelNow,
289  int eventSizeOldNow, int systemSizeOldNow, Event& event, int iDaughter,
290  int iRecoiler, int iMother, int iSister, int iNewRecoiler, int iNewOther,
291  double pT2, double xNew);
292 
293  // Select next pT in downwards evolution.
294  virtual double pTnext( Event& event, double pTbegAll, double pTendAll,
295  int nRadIn = -1, bool = false);
296  double newPoint( const Event& event);
297 
298  // Select next pT in downwards evolution, based only on dipole mass and
299  // incoming momentum fraction.
300  double pTnext( vector<DireSpaceEnd> dipEnds, Event event, double pTbegAll,
301  double pTendAll, double m2dip, int type, double s = -1., double x = -1.);
302  double noEmissionProbability( double pTbegAll, double pTendAll, double m2dip,
303  int id, int type, double s = -1., double x = -1.);
304 
305  // Setup branching kinematics.
306  virtual bool branch( Event& event);
307 
308  bool branch_II( Event& event, bool = false,
309  DireSplitInfo* split = nullptr);
310  bool branch_IF( Event& event, bool = false,
311  DireSplitInfo* split = nullptr);
312 
313  // Setup clustering kinematics.
314  virtual Event clustered( const Event& state, int iRad, int iEmt, int iRecAft,
315  string name) {
316  pair <Event, pair<int,int> > reclus
317  = clustered_internal(state, iRad, iEmt, iRecAft, name);
318  if (reclus.first.size() > 0)
319  reclus.first[0].mothers(reclus.second.first,reclus.second.second);
320  return reclus.first;
321  }
322  pair <Event, pair<int,int> > clustered_internal( const Event& state,
323  int iRad, int iEmt, int iRecAft, string name);
324  bool cluster_II( const Event& state, int iRad,
325  int iEmt, int iRecAft, int idRadBef, Particle& radBef, Particle& recBef,
326  Event& partialState);
327  bool cluster_IF( const Event& state, int iRad,
328  int iEmt, int iRecAft, int idRadBef, Particle& radBef, Particle& recBef,
329  Event& partialState);
330 
331  // Return ordering variable.
332  // From Pythia version 8.215 onwards no longer virtual.
333  double pT2Space ( const Particle& rad, const Particle& emt,
334  const Particle& rec) {
335  if (rec.isFinal()) return pT2_IF(rad,emt,rec);
336  return pT2_II(rad,emt,rec);
337  }
338 
339  double pT2_II ( const Particle& rad, const Particle& emt,
340  const Particle& rec);
341  double pT2_IF ( const Particle& rad, const Particle& emt,
342  const Particle& rec);
343 
344  // Return auxiliary variable.
345  // From Pythia version 8.215 onwards no longer virtual.
346  double zSpace ( const Particle& rad, const Particle& emt,
347  const Particle& rec) {
348  if (rec.isFinal()) return z_IF(rad,emt,rec);
349  return z_II(rad,emt,rec);
350  }
351 
352  double z_II ( const Particle& rad, const Particle& emt,
353  const Particle& rec);
354  double z_IF ( const Particle& rad, const Particle& emt,
355  const Particle& rec);
356 
357  double m2dipSpace ( const Particle& rad, const Particle& emt,
358  const Particle& rec) {
359  if (rec.isFinal()) return m2dip_IF(rad,emt,rec);
360  return m2dip_II(rad,emt,rec);
361  }
362  double m2dip_II ( const Particle& rad, const Particle& emt,
363  const Particle& rec);
364  double m2dip_IF ( const Particle& rad, const Particle& emt,
365  const Particle& rec);
366 
367  // From Pythia version 8.218 onwards.
368  // Return the evolution variable.
369  // Usage: getStateVariables( const Event& event, int iRad, int iEmt,
370  // int iRec, string name)
371  // Important note:
372  // - This map must contain an entry for the shower evolution variable,
373  // specified with key "t".
374  // - This map must contain an entry for the shower evolution variable from
375  // which the shower would be restarted after a branching. This entry
376  // must have key "tRS",
377  // - This map must contain an entry for the argument of \alpha_s used
378  // for the branching. This entry must have key "scaleAS".
379  // - This map must contain an entry for the argument of the PDFs used
380  // for the branching. This entry must have key "scalePDF".
381  virtual map<string, double> getStateVariables (const Event& state,
382  int rad, int emt, int rec, string name);
383 
384  // From Pythia version 8.215 onwards.
385  // Check if attempted clustering is handled by timelike shower
386  // Usage: isSpacelike( const Event& event, int iRad, int iEmt,
387  // int iRec, string name)
388  virtual bool isSpacelike(const Event& state, int iRad, int, int, string)
389  { return !state[iRad].isFinal(); }
390 
391  // From Pythia version 8.215 onwards.
392  // Return a string identifier of a splitting.
393  // Usage: getSplittingName( const Event& event, int iRad, int iEmt, int iRec)
394  virtual vector<string> getSplittingName( const Event& state, int iRad,
395  int iEmt,int) { return splittingsPtr->getSplittingName(state,iRad,iEmt); }
396 
397  // From Pythia version 8.215 onwards.
398  // Return the splitting probability.
399  // Usage: getSplittingProb( const Event& event, int iRad, int iEmt, int iRec)
400  virtual double getSplittingProb( const Event& state, int iRad,
401  int iEmt, int iRecAft, string);
402 
403  virtual bool allowedSplitting( const Event& state, int iRad, int iEmt);
404 
405  virtual vector<int> getRecoilers( const Event& state, int iRad, int iEmt,
406  string name);
407 
408  virtual double getCoupling( double mu2Ren, string name) {
409  if (splits.find(name) != splits.end())
410  return splits[name]->coupling(-1.,mu2Ren, 0., 1.);
411  return 1.;
412  }
413 
414  bool isSymmetric( string name, const Particle* rad, const Particle* emt) {
415  if (splits.find(name) != splits.end())
416  return splits[name]->isSymmetric(rad,emt);
417  return false;
418  }
419 
420  // Auxiliary function to return the position of a particle.
421  // Should go int Event class eventually!
422  int FindParticle( const Particle& particle, const Event& event,
423  bool checkStatus = true );
424 
425  // Print dipole list; for debug mainly.
426  virtual void list() const;
427 
428  Event makeHardEvent( int iSys, const Event& state, bool isProcess = false );
429 
430  // Check that particle has sensible momentum.
431  bool validMomentum( const Vec4& p, int id, int status);
432 
433  // Check colour/flavour correctness of state.
434  bool validEvent( const Event& state, bool isProcess = false );
435 
436  // Check that mother-daughter-relations are correctly set.
437  bool validMotherDaughter( const Event& state );
438 
439  // Find index colour partner for input colour.
440  int FindCol(int col, vector<int> iExc, const Event& event, int type,
441  int iSys = -1);
442 
443  // Pointers to the two incoming beams.
444  BeamParticle* getBeamA () { return beamAPtr; }
445  BeamParticle* getBeamB () { return beamBPtr; }
446 
447  // Pointer to Standard Model couplings.
448  CoupSM* getCoupSM () { return coupSMPtr; }
449 
450  // Function to calculate the correct alphaS/2*Pi value, including
451  // renormalisation scale variations + threshold matching.
452  double alphasNow( double pT2, double renormMultFacNow = 1., int iSys = 0 );
453 
454  bool isInit() { return isInitSave; }
455 
456  // Function to calculate the absolute phase-sace boundary for emissions.
457  double m2Max (int iDip, const Event& state) {
458  if ( state[dipEnd[iDip].iRecoiler].isFinal()
459  && state[dipEnd[iDip].iRadiator].isFinal())
460  return dipEnd[iDip].m2Dip;
461  int iSys = dipEnd[iDip].system;
462  int inA = partonSystemsPtr->getInA(iSys);
463  int inB = partonSystemsPtr->getInB(iSys);
464  double x = 1.;
465  int iRad(dipEnd[iDip].iRadiator), iRecNow(dipEnd[iDip].iRecoiler);
466  if (hasPDF(state[iRad].id()) && iRad == inA)
467  x *= state[inA].pPos() / state[0].m();
468  if (hasPDF(state[iRad].id()) && iRad == inB)
469  x *= state[inB].pNeg() / state[0].m();
470  if (hasPDF(state[iRecNow].id()) && iRecNow == inA)
471  x *= state[inA].pPos() / state[0].m();
472  if (hasPDF(state[iRecNow].id()) && iRecNow == inB)
473  x *= state[inB].pNeg() / state[0].m();
474  return dipEnd[iDip].m2Dip/x;
475  }
476 
477  bool dryrun;
478 
479 private:
480 
481  friend class DireTimes;
482 
483  // Number of times the same error message is repeated, unless overridden.
484  static const int TIMESTOPRINT;
485 
486  // Allow conversion from mb to pb.
487  static const double CONVERTMB2PB;
488 
489  // Colour factors.
490  //static const double CA, CF, TR, NC;
491  double CA, CF, TR, NC;
492 
493  // Store common beam quantities.
494  int idASave, idBSave;
495 
496 protected:
497 
498  // Store properties to be returned by methods.
499  int iSysSel;
500  double pTmaxFudge;
501 
502 private:
503 
504  // Constants: could only be changed in the code itself.
505  static const int MAXLOOPTINYPDF;
506  static const double MCMIN, MBMIN, CTHRESHOLD, BTHRESHOLD, EVALPDFSTEP,
507  TINYPDF, TINYKERNELPDF, TINYPT2, HEAVYPT2EVOL, HEAVYXEVOL,
508  EXTRASPACEQ, LAMBDA3MARGIN, PT2MINWARN, LEPTONXMIN, LEPTONXMAX,
509  LEPTONPT2MIN, LEPTONFUDGE, HEADROOMQ2Q, HEADROOMQ2G,
510  HEADROOMG2G, HEADROOMG2Q, TINYMASS,
511  PT2_INCREASE_OVERESTIMATE, KERNEL_HEADROOM;
512  static const double DPHI_II, DPHI_IF;
513  static const double G2QQPDFPOW1, G2QQPDFPOW2;
514 
515  // Initialization data, normally only set once.
516  bool isInitSave, doQCDshower, doQEDshowerByQ, doQEDshowerByL,
517  useSamePTasMPI, doMEcorrections, doMEafterFirst, doPhiPolAsym,
518  doPhiIntAsym, doRapidityOrder, useFixedFacScale, doSecondHard,
519  canVetoEmission, hasUserHooks, alphaSuseCMW, printBanner, doTrialNow;
520  int pTmaxMatch, pTdampMatch, alphaSorder, alphaSnfmax,
521  nQuarkIn, enhanceScreening, nFinalMax, nFinalMaxMECs, kernelOrder,
522  kernelOrderMPI, nWeightsSave, nMPI, asScheme;
523  double pTdampFudge, mc, mb, m2c, m2b, m2cPhys, m2bPhys, renormMultFac,
524  factorMultFac, fixedFacScale2, alphaSvalue, alphaS2pi, Lambda3flav,
525  Lambda4flav, Lambda5flav, Lambda3flav2, Lambda4flav2, Lambda5flav2,
526  pT0Ref, ecmRef, ecmPow, pTmin, sCM, eCM, pT0, pT20,
527  pT2min, m2min, mTolErr, pTmaxFudgeMPI, strengthIntAsym,
528  pT2minVariations, pT2minEnhance, pT2minMECs, Q2minMECs;
529  double alphaS2piOverestimate;
530  bool usePDFalphas, usePDFmasses, useSummedPDF, usePDF, useSystems,
531  useGlobalMapIF, forceMassiveMap, useMassiveBeams, suppressLargeMECs;
532 
533  unordered_map<int,double> pT2cutSave;
534  double pT2cut(int id) {
535  if (pT2cutSave.find(id) != pT2cutSave.end()) return pT2cutSave[id];
536  // Else return maximal value.
537  double ret = 0.;
538  for ( unordered_map<int,double>::iterator it = pT2cutSave.begin();
539  it != pT2cutSave.end(); ++it ) ret = max(ret, it->second);
540  return ret;
541  }
542  double pT2cutMax(DireSpaceEnd* dip) {
543  double ret = 0.;
544  for (int i=0; i < int(dip->allowedEmissions.size()); ++i)
545  ret = max( ret, pT2cut(dip->allowedEmissions[i]));
546  return ret;
547  }
548  double pT2cutMin(DireSpaceEnd* dip) {
549  double ret = 1e15;
550  for (int i=0; i < int(dip->allowedEmissions.size()); ++i)
551  ret = min( ret, pT2cut(dip->allowedEmissions[i]));
552  return ret;
553  }
554 
555  bool doDecaysAsShower;
556 
557  // alphaStrong and alphaEM calculations.
558  AlphaStrong alphaS;
559 
560  // Some current values.
561  bool sideA, dopTlimit1, dopTlimit2, dopTdamp;
562  int iNow, iRec, idDaughter, nRad, idResFirst, idResSecond;
563  double xDaughter, x1Now, x2Now, m2Dip, m2Rec, pT2damp, pTbegRef, pdfScale2;
564 
565  // List of emissions in different sides in different systems:
566  vector<int> nRadA,nRadB;
567 
568  // All dipole ends
569  vector<DireSpaceEnd> dipEnd;
570 
571  // Pointers to the current and hardest (so far) dipole ends.
572  int iDipNow, iSysNow;
573  DireSpaceEnd* dipEndNow;
574  DireSplitInfo splitInfoSel;
575  DireSplitting* splittingSel;
576  int iDipSel;
577  DireSpaceEnd* dipEndSel;
578  unordered_map<string,double> kernelSel, kernelNow;
579  double auxSel, overSel, boostSel, auxNow, overNow, boostNow;
580 
581  void setupQCDdip( int iSys, int side, int colTag, int colSign,
582  const Event& event, int MEtype, bool limitPTmaxIn);
583 
584  void getGenDip( int iSys, int side, const Event& event,
585  bool limitPTmaxIn, vector<DireSpaceEnd>& dipEnds );
586 
587  void getQCDdip( int iRad, int colTag, int colSign,
588  const Event& event, vector<DireSpaceEnd>& dipEnds);
589 
590  // Function to set up and append a new dipole.
591  bool appendDipole( const Event& state, int sys, int side,
592  int iRad, int iRecNow, double pTmax, int colType,
593  int chgType, int weakType, int MEtype, bool normalRecoil,
594  int weakPolIn, vector<int> iSpectatorIn, vector<double> massIn,
595  vector<DireSpaceEnd>& dipEnds);
596 
597  vector<int> sharedColor(const Particle& rad, const Particle& rec);
598 
599  // Function to set up and append a new dipole.
600  void saveSiblings(const Event& state, int iSys = -1);
601  void updateDipoles(const Event& state, int iSys = -1);
602  bool updateAllowedEmissions( const Event& state, DireSpaceEnd* dip);
603  bool appendAllowedEmissions( const Event& state, DireSpaceEnd* dip);
604 
605  // Flag for failure in branch(...) that will force a retry of parton level.
606  bool doRestart() const {return false;}
607  // Tell if latest scattering was a gamma->qqbar.
608  bool wasGamma2qqbar() { return false; }
609  // Tell whether ISR has done a weak emission.
610  bool getHasWeaklyRadiated() {return false;}
611  int system() const { return iSysSel;}
612 
613  // Evolve a QCD dipole end.
614  void pT2nextQCD( double pT2begDip, double pT2endDip,
615  DireSpaceEnd& dip, Event& event, double pT2endForce = -1.,
616  double pT2freeze = 0., bool forceBranching = false);
617  bool pT2nextQCD_II( double pT2begDip, double pT2endDip,
618  DireSpaceEnd& dip, Event& event, double pT2endForce = -1.,
619  double pT2freeze = 0., bool forceBranching = false);
620  bool pT2nextQCD_IF( double pT2begDip, double pT2endDip,
621  DireSpaceEnd& dip, Event& event, double pT2endForce = -1.,
622  double pT2freeze = 0., bool forceBranching = false);
623 
624  double tNextQCD( DireSpaceEnd*, double overestimateInt,
625  double tOld, double tMin, double tFreeze=0., int algoType = 0);
626  bool zCollNextQCD( DireSpaceEnd* dip, double zMin, double zMax,
627  double tMin = 0., double tMax = 0.);
628  bool virtNextQCD( DireSpaceEnd* dip, double tMin, double tMax,
629  double zMin =-1., double zMax =-1.);
630 
631  // Function to determine how often the integrated overestimate should be
632  // recalculated.
633  double evalpdfstep(int idRad, double pT2, double m2cp = -1.,
634  double m2bp = -1.) {
635  double ret = 0.2;
636  if (m2cp < 0.) m2cp = particleDataPtr->m0(4);
637  if (m2bp < 0.) m2bp = particleDataPtr->m0(5);
638  // More steps close to the thresholds.
639  if ( abs(idRad) == 4 && pT2 < 1.2*m2cp && pT2 > m2cp) ret = 1.0;
640  if ( abs(idRad) == 5 && pT2 < 1.2*m2bp && pT2 > m2bp) ret = 1.0;
641  return ret;
642  }
643 
644  double tinypdf( double x) {
645  double xref = 0.01;
646  return TINYPDF*log(1-x)/log(1-xref);
647  }
648 
649  // Function to check if id is part of the incoming hadron state.
650  bool hasPDF (int id) {
651  if ( !usePDF ) return false;
652  if ( particleDataPtr->colType(id) != 0) return true;
653  if ( particleDataPtr->isLepton(id)
654  && settingsPtr->flag("PDF:lepton")) return true;
655  return false;
656  }
657 
658  // Wrapper around PDF calls.
659  double getXPDF( int id, double x, double t, int iSys = 0,
660  BeamParticle* beam = nullptr, bool finalRec = false, double z = 0.,
661  double m2dip = 0.) {
662  // Return one if no PDF should be used.
663  if (!hasPDF(id)) return 1.0;
664  // Else get PDF from beam particle.
665  BeamParticle* b = beam;
666  if (b == nullptr) {
667  if (beamAPtr != nullptr || beamBPtr != nullptr) {
668  b = (beamAPtr != nullptr && particleDataPtr->isHadron(beamAPtr->id()))
669  ? beamAPtr
670  : (beamBPtr != nullptr && particleDataPtr->isHadron(beamBPtr->id()))
671  ? beamBPtr : nullptr;
672  }
673  if (b == nullptr && beamAPtr != nullptr) b = beamAPtr;
674  if (b == nullptr && beamBPtr != nullptr) b = beamBPtr;
675  }
676 
677  double scale2 = t;
678  if (asScheme == 2 && z != 0) {
679  if (!finalRec) {
680  double xcs = (z * (1.-z) - t/m2dip) / (1.-z);
681  double vcs = t/m2dip / (1.-z);
682  double sab = m2dip/xcs;
683  double saj = vcs*sab;
684  double sjb = sab-saj-m2dip;
685  scale2= abs(saj*sjb/sab);
686  } else {
687  double xcs = z;
688  double ucs = t/m2dip / (1.-z);
689  scale2 = (1-xcs)/xcs*ucs/(1-ucs)*m2dip;
690  }
691  }
692 
693  double ret = (useSummedPDF) ? b->xf(id, x, scale2)
694  : b->xfISR(iSys,id, x, scale2);
695  // Done.
696  return ret;
697  }
698 
699  // Functions to extract beams w/o requiring parton systems pointer.
700  int getInA ( int sys, const Event& state = Event() ) {
701  if (useSystems) return partonSystemsPtr->getInA(sys);
702  int inA = 0;
703  for (int i=0; i < state.size(); ++i)
704  if (state[i].mother1() == 1) {inA = i; break; }
705  return inA;
706  }
707  int getInB ( int sys, const Event& state = Event() ) {
708  if (useSystems) return partonSystemsPtr->getInB(sys);
709  int inB = 0;
710  for (int i=0; i < state.size(); ++i)
711  if (state[i].mother1() == 2) {inB = i; break; }
712  return inB;
713  }
714 
715 
716  DireSplittingLibrary* splittingsPtr;
717 
718  // Number of proposed splittings in hard scattering systems.
719  unordered_map<int,int> nProposedPT;
720 
721  // Return headroom factors for integrated/differential overestimates.
722  double overheadFactors( string, int, bool, double, double);
723  double enhanceOverestimateFurther( string, int, double );
724 
725  // Function to fill map of integrated overestimates.
726  void getNewOverestimates( int, DireSpaceEnd*, const Event&, double,
727  double, double, double, multimap<double,string>& );
728 
729  // Function to fill map of integrated overestimates.
730  double getPDFOverestimates( int, double, double, string, bool, double, int&,
731  int&);
732 
733  // Function to sum all integrated overestimates.
734  void addNewOverestimates( multimap<double,string>, double&);
735 
736  // Function to attach the correct alphaS weights to the kernels.
737  void alphasReweight(double t, double talpha, int iSys, bool forceFixedAs,
738  double& weight, double& fullWeight, double& overWeight,
739  double renormMultFacNow);
740 
741  // Function to evaluate the accept-probability, including picking of z.
742  void getNewSplitting( const Event&, DireSpaceEnd*, double, double, double,
743  double, double, int, string, bool, int&, int&, double&, double&,
744  unordered_map<string,double>&, double&);
745 
746  pair<bool, pair<double,double> > getMEC ( const Event& state,
747  DireSplitInfo* splitInfo);
748  bool applyMEC ( const Event& state, DireSplitInfo* splitInfo,
749  vector<Event> auxEvent = vector<Event>() );
750 
751  // Get particle masses.
752  double getMass(int id, int strategy, double mass = 0.) {
753  BeamParticle& beam = ( particleDataPtr->isHadron(beamAPtr->id()) )
754  ? *beamAPtr : *beamBPtr;
755  bool usePDFmass = usePDFmasses
756  && (toLower(settingsPtr->word("PDF:pSet")).find("lhapdf")
757  != string::npos);
758  double mRet = 0.;
759  // Parton masses.
760  if ( particleDataPtr->colType(id) != 0) {
761  if (strategy == 1) mRet = particleDataPtr->m0(id);
762  if (strategy == 2 && usePDFmass) mRet = beam.mQuarkPDF(id);
763  if (strategy == 2 && !usePDFmass) mRet = particleDataPtr->m0(id);
764  if (strategy == 3) mRet = mass;
765  if (mRet < TINYMASS) mRet = 0.;
766  // Masses of other particles.
767  } else {
768  mRet = particleDataPtr->m0(id);
769  if (strategy == 3) mRet = mass;
770  if (mRet < TINYMASS) mRet = 0.;
771  }
772  return pow2(max(0.,mRet));
773  }
774 
775  // Check if variables are in allowed phase space.
776  bool inAllowedPhasespace(int kinType, double z, double pT2, double m2dip,
777  double xOld, int splitType = 0, double m2RadBef = 0.,
778  double m2r = 0., double m2s = 0., double m2e = 0.,
779  vector<double> aux = vector<double>());
780 
781  // Function to attach the correct alphaS weights to the kernels.
782  // Auxiliary function to get number of flavours.
783  double getNF(double pT2);
784 
785  // Auxiliary functions to get beta function coefficients.
786  double beta0 (double NF)
787  { return 11./6.*CA - 2./3.*NF*TR; }
788  double beta1 (double NF)
789  { return 17./6.*pow2(CA) - (5./3.*CA+CF)*NF*TR; }
790  double beta2 (double NF)
791  { return 2857./432.*pow(CA,3)
792  + (-1415./216.*pow2(CA) - 205./72.*CA*CF + pow2(CF)/4.) *TR*NF
793  + ( 79.*CA + 66.*CF)/108.*pow2(TR*NF); }
794 
795  // Identifier of the splitting
796  string splittingNowName, splittingSelName;
797 
798  // Weighted shower book-keeping.
799  unordered_map<string, map<double,double> > acceptProbability;
800  unordered_map<string, multimap<double,double> > rejectProbability;
801 
802 public:
803 
804  DireWeightContainer* weights;
805  DireInfo* direInfoPtr;
806  // List of splitting kernels.
807  unordered_map<string, DireSplitting* > splits;
808 
809 private:
810 
811  bool doVariations;
812 
813  // Dynamically adjustable overestimate factors.
814  unordered_map<string, double > overhead;
815  void scaleOverheadFactor(string name, double scale) {
816  overhead[name] *= scale;
817  return;
818  }
819  void resetOverheadFactors() {
820  for ( unordered_map<string,double>::iterator it = overhead.begin();
821  it != overhead.end(); ++it )
822  it->second = 1.0;
823  return;
824  }
825 
826  // Map to store some settings, to be passes to splitting kernels.
827  unordered_map<string,bool> bool_settings;
828 
829 };
830 
831 //==========================================================================
832 
833 } // end namespace Pythia8
834 
835 #endif // Pythia8_DireSpace_H
constexpr double pow2(const double &x)
Powers of small integers - for balance speed/code clarity.
Definition: PythiaStdlib.h:182
CoupSM * getCoupSM()
Pointer to Standard Model couplings.
Definition: DireSpace.h:448
virtual vector< string > getSplittingName(const Event &state, int iRad, int iEmt, int)
Definition: DireSpace.h:394
Definition: StandardModel.h:23
BeamParticle * getBeamA()
Pointers to the two incoming beams.
Definition: DireSpace.h:444
void appendAllowedEmt(int id)
Definition: DireSpace.h:137
Data on radiating dipole ends, only used inside DireSpace.
Definition: DireSpace.h:38
Definition: Info.h:45
The Event class holds all info on the generated event.
Definition: Event.h:453
Definition: BeamParticle.h:133
vector< int > iSpectator
Extended list of recoilers.
Definition: DireSpace.h:131
int system
Basic properties related to evolution and matrix element corrections.
Definition: DireSpace.h:112
virtual Event clustered(const Event &state, int iRad, int iEmt, int iRecAft, string name)
Setup clustering kinematics.
Definition: DireSpace.h:314
virtual double enhancePTmax() const
Potential enhancement factor of pTmax scale for hardest emission.
Definition: DireSpace.h:272
DireSpaceEnd(const DireSpaceEnd &dip)
Explicit copy constructor.
Definition: DireSpace.h:65
Definition: DireBasics.h:374
string toLower(const string &name, bool trim=true)
Definition: PythiaStdlib.cc:17
double pT2Space(const Particle &rad, const Particle &emt, const Particle &rec)
Definition: DireSpace.h:333
The DireSpace class does spacelike showers.
Definition: DireSpace.h:183
vector< double > mass
Stored masses.
Definition: DireSpace.h:128
DireSpace()
Constructor.
Definition: DireSpace.h:188
DireSpaceEnd(int systemIn=0, int sideIn=0, int iRadiatorIn=0, int iRecoilerIn=0, double pTmaxIn=0., int colTypeIn=0, int chgTypeIn=0, int weakTypeIn=0, int MEtypeIn=0, bool normalRecoilIn=true, int weakPolIn=0, DireSingleColChain iSiblingsIn=DireSingleColChain(), vector< int > iSpectatorIn=vector< int >(), vector< double > massIn=vector< double >(), vector< int > allowedIn=vector< int >())
Constructor.
Definition: DireSpace.h:43
void reinitPtr(Info *infoPtrIn, MergingHooksPtr mergingHooksPtrIn, DireSplittingLibrary *splittingsPtrIn, DireInfo *direInfoPtrIn)
Definition: DireSpace.h:246
int iSysSel
Store properties to be returned by methods.
Definition: DireSpace.h:499
Definition: DireSplittings.h:53
double sa1
Properties of 1->3 splitting.
Definition: DireSpace.h:125
Definition of color chains.
Definition: DireSplitInfo.h:28
void store(int idDaughterIn, int idMotherIn, int idSisterIn, double x1In, double x2In, double m2DipIn, double pT2In, double zIn, double sa1In, double xaIn, double xMoIn, double Q2In, double mSisterIn, double m2SisterIn, double pT2corrIn, double phiIn=-1., double phia1In=1.)
Store values for trial emission.
Definition: DireSpace.h:99
double xf(int idIn, double x, double Q2)
Standard parton distributions.
Definition: BeamParticle.h:240
unordered_map< string, DireSplitting * > splits
List of splitting kernels.
Definition: DireSpace.h:807
double mQuarkPDF(int idIn)
Return quark masses used in the PDF fit (LHAPDF6 only).
Definition: BeamParticle.h:270
Container for all shower weights, including handling.
Definition: DireWeightContainer.h:82
virtual ~DireSpace()
Destructor.
Definition: DireSpace.h:234
Definition: StandardModel.h:135
int size() const
Event record size.
Definition: Event.h:504
DireSpaceEnd & operator=(const DireSpaceEnd &s)
Assignment operator.
Definition: DireSpace.h:82
Definition: Event.h:32
int nBranch
Properties specific to current trial emission.
Definition: DireSpace.h:119
double m2Max(int iDip, const Event &state)
Function to calculate the absolute phase-sace boundary for emissions.
Definition: DireSpace.h:457
virtual bool isSpacelike(const Event &state, int iRad, int, int, string)
Definition: DireSpace.h:388
double zSpace(const Particle &rad, const Particle &emt, const Particle &rec)
Definition: DireSpace.h:346
Header for classes to set beam momentum and interaction vertex spread.
Definition: Analysis.h:20
Definition: DireSplitInfo.h:212
Definition: DireSplittingLibrary.h:35
The SpaceShower class does spacelike showers.
Definition: SpaceShower.h:33
const double CA
Colour factors in Vincia normalization.
Definition: VinciaCommon.h:40
bool validEvent(const Event &event)
The Merging class.
Definition: DireMerging.cc:24
The DireTimes class does timelike showers.
Definition: DireTimes.h:209
Settings * settingsPtr
Pointer to the settings database.
Definition: Info.h:80
Definition: Basics.h:32
void list() const
Definition: DireSpace.h:153