StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
DireSplittings.h
1 // DireSplittings.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 OverheadInfo and DireSplitting classes.
7 
8 #ifndef Pythia8_DireSplittings_H
9 #define Pythia8_DireSplittings_H
10 
11 #define DIRE_SPLITTINGS_VERSION "2.002"
12 
13 #include "Pythia8/Basics.h"
14 #include "Pythia8/BeamParticle.h"
15 #include "Pythia8/ParticleData.h"
16 #include "Pythia8/PythiaStdlib.h"
17 #include "Pythia8/Settings.h"
18 #include "Pythia8/StandardModel.h"
19 #include "Pythia8/Info.h"
20 #include "Pythia8/DireSplitInfo.h"
21 #include "Pythia8/DireBasics.h"
22 
23 namespace Pythia8 {
24 
25 class DireSpace;
26 class DireTimes;
27 
28 //==========================================================================
29 
30 class OverheadInfo {
31 
32 public:
33 
34  OverheadInfo(int nFinalIn, int idIn, double valIn, double xIn, double pT2In)
35  : nFinal(nFinalIn), id(idIn), val(valIn), x(xIn), pT2(pT2In) {}
36 
37  int nFinal, id;
38  double val, x, pT2;
39 
40  bool match(int idIn, int nfIn) { return (idIn==id && nfIn==nFinal); }
41 
42  string list () const {
43  ostringstream os;
44  os << scientific << setprecision(6)
45  << "pT2 " << setw(10) << pT2 << " x " << setw(10) << x
46  << " id " << setw(4) << id << " nf " << setw(4) << nFinal
47  << " val=" << val;
48  return os.str();
49  }
50 
51 };
52 
54 
55 public:
56 
57  // Constructor and destructor.
58  DireSplitting() : renormMultFac(0),
59  id("void"), correctionOrder(0), settingsPtr(0),
60  particleDataPtr(0), rndmPtr(0), beamAPtr(0),
61  beamBPtr(0), coupSMPtr(0), infoPtr(0), direInfoPtr(0),
62  is_qcd(false), is_qed(false), is_ewk(false), is_fsr(false),
63  is_isr(false), is_dire(false), nameHash(0) {}
64  DireSplitting(string idIn, int softRS, Settings* settings,
65  ParticleData* particleData, Rndm* rndm, BeamParticle* beamA,
66  BeamParticle* beamB, CoupSM* coupSMPtrIn, Info* infoPtrIn,
67  DireInfo* direInfo) :
68  renormMultFac(0), id(idIn), correctionOrder(softRS),
69  settingsPtr(settings), particleDataPtr(particleData), rndmPtr(rndm),
70  beamAPtr(beamA), beamBPtr(beamB), coupSMPtr(coupSMPtrIn),
71  infoPtr(infoPtrIn), direInfoPtr(direInfo), is_qcd(false), is_qed(false),
72  is_ewk(false), is_fsr(false), is_isr(false), is_dire(false),
73  nameHash(0) { init(); splitInfo.storeName(name()); }
74  virtual ~DireSplitting() {}
75 
76  void init();
77 
78 public:
79 
80  double renormMultFac;
81 
82  string id;
83  int correctionOrder;
84  Settings* settingsPtr;
85  ParticleData* particleDataPtr;
86  Rndm* rndmPtr;
87  BeamParticle* beamAPtr;
88  BeamParticle* beamBPtr;
89  CoupSM* coupSMPtr;
90  Info* infoPtr;
91  DireInfo* direInfoPtr;
92 
93  // Some short-cuts and string hashes to help avoid string comparisons.
94  bool is_qcd, is_qed, is_ewk, is_fsr, is_isr, is_dire;
95  ulong nameHash;
96  bool is( ulong pattern ) {
97  if (pattern == nameHash) return true;
98  return false;
99  }
100 
101  unordered_map<string,double> kernelVals;
102 
103  string name () {return id;}
104 
105  virtual bool canRadiate ( const Event&, pair<int,int>,
106  unordered_map<string,bool> = unordered_map<string,bool>(),
107  Settings* = NULL, PartonSystems* = NULL, BeamParticle* = NULL)
108  {return false;}
109 
110  // Discard below the cut-off for the splitting.
111  virtual bool aboveCutoff( double, const Particle&, const Particle&, int,
112  PartonSystems* = NULL) { return true; }
113 
114  virtual bool useFastFunctions() { return false; }
115  virtual bool canRadiate ( const Event&, int, int,
116  Settings* = NULL, PartonSystems* = NULL, BeamParticle* = NULL)
117  {return false;}
118 
119  // Function to return an identifier for the phase space mapping
120  // that is used for setting up this splitting.
121  // return values: 1 --> Default Dire mapping.
122  // 2 --> Dire 1->3 mapping.
123  virtual int kinMap () {return 1;}
124 
125  // Return id of mother after splitting.
126  virtual int motherID(int) {return 0;}
127 
128  // Return id of emission.
129  virtual int sisterID(int) {return 0;}
130 
131  // Return a pair of ids for the radiator and emission after
132  // the splitting.
133  virtual vector <int> radAndEmt(int, int) { return vector<int>(); }
134  virtual vector < pair<int,int> > radAndEmtCols( int, int, Event)
135  { return vector<pair<int,int> >(); }
136  virtual bool canUseForBranching() { return false; }
137  virtual bool isPartial() { return false; }
138  virtual int nEmissions() { return 0; }
139 
140  virtual bool swapRadEmt() { return false; }
141  virtual bool isSymmetric( const Particle* = NULL, const Particle* = NULL)
142  { return false; }
143 
144  // Return a vector of all possible recoiler positions, given the
145  // positions of the radiator and emission after the splitting.
146  virtual vector <int> recPositions( const Event&, int, int)
147  {return vector<int>();}
148 
149  // Return id of recombined radiator (before splitting!)
150  virtual int radBefID(int, int) {return 0;}
151 
152  // Return colours of recombined radiator (before splitting!)
153  virtual pair<int,int> radBefCols(int, int, int, int)
154  {return make_pair(0,0);}
155 
156  // Return color factor for splitting.
157  virtual double gaugeFactor (int, int) {return 1.;}
158 
159  // Return symmetry factor for splitting.
160  virtual double symmetryFactor (int, int) {return 1.;}
161 
162  // Return an identifier for the interaction that causes the
163  // branching.
164  // return values -1 --> Type not defined.
165  // 1 --> QCD splitting (i.e. proportional to alphaS)
166  // 2 --> QED splitting (i.e. proportional to alphaEM)
167  // 3 --> EW splitting (i.e. proportional to sinThetaW)
168  // 4 --> Yukawa splitting (i.e. proportional to y)
169  virtual int couplingType (int, int) {return -1;}
170 
171  // Return the value of the coupling that should be used for this branching.
172  // Note that the last input allows easy access to the PS evolution variable.
173  // return values -1 --> Coupling value not defined.
174  // double > 0 --> Value to be used for this branching.
175  virtual double coupling (double = 0., double = 0., double = 0.,
176  double = -1., pair<int,bool> = pair<int,bool>(),
177  pair<int,bool> = pair<int,bool>()) {
178  return -1.;
179  }
180  virtual double couplingScale2 (double = 0., double = 0., double = 0.,
181  pair<int,bool> = pair<int,bool>(), pair<int,bool> = pair<int,bool>()) {
182  return -1.;
183  }
184 
185  // Pick z for new splitting.
186  virtual double zSplit(double, double, double) {return 0.5;}
187 
188  // New overestimates, z-integrated versions.
189  virtual double overestimateInt(double, double, double, double, int = -1)
190  { return 0.;}
191 
192  // Return kernel for new splitting.
193  virtual double overestimateDiff(double, double, int = -1) {return 1.;}
194 
195  // Functions to store and retrieve all the variants of the kernel.
196  virtual double getKernel(string = "");
197  virtual unordered_map<string,double> getKernelVals() { return kernelVals; }
198  virtual void clearKernels() { kernelVals.clear(); }
199 
200  DireSplitInfo splitInfo;
201 
202  // Functions to calculate the kernel from SplitInfo information.
203  virtual bool calc(const Event& = Event(), int = -1) { return false; }
204 
205  shared_ptr<DireSpace> isr;
206  shared_ptr<DireTimes> fsr;
207  shared_ptr<DireTimes> fsrDec;
208  void setTimesPtr(shared_ptr<DireTimes> fsrIn) { fsr=fsrIn;}
209  void setTimesDecPtr(shared_ptr<DireTimes> fsrIn) { fsrDec=fsrIn;}
210  void setSpacePtr(shared_ptr<DireSpace> isrIn) { isr=isrIn;}
211 
212  // Functions that allow different ordering variables for emissions.
213  // Note: Only works after splitInfo has been properly filled.
214  virtual double getJacobian( const Event& = Event(),
215  PartonSystems* = 0) { return 0.;}
216  // Map filled identical to shower state variables map.
217  virtual unordered_map<string, double> getPhasespaceVars(
218  const Event& = Event(),
219  PartonSystems* = 0) { return unordered_map<string,double>(); }
220 
221  // Treatment of additional virtual corrections.
222  virtual bool allow_z_endpoint_for_kinematics() { return false; }
223  virtual bool allow_pT2_endpoint_for_kinematics() { return false; }
224  virtual bool allow_sai_endpoint_for_kinematics() { return false; }
225  virtual bool allow_xa_endpoint_for_kinematics() { return false; }
226  // Functions to set if kernel should contribute to a kinematical endpoint.
227  virtual void try_z_endpoint() { return; }
228  virtual void try_pT2_endpoint() { return; }
229  virtual void try_sai_endpoint() { return; }
230  virtual void try_xa_endpoint() { return; }
231  // Return endpoint information.
232  virtual bool is_z_endpoint() { return false; }
233  virtual bool is_pT2_endpoint() { return false; }
234  virtual bool is_sai_endpoint() { return false; }
235  virtual bool is_xa_endpoint() { return false; }
236 
237  // Functions to calculate Dire variables from the evolution variables.
238  virtual double tdire_ff(double, double t, double) { return t; }
239  virtual double zdire_ff(double z, double, double) { return z; }
240  virtual double tdire_fi(double, double t, double) { return t; }
241  virtual double zdire_fi(double z, double, double) { return z; }
242  virtual double tdire_if(double, double t, double) { return t; }
243  virtual double zdire_if(double z, double, double) { return z; }
244  virtual double tdire_ii(double, double t, double) { return t; }
245  virtual double zdire_ii(double z, double, double) { return z; }
246 
247  virtual bool hasMECBef(const Event&, double) { return false; }
248  virtual bool hasMECAft(const Event&, double) { return false; }
249 
250  virtual void storeOverhead(double pT2, double x, int radid, int nf,
251  double val) { overhead_map.insert(make_pair(pT2, OverheadInfo(nf, radid,
252  val, x, pT2))); }
253  multimap<double,OverheadInfo> overhead_map;
254 
255  virtual double overhead(double pT2, int idd, int nf) {
256 
257  if (overhead_map.empty()) return 1.;
258 
259  multimap<double,OverheadInfo>::iterator lo = overhead_map.lower_bound(pT2);
260  if (lo != overhead_map.begin()) lo--;
261  if (lo != overhead_map.begin()) lo--;
262  multimap<double,OverheadInfo>::iterator hi = overhead_map.upper_bound(pT2);
263  if (hi != overhead_map.end()) hi++;
264  if (hi == overhead_map.end()) hi--;
265 
266  int n(0);
267  double sum = 0.;
268  for ( multimap<double,OverheadInfo>::iterator it = lo;
269  it != hi; ++it ){
270  if (!it->second.match(idd,nf)) continue;
271  sum += it->second.val;
272  n++;
273  }
274 
275  if (hi->second.match(idd,nf)) {
276  sum += hi->second.val;
277  n++;
278  }
279 
280  return max(sum/max(1,n),1.);
281 
282  }
283 
284 };
285 
286 //==========================================================================
287 
288 } // end namespace Pythia8
289 
290 #endif // end Pythia8_DireSplittings_H