StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
DireTimes.h
1 // DireTimes.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2020 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 timelike final-state showers.
7 // DireTimesEnd: data on a radiating dipole end.
8 // DireTimes: handles the showering description.
9 
10 #ifndef Pythia8_DireTimes_H
11 #define Pythia8_DireTimes_H
12 
13 #define DIRE_TIMES_VERSION "2.002"
14 
15 #include "Pythia8/Basics.h"
16 #include "Pythia8/TimeShower.h"
17 #include "Pythia8/BeamParticle.h"
18 #include "Pythia8/ProcessLevel.h"
19 #include "Pythia8/Event.h"
20 #include "Pythia8/Info.h"
21 #include "Pythia8/ParticleData.h"
22 #include "Pythia8/PartonSystems.h"
23 #include "Pythia8/PythiaStdlib.h"
24 #include "Pythia8/Settings.h"
25 #include "Pythia8/StandardModel.h"
26 #include "Pythia8/UserHooks.h"
27 #include "Pythia8/MergingHooks.h"
28 
29 #include "Pythia8/DireBasics.h"
30 #include "Pythia8/DireSplittingLibrary.h"
31 #include "Pythia8/DireWeightContainer.h"
32 
33 namespace Pythia8 {
34 
35 //==========================================================================
36 
37 // Data on radiating dipole ends; only used inside DireTimes class.
38 
39 class DireTimesEnd {
40 
41 public:
42 
43  // Constructors.
44  DireTimesEnd() : iRadiator(-1), iRecoiler(-1), pTmax(0.), colType(0),
45  chgType(0), gamType(0), weakType(0), isrType(0), system(0), systemRec(0),
46  MEtype(0), iMEpartner(-1), weakPol(0), isOctetOnium(false),
47  isHiddenValley(false), colvType(0), MEmix(0.), MEorder(true),
48  MEsplit(true), MEgluinoRec(false), isFlexible(false), flavour(0), iAunt(0),
49  idRadAft(0), idEmtAft(0) {
50  mRad = m2Rad = mRec = m2Rec = mDip = m2Dip = m2DipCorr = pT2 = m2 = z
51  = mFlavour = asymPol = flexFactor = phi = 0.;
52  sa1 = xa = phia1 = pT2start = pT2stop = pT2Old = 0.;
53  }
54 
55  DireTimesEnd(int iRadiatorIn, int iRecoilerIn, double pTmaxIn = 0.,
56  int colIn = 0, int chgIn = 0, int gamIn = 0, int weakTypeIn = 0,
57  int isrIn = 0, int systemIn = 0, int MEtypeIn = 0, int iMEpartnerIn = -1,
58  int weakPolIn = 0, bool isOctetOniumIn = false,
59  DireSingleColChain iSiblingsIn = DireSingleColChain(),
60  bool isHiddenValleyIn = false,
61  int colvTypeIn = 0, double MEmixIn = 0.,
62  bool MEorderIn = true, bool MEsplitIn = true, bool MEgluinoRecIn = false,
63  bool isFlexibleIn = false, int idRadAftIn = 0, int idEmtAftIn = 0,
64  vector<int> iSpectatorIn = vector<int>(),
65  vector<double> massIn = vector<double>(),
66  vector<int> allowedIn = vector<int>() ) :
67  iRadiator(iRadiatorIn), iRecoiler(iRecoilerIn), pTmax(pTmaxIn),
68  colType(colIn), chgType(chgIn), gamType(gamIn), weakType(weakTypeIn),
69  isrType(isrIn), system(systemIn), systemRec(systemIn), MEtype(MEtypeIn),
70  iMEpartner(iMEpartnerIn), weakPol(weakPolIn), isOctetOnium(isOctetOniumIn),
71  isHiddenValley(isHiddenValleyIn), colvType(colvTypeIn), MEmix(MEmixIn),
72  MEorder (MEorderIn), MEsplit(MEsplitIn), MEgluinoRec(MEgluinoRecIn),
73  isFlexible(isFlexibleIn), flavour(0), iAunt(0), mass(massIn),
74  idRadAft(idRadAftIn), idEmtAft(idEmtAftIn), iSpectator(iSpectatorIn),
75  allowedEmissions(allowedIn), iSiblings(iSiblingsIn) {
76  mRad = m2Rad = mRec = m2Rec = mDip = m2Dip = m2DipCorr = pT2 = m2 = z
77  = mFlavour = asymPol = flexFactor = phi = 0.;
78  sa1 = xa = phia1 = 0.;
79  pT2start = pT2stop = pT2Old = 0.;
80  }
81 
82  DireTimesEnd( const DireTimesEnd& dip )
83  : iRadiator(dip.iRadiator), iRecoiler(dip.iRecoiler), pTmax(dip.pTmax),
84  colType(dip.colType), chgType(dip.chgType), gamType(dip.gamType),
85  weakType(dip.weakType), isrType(dip.isrType), system(dip.system),
86  systemRec(dip.systemRec), MEtype(dip.MEtype), iMEpartner(dip.iMEpartner),
87  weakPol(dip.weakPol), isOctetOnium(dip.isOctetOnium),
88  isHiddenValley(dip.isHiddenValley), colvType(dip.colvType),
89  MEmix(dip.MEmix), MEorder(dip.MEorder), MEsplit(dip.MEsplit),
90  MEgluinoRec(dip.MEgluinoRec), isFlexible(dip.isFlexible),
91  flavour(dip.flavour), iAunt(dip.iAunt),
92  mRad(dip.mRad), m2Rad(dip.m2Rad), mRec(dip.mRec), m2Rec(dip.m2Rec),
93  mDip(dip.mDip), m2Dip(dip.m2Dip), m2DipCorr(dip.m2DipCorr), pT2(dip.pT2),
94  m2(dip.m2), z(dip.z), mFlavour(dip.mFlavour), asymPol(dip.asymPol),
95  flexFactor(dip.flexFactor), phi(dip.phi), pT2start(dip.pT2start),
96  pT2stop(dip.pT2stop), pT2Old(dip.pT2Old), sa1(dip.sa1), xa(dip.xa),
97  phia1(dip.phia1), mass(dip.mass), idRadAft(dip.idRadAft),
98  idEmtAft(dip.idEmtAft), iSpectator(dip.iSpectator),
99  allowedEmissions(dip.allowedEmissions), iSiblings(dip.iSiblings) {}
100 
101  DireTimesEnd & operator=(const DireTimesEnd& t) { if (this != &t)
102  { iRadiator = t.iRadiator; iRecoiler = t.iRecoiler; pTmax = t.pTmax;
103  colType = t.colType; chgType = t.chgType; gamType = t.gamType;
104  weakType = t.weakType; isrType = t.isrType; system = t.system;
105  systemRec = t.systemRec; MEtype = t.MEtype; iMEpartner = t.iMEpartner;
106  weakPol = t.weakPol; isOctetOnium = t.isOctetOnium;
107  isHiddenValley = t.isHiddenValley; colvType = t.colvType;
108  MEmix = t.MEmix; MEorder = t.MEorder; MEsplit = t.MEsplit;
109  MEgluinoRec = t.MEgluinoRec; isFlexible = t.isFlexible;
110  flavour = t.flavour; iAunt = t.iAunt;
111  mRad = t.mRad; m2Rad = t.m2Rad; mRec = t.mRec; m2Rec = t.m2Rec;
112  mDip = t.mDip; m2Dip = t.m2Dip; m2DipCorr = t.m2DipCorr; pT2 = t.pT2;
113  m2 = t.m2; z = t.z; mFlavour = t.mFlavour; asymPol = t.asymPol;
114  flexFactor = t.flexFactor; phi = t.phi; pT2start = t.pT2start;
115  pT2stop = t.pT2stop; pT2Old = t.pT2Old; sa1 = t.sa1; xa = t.xa;
116  phia1 = t.phia1; mass = t.mass; idRadAft = t.idRadAft;
117  idEmtAft = t.idEmtAft; iSpectator = t.iSpectator;
118  allowedEmissions = t.allowedEmissions; iSiblings = t.iSiblings;}
119  return *this; }
120 
121  // Basic properties related to dipole and matrix element corrections.
122  int iRadiator, iRecoiler;
123  double pTmax;
124  int colType, chgType, gamType, weakType, isrType, system, systemRec,
125  MEtype, iMEpartner, weakPol;
126  bool isOctetOnium, isHiddenValley;
127  int colvType;
128  double MEmix;
129  bool MEorder, MEsplit, MEgluinoRec, isFlexible;
130 
131  // Properties specific to current trial emission.
132  int flavour, iAunt;
133  double mRad, m2Rad, mRec, m2Rec, mDip, m2Dip, m2DipCorr,
134  pT2, m2, z, mFlavour, asymPol, flexFactor, phi, pT2start, pT2stop,
135  pT2Old;
136 
137  // Properties of 1->3 splitting.
138  double sa1, xa, phia1;
139 
140  // Stored masses.
141  vector<double> mass;
142 
143  int idRadAft, idEmtAft;
144 
145  // Extended list of recoilers.
146  vector<int> iSpectator;
147  // List of allowed emissions (to avoid double-counting, since one
148  // particle can be part of many different dipoles.
149  void appendAllowedEmt( int id) {
150  if ( find(allowedEmissions.begin(), allowedEmissions.end(), id)
151  == allowedEmissions.end() ) allowedEmissions.push_back(id);
152  }
153  void removeAllowedEmt( int id) {
154  if ( find(allowedEmissions.begin(), allowedEmissions.end(), id)
155  != allowedEmissions.end() ) allowedEmissions.erase (
156  remove(allowedEmissions.begin(), allowedEmissions.end(), id),
157  allowedEmissions.end());
158  }
159  void clearAllowedEmt() { allowedEmissions.resize(0); }
160  vector<int> allowedEmissions;
161  bool canEmit() { return int(allowedEmissions.size() > 0); }
162 
163  void init(const Event& state) {
164  mRad = state[iRadiator].m();
165  mRec = state[iRecoiler].m();
166  mDip = sqrt( abs(2. * state[iRadiator].p() * state[iRecoiler].p()));
167  m2Rad = pow2(mRad);
168  m2Rec = pow2(mRec);
169  m2Dip = pow2(mDip);
170  }
171 
172  void list() const {
173  // Header.
174  cout << "\n -------- Begin DireTimesEnd Listing ----------------"
175  << "------------------------------------------------------- \n \n "
176  << " rad rec pTmax col parent1 parent2 isr"
177  << " sys sysR type MErec pol m2 allowedIds\n"
178  << fixed << setprecision(3);
179  cout << scientific << setprecision(4)
180  << setw(7) << iRadiator
181  << setw(7) << iRecoiler
182  << setw(12)<< pTmax
183  << setw(5) << colType
184  << setw(5) << isrType
185  << setw(5) << system << setw(5) << systemRec
186  << setw(5) << MEtype << setw(7) << iMEpartner
187  << setw(5) << weakPol
188  << setw(12) << m2Dip;
189  for (int j = 0; j < int(allowedEmissions.size()); ++j)
190  cout << setw(5) << allowedEmissions[j] << " ";
191  cout << endl;
192  // Done.
193  cout << "\n -------- End DireTimesEnd Listing ------------"
194  << "-------------------------------------------------------" << endl;
195  }
196 
197  DireSingleColChain iSiblings;
198  void setSiblings(DireSingleColChain s) { clearSiblings(); iSiblings = s; }
199  void clearSiblings() { iSiblings.clear(); }
200 
201  friend bool operator==(const DireTimesEnd& dip1, const DireTimesEnd& dip2);
202 
203 };
204 
205 //==========================================================================
206 
207 // The DireTimes class does timelike showers.
208 
209 class DireTimes : public TimeShower {
210 
211 public:
212 
213  // Constructor.
214  DireTimes() {}
215 
216  DireTimes( MergingHooksPtr mergingHooksPtrIn, PartonVertexPtr ) {
217  mergingHooksPtr = mergingHooksPtrIn;
218  beamOffset = 0;
219  userHooksPtr = 0;
220  splittingsPtr = 0;
221  weights = 0;
222  direInfoPtr = 0;
223  printBanner = true;
224  isInitSave = false;
225  usePDFalphas = false;
226  usePDF = true;
227  useSystems = true;
228  suppressLargeMECs = false;
229  }
230 
231  // Destructor.
232  virtual ~DireTimes() {}
233 
234  // Initialize alphaStrong and related pTmin parameters.
235  virtual void init( BeamParticle* beamAPtrIn = 0,
236  BeamParticle* beamBPtrIn = 0);
237 
238  bool initSplits() {
239  if (splittingsPtr) splits = splittingsPtr->getSplittings();
240  return (splits.size() > 0);
241  }
242 
243  // Initialize various pointers.
244  // (Separated from rest of init since not virtual.)
245  void reinitPtr(Info* infoPtrIn, MergingHooksPtr mergingHooksPtrIn,
246  DireSplittingLibrary* splittingsPtrIn, DireInfo* direInfoPtrIn) {
247  infoPtr = infoPtrIn;
248  settingsPtr = infoPtr->settingsPtr;
249  particleDataPtr = infoPtr->particleDataPtr;
250  rndmPtr = infoPtr->rndmPtr;
251  partonSystemsPtr = infoPtr->partonSystemsPtr;
252  userHooksPtr = infoPtr->userHooksPtr;
253  mergingHooksPtr = mergingHooksPtrIn;
254  splittingsPtr = splittingsPtrIn;
255  direInfoPtr = direInfoPtrIn;
256  }
257 
258  void initVariations();
259 
260  // Reset parton shower.
261  void clear();
262 
263  void setWeightContainerPtr(DireWeightContainer* weightsIn) {
264  weights = weightsIn;}
265 
266  // Find whether to limit maximum scale of emissions, and whether to dampen.
267  virtual bool limitPTmax( Event& event, double Q2Fac = 0.,
268  double Q2Ren = 0.);
269 
270  // Potential enhancement factor of pTmax scale for hardest emission.
271  virtual double enhancePTmax() { return pTmaxFudge;}
272 
273  // Top-level routine to do a full time-like shower in resonance decay.
274  virtual int shower( int iBeg, int iEnd, Event& event, double pTmax,
275  int nBranchMax = 0);
276 
277  // Top-level routine for QED radiation in hadronic decay to two leptons.
278  virtual int showerQED( int i1, int i2, Event& event, double pTmax);
279 
280  // Global recoil: reset counters and store locations of outgoing partons.
281  virtual void prepareGlobal( Event&);
282 
283  // Prepare system for evolution after each new interaction; identify ME.
284  virtual void prepare( int iSys, Event& event, bool limitPTmaxIn = true);
285 
286  // Finalize event after evolution.
287  void finalize( Event& event);
288 
289  // Update dipole list after a multiparton interactions rescattering.
290  virtual void rescatterUpdate( int iSys, Event& event);
291 
292  // Update dipole list after each ISR emission.
293  virtual void update( int iSys, Event& event, bool = false);
294 
295  // Update dipole list after final-final splitting.
296  void updateAfterFF( int iSysSelNow, int iSysSelRec,
297  Event& event, int iRadBef, int iRecBef, int iRad, int iEmt, int iRec,
298  int flavour, int colType, double pTsel);
299 
300  // Update dipole list after final-final splitting.
301  void updateAfterFI( int iSysSelNow, int iSysSelRec,
302  Event& event, int iRadBef, int iRecBef, int iRad, int iEmt, int iRec,
303  int flavour, int colType, double pTsel, double xNew);
304 
305  // Select next pT in downwards evolution. Wrapper function inherited from
306  // Pythia.
307  virtual double pTnext( Event& event, double pTbegAll, double pTendAll,
308  bool = false, bool = false);
309  double newPoint( const Event& event);
310  // Setup branching kinematics.
311  virtual bool branch( Event& event, bool = false);
312  bool branch_FF( Event& event, bool = false,
313  DireSplitInfo* split = NULL);
314  bool branch_FI( Event& event, bool = false,
315  DireSplitInfo* split = NULL);
316 
317  pair < Vec4, Vec4 > decayWithOnshellRec( double zCS, double yCS, double phi,
318  double m2Rec, double m2RadAft, double m2EmtAft,
319  Vec4 pRadBef, Vec4 pRecBef);
320  pair < Vec4, Vec4 > decayWithOffshellRec( double zCS, double yCS, double phi,
321  double m2RadBef, double m2RadAft, double m2EmtAft,
322  Vec4 pRadBef, Vec4 pRecBef);
323 
324  bool getHasWeaklyRadiated() {return false;}
325  // Tell which system was the last processed one.
326  int system() const {return iSysSel;};
327 
328  // Setup clustering kinematics.
329  virtual Event clustered( const Event& state, int iRad, int iEmt, int iRec,
330  string name) {
331  pair <Event, pair<int,int> > reclus
332  = clustered_internal(state, iRad, iEmt, iRec, name);
333  if (reclus.first.size() > 0)
334  reclus.first[0].mothers(reclus.second.first,reclus.second.second);
335  return reclus.first;
336  }
337  pair <Event, pair<int,int> > clustered_internal( const Event& state,
338  int iRad, int iEmt, int iRec, string name);
339  bool cluster_FF( const Event& state, int iRad,
340  int iEmt, int iRec, int idRadBef, Particle& radBef, Particle& recBef);
341  bool cluster_FI( const Event& state, int iRad,
342  int iEmt, int iRec, int idRadBef, Particle& radBef, Particle& recBef);
343 
344  // From Pythia version 8.215 onwards no longer virtual.
345  double pT2Times ( const Particle& rad, const Particle& emt,
346  const Particle& rec) {
347  if (rec.isFinal()) return pT2_FF(rad,emt,rec);
348  return pT2_FI(rad,emt,rec);
349  }
350 
351  double pT2_FF ( const Particle& rad, const Particle& emt,
352  const Particle& rec);
353  double pT2_FI ( const Particle& rad, const Particle& emt,
354  const Particle& rec);
355 
356  // From Pythia version 8.215 onwards no longer virtual.
357  double zTimes ( const Particle& rad, const Particle& emt,
358  const Particle& rec) {
359  if (rec.isFinal()) return z_FF(rad,emt,rec);
360  return z_FI(rad,emt,rec);
361  }
362 
363  double z_FF ( const Particle& rad, const Particle& emt,
364  const Particle& rec);
365  double z_FI ( const Particle& rad, const Particle& emt,
366  const Particle& rec);
367  double z_FF_fromVec ( const Vec4& rad, const Vec4& emt, const Vec4& rec);
368 
369  double m2dipTimes ( const Particle& rad, const Particle& emt,
370  const Particle& rec) {
371  if (rec.isFinal()) return m2dip_FF(rad,emt,rec);
372  return m2dip_FI(rad,emt,rec);
373  }
374 
375  double m2dip_FF ( const Particle& rad, const Particle& emt,
376  const Particle& rec);
377  double m2dip_FI ( const Particle& rad, const Particle& emt,
378  const Particle& rec);
379 
380  // From Pythia version 8.218 onwards.
381  // Return the evolution variable.
382  // Usage: getStateVariables( const Event& event, int iRad, int iEmt,
383  // int iRec, string name)
384  // Important note:
385  // - This map must contain an entry for the shower evolution variable,
386  // specified with key "t".
387  // - This map must contain an entry for the shower evolution variable from
388  // which the shower would be restarted after a branching. This entry
389  // must have key "tRS",
390  // - This map must contain an entry for the argument of \alpha_s used
391  // for the branching. This entry must have key "scaleAS".
392  // - This map must contain an entry for the argument of the PDFs used
393  // for the branching. This entry must have key "scalePDF".
394  virtual map<string, double> getStateVariables (const Event& state,
395  int rad, int emt, int rec, string name);
396 
397  // From Pythia version 8.215 onwards.
398  // Check if attempted clustering is handled by timelike shower
399  // Usage: isTimelike( const Event& event, int iRad, int iEmt,
400  // int iRec, string name)
401  virtual bool isTimelike(const Event& state, int iRad, int, int, string)
402  { return state[iRad].isFinal(); }
403 
404  // From Pythia version 8.215 onwards.
405  // Return a string identifier of a splitting.
406  // Usage: getSplittingName( const Event& event, int iRad, int iEmt, int iRec)
407  virtual vector<string> getSplittingName( const Event& state, int iRad,
408  int iEmt, int)
409  { return splittingsPtr->getSplittingName(state,iRad,iEmt); }
410 
411  // From Pythia version 8.215 onwards.
412  // Return the splitting probability.
413  // Usage: getSplittingProb( const Event& event, int iRad, int iEmt, int iRec)
414  virtual double getSplittingProb( const Event& state, int iRad,
415  int iEmt, int iRec, string name);
416 
417  virtual bool allowedSplitting( const Event& state, int iRad, int iEmt);
418 
419  virtual vector<int> getRecoilers( const Event& state, int iRad, int iEmt,
420  string name);
421 
422  virtual double getCoupling( double mu2Ren, string name) {
423  if (splits.find(name) != splits.end())
424  return splits[name]->coupling(-1.,mu2Ren, 0., 1.);
425  return 1.;
426  }
427 
428  bool isSymmetric( string name, const Particle* rad, const Particle* emt) {
429  if (splits.find(name) != splits.end())
430  return splits[name]->isSymmetric(rad,emt);
431  return false;
432  }
433 
434  // Auxiliary function to return the position of a particle.
435  // Should go int Event class eventually!
436  int FindParticle( const Particle& particle, const Event& event,
437  bool checkStatus = true );
438 
439  // Print dipole list; for debug mainly.
440  virtual void list() const;
441 
442  Event makeHardEvent( int iSys, const Event& state, bool isProcess = false );
443 
444  // Check that particle has sensible momentum.
445  bool validMomentum( const Vec4& p, int id, int status);
446 
447  // Check colour/flavour correctness of state.
448  bool validEvent( const Event& state, bool isProcess = false,
449  int iSysCheck = -1 );
450 
451  // Check that mother-daughter-relations are correctly set.
452  bool validMotherDaughter( const Event& state );
453 
454  // Find index of colour partner for input colour.
455  int FindCol(int col, vector<int> iExclude, const Event& event, int type,
456  int iSys = -1);
457 
458  // Pointers to the two incoming beams.
459  BeamParticle* getBeamA () { return beamAPtr; }
460  BeamParticle* getBeamB () { return beamBPtr; }
461 
462  // Function to calculate the correct alphaS/2*Pi value, including
463  // renormalisation scale variations + threshold matching.
464  double alphasNow( double pT2, double renormMultFacNow = 1., int iSys = 0 );
465 
466  // Function to calculate the correct alphaEM/2*Pi value.
467  double alphaemNow( double pT2, double renormMultFacNow = 1., int iSys = 0 );
468 
469  bool isInit() { return isInitSave; }
470 
471  // Function to calculate the absolute phase-sace boundary for emissions.
472  double m2Max (int iDip, const Event& state) {
473  if ( state[dipEnd[iDip].iRecoiler].isFinal()
474  && state[dipEnd[iDip].iRadiator].isFinal())
475  return dipEnd[iDip].m2Dip;
476  int iSys = dipEnd[iDip].system;
477  int inA = partonSystemsPtr->getInA(iSys);
478  int inB = partonSystemsPtr->getInB(iSys);
479  double x = 1.;
480  int iRad(dipEnd[iDip].iRadiator), iRec(dipEnd[iDip].iRecoiler);
481  if (hasPDF(state[iRad].id()) && iRad == inA)
482  x *= state[inA].pPos() / state[0].m();
483  if (hasPDF(state[iRad].id()) && iRad == inB)
484  x *= state[inB].pNeg() / state[0].m();
485  if (hasPDF(state[iRec].id()) && iRec == inA)
486  x *= state[inA].pPos() / state[0].m();
487  if (hasPDF(state[iRec].id()) && iRec == inB)
488  x *= state[inB].pNeg() / state[0].m();
489  return dipEnd[iDip].m2Dip/x;
490  }
491 
492  bool dryrun;
493 
494 private:
495 
496  friend class DireSplitting;
497  friend class DireSpace;
498 
499  // Number of times the same error message is repeated, unless overridden.
500  static const int TIMESTOPRINT;
501 
502  // Allow conversion from mb to pb.
503  static const double CONVERTMB2PB;
504 
505  // Colour factors.
506  //static const double CA, CF, TR, NC, LEPTONZMAX;
507  static const double LEPTONZMAX;
508  double CA, CF, TR, NC;
509 
510  // Store common beam quantities.
511  int idASave, idBSave;
512 
513 protected:
514 
515  // Store properties to be returned by methods.
516  int iSysSel;
517  double pTmaxFudge, pTLastBranch;
518 
519 private:
520 
521  // Constants: could only be changed in the code itself.
522  static const int MAXLOOPTINYPDF;
523  static const double MCMIN, MBMIN, SIMPLIFYROOT, XMARGIN, XMARGINCOMB,
524  TINYPDF, LARGEM2, THRESHM2, LAMBDA3MARGIN, TINYMASS, TINYOVERESTIMATE,
525  PT2MIN_PDF_IN_OVERESTIMATE, PT2_INCREASE_OVERESTIMATE,
526  KERNEL_HEADROOM;
527 
528  // Initialization data, normally only set once.
529  bool isInitSave, doQCDshower, doQEDshowerByQ, doQEDshowerByL,
530  doMEcorrections, doMEafterFirst, doPhiPolAsym,
531  doInterleave, allowBeamRecoil, dampenBeamRecoil, recoilToColoured,
532  useFixedFacScale, allowRescatter, canVetoEmission, hasUserHooks,
533  doSecondHard, alphaSuseCMW, printBanner, doTrialNow;
534  int pTmaxMatch, pTdampMatch, alphaSorder, alphaSnfmax, alphaEMorder,
535  nGluonToQuark, nGammaToQuark, nGammaToLepton, nFinalMax,
536  nFinalMaxMECs,kernelOrder, kernelOrderMPI, nMPI, asScheme;
537  double pTdampFudge, mc, mb, m2c, m2b, renormMultFac, factorMultFac,
538  fixedFacScale2, alphaSvalue, alphaS2pi, Lambda3flav, Lambda4flav,
539  Lambda5flav, Lambda3flav2, Lambda4flav2, Lambda5flav2,
540  pTcolCutMin, pTcolCut,
541  pT2colCut, m2colCut, mTolErr, mZ, gammaZ, thetaW, mW, gammaW,
542  pTmaxFudgeMPI, sumCharge2L, sumCharge2Q, sumCharge2Tot,
543  pT2minVariations, pT2minEnhance, pT2minMECs, Q2minMECs, pT2recombine,
544  m2cPhys, m2bPhys;
545  double alphaS2piOverestimate;
546  bool usePDFalphas, usePDFmasses, useSummedPDF, usePDF, useSystems,
547  useMassiveBeams, suppressLargeMECs;
548 
549  double pTchgQCut, pT2chgQCut, pTchgLCut, pT2chgLCut;
550 
551  unordered_map<int,double> pT2cutSave;
552  double pT2cut(int id) {
553  if (pT2cutSave.find(id) != pT2cutSave.end()) return pT2cutSave[id];
554  // Else return maximal value.
555  double ret = 0.;
556  for ( unordered_map<int,double>::iterator it = pT2cutSave.begin();
557  it != pT2cutSave.end(); ++it ) ret = max(ret, it->second);
558  return ret;
559  }
560  double pT2cutMax(DireTimesEnd* dip) {
561  double ret = 0.;
562  for (int i=0; i < int(dip->allowedEmissions.size()); ++i)
563  ret = max( ret, pT2cut(dip->allowedEmissions[i]));
564  return ret;
565  }
566  double pT2cutMin(DireTimesEnd* dip) {
567  double ret = 1e15;
568  for (int i=0; i < int(dip->allowedEmissions.size()); ++i)
569  ret = min( ret, pT2cut(dip->allowedEmissions[i]));
570  return ret;
571  }
572 
573  bool doDecaysAsShower;
574 
575  // alphaStrong and alphaEM calculations.
576  AlphaStrong alphaS;
577  AlphaEM alphaEM;
578 
579  // Some current values.
580  bool dopTlimit1, dopTlimit2, dopTdamp;
581  double pT2damp, kRad, kEmt, pdfScale2;
582 
583  // All dipole ends and a pointer to the selected hardest dipole end.
584  vector<DireTimesEnd> dipEnd;
585  DireTimesEnd* dipSel;
586  DireSplitInfo splitInfoSel;
587  DireSplitting* splittingSel;
588  int iDipSel;
589  unordered_map<string,double> kernelSel, kernelNow;
590  double auxSel, overSel, boostSel, auxNow, overNow, boostNow;
591 
592  double tinypdf( double x) {
593  double xref = 0.01;
594  return TINYPDF*log(1-x)/log(1-xref);
595  }
596 
597  // Function to check if id is part of the incoming hadron state.
598  bool hasPDF (int id) {
599  if ( !usePDF ) return false;
600  if ( particleDataPtr->colType(id) != 0) return true;
601  if ( particleDataPtr->isLepton(id)
602  && settingsPtr->flag("PDF:lepton")) return true;
603  return false;
604  }
605 
606  // Wrapper around PDF calls.
607  double getXPDF( int id, double x, double t, int iSys = 0,
608  BeamParticle* beam = NULL, bool finalRec = true, double z = 0.,
609  double m2dip = 0.) {
610  // Return one if no PDF should be used.
611  if (!hasPDF(id)) return 1.0;
612  // Else get PDF from beam particle.
613  BeamParticle* b = beam;
614  if (b == NULL) {
615  if (beamAPtr != NULL || beamBPtr != NULL) {
616  b = (beamAPtr != NULL && particleDataPtr->isHadron(beamAPtr->id()))
617  ? beamAPtr
618  : (beamBPtr != NULL && particleDataPtr->isHadron(beamBPtr->id()))
619  ? beamBPtr : NULL;
620  }
621  if (b == NULL && beamAPtr != 0) beam = beamAPtr;
622  if (b == NULL && beamBPtr != 0) beam = beamBPtr;
623  }
624 
625  double scale2 = t;
626  if (asScheme == 2 && z != 0. && finalRec) {
627  double zcs = z;
628  double xcs = m2dip * zcs * (1.-zcs) / (t + m2dip * zcs * (1.-zcs));
629  scale2 = (1-zcs)*(1-xcs)/xcs/zcs*m2dip;
630  }
631 
632  double ret = (useSummedPDF) ? b->xf(id, x, scale2)
633  : b->xfISR(iSys,id, x, scale2);
634  // Done.
635  return ret;
636  }
637 
638  // Evolve a QCD dipole end near heavy quark threshold region.
639  // Setup a dipole end, either QCD, QED/photon, weak or Hidden Valley one.
640  void setupQCDdip( int iSys, int i, int colTag, int colSign, Event& event,
641  bool isOctetOnium = false, bool limitPTmaxIn = true);
642  void getGenDip( int iSys, int i, int iRad, const Event& event,
643  bool limitPTmaxIn, vector<DireTimesEnd>& dipEnds );
644  void getQCDdip( int iRad, int colTag, int colSign, const Event& event,
645  vector<DireTimesEnd>& dipEnds );
646  void setupDecayDip( int iSys, int iRad, const Event& event,
647  vector<DireTimesEnd>& dipEnds);
648 
649  // Function to set up and append a new dipole.
650  bool appendDipole( const Event& state, int iRad, int iRec, double pTmax,
651  int colType, int chgType, int gamType, int weakType, int isrType, int iSys,
652  int MEtype, int iMEpartner, int weakPol, bool isOctetOnium,
653  vector<DireTimesEnd>& dipEnds);
654 
655  vector<int> sharedColor(const Particle& rad, const Particle& rec);
656 
657  // Function to set up and append a new dipole.
658  void updateDipoles(const Event& state, int iSys = -1);
659  void checkDipoles(const Event& state);
660  void saveSiblings(const Event& state, int iSys = -1);
661  bool updateAllowedEmissions( const Event& state, DireTimesEnd* dip);
662  bool appendAllowedEmissions( const Event& state, DireTimesEnd* dip);
663 
664  // Evolve a QCD dipole end.
665  void pT2nextQCD( double pT2begDip, double pT2sel, DireTimesEnd& dip,
666  Event& event, double pT2endForce = -1., double pT2freeze = 0.,
667  bool forceBranching = false);
668  bool pT2nextQCD_FF( double pT2begDip, double pT2sel, DireTimesEnd& dip,
669  const Event& event, double pT2endForce = -1., double pT2freeze = 0.,
670  bool forceBranching = false);
671  bool pT2nextQCD_FI( double pT2begDip, double pT2sel, DireTimesEnd& dip,
672  const Event& event, double pT2endForce = -1., double pT2freeze = 0.,
673  bool forceBranching = false);
674 
675  double tNextQCD( DireTimesEnd*, double overestimateInt,
676  double tOld, double tMin, double tFreeze=0., int algoType = 0);
677  bool zCollNextQCD( DireTimesEnd* dip, double zMin, double zMax,
678  double tMin = 0., double tMax = 0.);
679  bool virtNextQCD( DireTimesEnd* dip, double tMin, double tMax,
680  double zMin =-1., double zMax =-1.);
681 
682  // Function to determine how often the integrated overestimate should be
683  // recalculated.
684  double evalpdfstep(int idRad, double pT2, double m2cp = -1.,
685  double m2bp = -1.) {
686  double ret = 0.2;
687  if (m2cp < 0.) m2cp = particleDataPtr->m0(4);
688  if (m2bp < 0.) m2bp = particleDataPtr->m0(5);
689  // More steps close to the thresholds.
690  if ( abs(idRad) == 4 && pT2 < 1.2*m2cp && pT2 > m2cp) ret = 1.0;
691  if ( abs(idRad) == 5 && pT2 < 1.2*m2bp && pT2 > m2bp) ret = 1.0;
692  return ret;
693  }
694 
695  DireSplittingLibrary* splittingsPtr;
696 
697  // Number of proposed splittings in hard scattering systems.
698  unordered_map<int,int> nProposedPT;
699 
700  // Return headroom factors for integrated/differential overestimates.
701  double overheadFactors(DireTimesEnd*, const Event&, string, double,
702  double, double);
703  double enhanceOverestimateFurther( string, int, double );
704  double overheadFactorsMEC(const Event&, DireSplitInfo*, string);
705 
706  // Function to fill map of integrated overestimates.
707  void getNewOverestimates( DireTimesEnd*, const Event&, double, double,
708  double, double, multimap<double,string>&);
709 
710  // Function to sum all integrated overestimates.
711  void addNewOverestimates( multimap<double,string>, double&);
712 
713  // Function to attach the correct alphaS weights to the kernels.
714  void alphasReweight(double t, double talpha, int iSys, bool forceFixedAs,
715  double& weight, double& fullWeight, double& overWeight,
716  double renormMultFacNow);
717 
718  // Function to evaluate the accept-probability, including picking of z.
719  void getNewSplitting( const Event&, DireTimesEnd*, double, double, double,
720  double, double, int, string, bool, int&, int&, double&, double&,
721  unordered_map<string,double>&, double&);
722 
723  pair<bool, pair<double,double> > getMEC ( const Event& state,
724  DireSplitInfo* splitInfo);
725  bool applyMEC ( const Event& state, DireSplitInfo* splitInfo,
726  vector<Event> auxEvent = vector<Event>() );
727 
728  // Get particle masses.
729  double getMass(int id, int strategy, double mass = 0.) {
730  BeamParticle* beam = NULL;
731  if (beamAPtr != NULL || beamBPtr != NULL) {
732  beam = (beamAPtr != NULL && particleDataPtr->isHadron(beamAPtr->id()))
733  ? beamAPtr
734  : (beamBPtr != NULL && particleDataPtr->isHadron(beamBPtr->id()))
735  ? beamBPtr : NULL;
736  }
737  bool usePDFmass = usePDFmasses
738  && (toLower(settingsPtr->word("PDF:pSet")).find("lhapdf")
739  != string::npos);
740  double mRet = 0.;
741  // Parton masses.
742  if ( particleDataPtr->colType(id) != 0) {
743  if (strategy == 1) mRet = particleDataPtr->m0(id);
744  if (strategy == 2 && usePDFmass && beam != NULL)
745  mRet = beam->mQuarkPDF(id);
746  if (strategy == 2 && (!usePDFmass || beam == NULL))
747  mRet = particleDataPtr->m0(id);
748  if (strategy == 3) mRet = mass;
749  if (mRet < TINYMASS) mRet = 0.;
750  // Masses of other particles.
751  } else {
752  mRet = particleDataPtr->m0(id);
753  if (strategy == 3) mRet = mass;
754  if (mRet < TINYMASS) mRet = 0.;
755  }
756  return pow2(max(0.,mRet));
757  }
758 
759  // Check if variables are in allowed phase space.
760  bool inAllowedPhasespace(int kinType, double z, double pT2, double m2dip,
761  double q2, double xOld, int splitType = 0, double m2RadBef = 0.,
762  double m2r = 0., double m2s = 0., double m2e = 0.,
763  vector<double> aux = vector<double>());
764 
765  // Auxiliary function to get number of flavours.
766  double getNF(double pT2);
767 
768  // Auxiliary functions to get beta function coefficients.
769  double beta0 (double NF)
770  { return 11./6.*CA - 2./3.*NF*TR; }
771  double beta1 (double NF)
772  { return 17./6.*pow2(CA) - (5./3.*CA+CF)*NF*TR; }
773  double beta2 (double NF)
774  { return 2857./432.*pow(CA,3)
775  + (-1415./216.*pow2(CA) - 205./72.*CA*CF + pow2(CF)/4.) *TR*NF
776  + ( 79.*CA + 66.*CF)/108.*pow2(TR*NF); }
777 
778  // Identifier of the splitting
779  string splittingNowName, splittingSelName;
780 
781  // Weighted shower book-keeping.
782  unordered_map<string, map<double,double> > acceptProbability;
783  unordered_map<string, multimap<double,double> > rejectProbability;
784 
785 public:
786 
787  DireWeightContainer* weights;
788  DireInfo* direInfoPtr;
789  ProcessLevel processLevel;
790  unordered_map<string, DireSplitting* > splits;
791  vector<int> bornColors;
792 
793 private:
794 
795  bool doVariations;
796 
797  // List of splitting kernels.
798  //map<string, DireSplitting* > splits;
799  unordered_map<string, double > overhead;
800  void scaleOverheadFactor(string name, double scale) {
801  overhead[name] *= scale;
802  return;
803  }
804  void resetOverheadFactors() {
805  for ( unordered_map<string,double>::iterator it = overhead.begin();
806  it != overhead.end(); ++it )
807  it->second = 1.0;
808  return;
809  }
810 
811  double octetOniumColFac;
812  bool useLocalRecoilNow;
813 
814  // Map to store some settings, to be passes to splitting kernels.
815  unordered_map<string,bool> bool_settings;
816 
817 };
818 
819 //==========================================================================
820 
821 } // end namespace Pythia8
822 
823 #endif // end Pythia8_DireTimes_H
Definition: beam.h:43