PYTHIA  8.312
aMCatNLOHooks.h
1 // aMCatNLOHooks.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 // This program is written by Stefan Prestel.
7 // It illustrates how to do run PYTHIA with LHEF input, allowing a
8 // sample-by-sample generation of
9 // a) Non-matched/non-merged events
10 // b) MLM jet-matched events (kT-MLM, shower-kT, FxFx)
11 // c) CKKW-L and UMEPS-merged events
12 // d) UNLOPS NLO merged events
13 // see the respective sections in the online manual for details.
14 
15 #ifndef Pythia8_aMCatNLOHooks_H
16 #define Pythia8_aMCatNLOHooks_H
17 
18 #include "Pythia8/Pythia.h"
19 
20 namespace Pythia8 {
21 
22 //==========================================================================
23 
24 // Use userhooks to set the number of requested partons dynamically, as
25 // needed when running CKKW-L or UMEPS on a single input file that contains
26 // all parton multiplicities.
27 
29 
30 public:
31 
32  // Constructor and destructor.
33  amcnlo_unitarised_interface() : mergingScheme(0), normFactor(1.) {}
34  amcnlo_unitarised_interface(int mergingSchemeIn)
35  : mergingScheme(mergingSchemeIn), normFactor(1.) {}
37 
38  double getNormFactor(){return normFactor;}
39 
40  // Allow to set the number of partons.
41  bool canVetoProcessLevel() { return true;}
42  // Set the number of partons.
43  bool doVetoProcessLevel( Event& process) {
44 
45  int nPartons = 0;
46  normFactor = 1.;
47 
48  // Do not include resonance decay products in the counting.
49  omitResonanceDecays(process);
50 
51  // Get the maximal quark flavour counted as "additional" parton.
52  int nQuarksMerge = settingsPtr->mode("Merging:nQuarksMerge");
53 
54  // Save information on the process string.
55  string procSave = settingsPtr->word("Merging:Process");
56  bool hasIN = (procSave.find(">", 0) != string::npos);
57  string incoming = procSave.substr(0,procSave.find(">", 0));
58  string outgoing = procSave.substr(procSave.find(">", 0)+1,procSave.size());
59  string outgoingSave = outgoing;
60 
61  // Count number of user particles.
62  int nCommas = 0;
63  for(int n = outgoing.find(",", 0); n != int(string::npos);
64  n = outgoing.find(",", n)) { nCommas++; n++; }
65  vector <string> out;
66  for(int i =0; i < nCommas;++i) {
67  int n = outgoing.find(",", 0);
68  out.push_back(outgoing.substr(0,n));
69  outgoing = outgoing.substr(n+1,outgoing.size());
70  }
71  if (outgoing.size()>0) out.push_back(outgoing);
72 
73  // Dynamically set the process string.
74  string proc;
75  if ( settingsPtr->word("Merging:Process") == "guess" ) {
76  string processString = "";
77  // Set incoming particles.
78  if (hasIN) {
79  processString += incoming + ">";
80  } else {
81  int beamAid = beamAPtr->id();
82  int beamBid = beamBPtr->id();
83  if (abs(beamAid) == 2212) processString += "p";
84  if (beamAid == 11) processString += "e-";
85  if (beamAid ==-11) processString += "e+";
86  if (abs(beamBid) == 2212) processString += "p";
87  if (beamBid == 11) processString += "e-";
88  if (beamBid ==-11) processString += "e+";
89  processString += ">";
90  }
91  // Set outgoing particles.
92  bool foundOutgoing = false;
93  for(int i=0; i < int(workEvent.size()); ++i) {
94  if ( workEvent[i].isFinal()
95  && ( workEvent[i].colType() == 0
96  || workEvent[i].idAbs() > 21
97  || ( workEvent[i].id() != 21
98  && workEvent[i].idAbs() > nQuarksMerge) ) ) {
99  foundOutgoing = true;
100  ostringstream procOSS;
101  procOSS << "{" << workEvent[i].name() << "," << workEvent[i].id()
102  << "}";
103  processString += procOSS.str();
104  }
105  }
106 
107  for (int i =0; i < int(out.size()); ++i) {
108  if (out[i].find("guess") != string::npos) continue;
109  processString += "," + out[i];
110  }
111 
112  if (foundOutgoing) proc = processString;
113  }
114 
115  if (proc.size()>0) settingsPtr->word("Merging:Process", proc);
116  else proc = procSave;
117 
118  // Loop through event and count.
119  for(int i=0; i < int(workEvent.size()); ++i)
120  if ( workEvent[i].isFinal()
121  && workEvent[i].colType()!= 0
122  && ( workEvent[i].id() == 21 || workEvent[i].idAbs() <= nQuarksMerge))
123  nPartons++;
124 
125  // Store merging scheme.
126  bool isumeps = (mergingScheme == 1);
127  bool isunlops = (mergingScheme == 2);
128 
129  // Get number of requested partons.
130  string nps_nlo = infoPtr->getEventAttribute("npNLO",true);
131  int np_nlo = (nps_nlo != "") ? atoi((char*)nps_nlo.c_str()) : -1;
132  string nps_lo = infoPtr->getEventAttribute("npLO",true);
133  int np_lo = (nps_lo != "") ? atoi((char*)nps_lo.c_str()) : -1;
134 
135  if ( (settingsPtr->flag("Merging:doUNLOPSTree")
136  || settingsPtr->flag("Merging:doUNLOPSSubt")) && np_lo == 0)
137  return true;
138 
139  int nj = 0;
140  for(int n = proc.find("j", 0); n != int(string::npos);
141  n = proc.find("j", n)) { nj++; n++; }
142  nPartons -= nj;
143  if (settingsPtr->word("Merging:process").compare("e+e->jj") == 0) {
144  np_lo -= 2;
145  np_nlo -= 2;
146  }
147 
148  // Set number of requested partons.
149  if (np_nlo > -1){
150  settingsPtr->mode("Merging:nRequested", np_nlo);
151  np_lo = -1;
152  } else if (np_lo > -1){
153  settingsPtr->mode("Merging:nRequested", np_lo);
154  np_nlo = -1;
155  } else {
156  settingsPtr->mode("Merging:nRequested", nPartons);
157  np_nlo = -1;
158  np_lo = nPartons;
159  }
160 
161  // Reset the event weight to incorporate corrective factor.
162  bool updateWgt = true;
163 
164  // Choose randomly if this event should be treated as subtraction or
165  // as regular event. Put the correct settings accordingly.
166  if (isunlops && np_nlo == 0 && np_lo == -1) {
167  settingsPtr->flag("Merging:doUNLOPSTree", false);
168  settingsPtr->flag("Merging:doUNLOPSSubt", false);
169  settingsPtr->flag("Merging:doUNLOPSLoop", true);
170  settingsPtr->flag("Merging:doUNLOPSSubtNLO", false);
171  settingsPtr->mode("Merging:nRecluster",0);
172  normFactor *= 1.;
173  } else if (isunlops && np_nlo > 0 && np_lo == -1) {
174  normFactor *= 2.;
175  if (rndmPtr->flat() < 0.5) {
176  normFactor *= -1.;
177  settingsPtr->flag("Merging:doUNLOPSTree", false);
178  settingsPtr->flag("Merging:doUNLOPSSubt", false);
179  settingsPtr->flag("Merging:doUNLOPSLoop", false);
180  settingsPtr->flag("Merging:doUNLOPSSubtNLO", true);
181  settingsPtr->mode("Merging:nRecluster",1);
182  } else {
183  settingsPtr->flag("Merging:doUNLOPSTree", false);
184  settingsPtr->flag("Merging:doUNLOPSSubt", false);
185  settingsPtr->flag("Merging:doUNLOPSLoop", true);
186  settingsPtr->flag("Merging:doUNLOPSSubtNLO", false);
187  settingsPtr->mode("Merging:nRecluster",0);
188  }
189  } else if (isunlops && np_nlo == -1 && np_lo > -1) {
190  normFactor *= 2.;
191  if (rndmPtr->flat() < 0.5) {
192  normFactor *= -1.;
193  settingsPtr->flag("Merging:doUNLOPSTree", false);
194  settingsPtr->flag("Merging:doUNLOPSSubt", true);
195  settingsPtr->flag("Merging:doUNLOPSLoop", false);
196  settingsPtr->flag("Merging:doUNLOPSSubtNLO", false);
197  settingsPtr->mode("Merging:nRecluster",1);
198 
199  // Double reclustering for exclusive NLO cross sections.
200  bool isnlotilde = settingsPtr->flag("Merging:doUNLOPSTilde");
201  int nmaxNLO = settingsPtr->mode("Merging:nJetMaxNLO");
202  if ( isnlotilde
203  && nmaxNLO > 0
204  && np_lo <= nmaxNLO+1
205  && np_lo > 1 ){
206  normFactor *= 2.;
207  if (rndmPtr->flat() < 0.5)
208  settingsPtr->mode("Merging:nRecluster",2);
209  }
210  } else {
211  settingsPtr->flag("Merging:doUNLOPSTree", true);
212  settingsPtr->flag("Merging:doUNLOPSSubt", false);
213  settingsPtr->flag("Merging:doUNLOPSLoop", false);
214  settingsPtr->flag("Merging:doUNLOPSSubtNLO", false);
215  settingsPtr->mode("Merging:nRecluster",0);
216  }
217  } else if (isumeps && np_lo == 0) {
218  settingsPtr->flag("Merging:doUMEPSTree", true);
219  settingsPtr->flag("Merging:doUMEPSSubt", false);
220  settingsPtr->mode("Merging:nRecluster",0);
221  } else if (isumeps && np_lo > 0) {
222  normFactor *= 2.;
223  if (rndmPtr->flat() < 0.5) {
224  normFactor *= -1.;
225  settingsPtr->flag("Merging:doUMEPSTree", false);
226  settingsPtr->flag("Merging:doUMEPSSubt", true);
227  settingsPtr->mode("Merging:nRecluster",1);
228  } else {
229  settingsPtr->flag("Merging:doUMEPSTree", true);
230  settingsPtr->flag("Merging:doUMEPSSubt", false);
231  settingsPtr->mode("Merging:nRecluster",0);
232  }
233  }
234  // Reset the event weight to incorporate corrective factor.
235  if ( updateWgt) {
236  infoPtr->weightContainerPtr->weightNominal *= normFactor;
237  normFactor = 1.;
238  }
239 
240  // Done
241  return false;
242  }
243 
244 private:
245 
246  int mergingScheme;
247  double normFactor;
248 };
249 
250 //==========================================================================
251 
252 } // end namespace Pythia8
253 
254 #endif // end Pythia8_aMCatNLOHooks_H
double weightNominal
The nominal Pythia weight, in pb for lha strategy 4 and -4.
Definition: Weights.h:403
bool flag(string keyIn)
Give back current value, with check that key exists.
Definition: Settings.cc:1478
bool canVetoProcessLevel()
Allow to set the number of partons.
Definition: aMCatNLOHooks.h:41
Settings * settingsPtr
Pointer to the settings database.
Definition: PhysicsBase.h:81
The Event class holds all info on the generated event.
Definition: Event.h:453
Rndm * rndmPtr
Pointer to the random number generator.
Definition: PhysicsBase.h:93
void omitResonanceDecays(const Event &process, bool finalOnly=false)
omitResonanceDecays omits resonance decay chains from process record.
Definition: UserHooks.cc:132
int size() const
Event record size.
Definition: Event.h:504
int id() const
Member functions for output.
Definition: BeamParticle.h:195
bool doVetoProcessLevel(Event &process)
Set the number of partons.
Definition: aMCatNLOHooks.h:43
Header for classes to set beam momentum and interaction vertex spread.
Definition: Analysis.h:20
Event workEvent
Have one event object around as work area.
Definition: UserHooks.h:248
Info * infoPtr
Definition: PhysicsBase.h:78
UserHooks is base class for user access to program execution.
Definition: UserHooks.h:32
double flat()
Generate next random number uniformly between 0 and 1.
Definition: Basics.cc:189
amcnlo_unitarised_interface()
Constructor and destructor.
Definition: aMCatNLOHooks.h:33
string getEventAttribute(string key, bool doRemoveWhitespace=false) const
Retrieve events tag information.
Definition: Info.cc:307
Definition: aMCatNLOHooks.h:28