StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
aMCatNLOHooks.h
1 // aMCatNLOHooks.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2020 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 
28 class amcnlo_unitarised_interface : public UserHooks {
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.) {}
36  ~amcnlo_unitarised_interface() {}
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  // Dynamically set the process string.
55  if ( settingsPtr->word("Merging:Process") == "guess" ) {
56  string processString = "";
57  // Set incoming particles.
58  int beamAid = beamAPtr->id();
59  int beamBid = beamBPtr->id();
60  if (abs(beamAid) == 2212) processString += "p";
61  if (beamAid == 11) processString += "e-";
62  if (beamAid ==-11) processString += "e+";
63  if (abs(beamBid) == 2212) processString += "p";
64  if (beamBid == 11) processString += "e-";
65  if (beamBid ==-11) processString += "e+";
66  processString += ">";
67  // Set outgoing particles.
68  bool foundOutgoing = false;
69  for(int i=0; i < int(workEvent.size()); ++i)
70  if ( workEvent[i].isFinal()
71  && ( workEvent[i].colType() == 0
72  || workEvent[i].idAbs() > 21
73  || ( workEvent[i].id() != 21
74  && workEvent[i].idAbs() > nQuarksMerge) ) ) {
75  foundOutgoing = true;
76  ostringstream proc;
77  proc << "{" << workEvent[i].name() << "," << workEvent[i].id()
78  << "}";
79  processString += proc.str();
80  }
81  // Set the process string.
82  if (foundOutgoing) settingsPtr->word("Merging:Process", processString);
83  }
84 
85  // Loop through event and count.
86  for(int i=0; i < int(workEvent.size()); ++i)
87  if ( workEvent[i].isFinal()
88  && workEvent[i].colType()!= 0
89  && ( workEvent[i].id() == 21 || workEvent[i].idAbs() <= nQuarksMerge))
90  nPartons++;
91 
92  // Store merging scheme.
93  bool isumeps = (mergingScheme == 1);
94  bool isunlops = (mergingScheme == 2);
95 
96  // Get number of requested partons.
97  string nps_nlo = infoPtr->getEventAttribute("npNLO",true);
98  int np_nlo = (nps_nlo != "") ? atoi((char*)nps_nlo.c_str()) : -1;
99  string nps_lo = infoPtr->getEventAttribute("npLO",true);
100  int np_lo = (nps_lo != "") ? atoi((char*)nps_lo.c_str()) : -1;
101 
102  if ( (settingsPtr->flag("Merging:doUNLOPSTree")
103  || settingsPtr->flag("Merging:doUNLOPSSubt")) && np_lo == 0)
104  return true;
105 
106  if (settingsPtr->word("Merging:process").compare("pp>aj") == 0)
107  nPartons -= 1;
108  if (settingsPtr->word("Merging:process").compare("pp>jj") == 0)
109  nPartons -= 2;
110  if (settingsPtr->word("Merging:process").compare("e+e->jj") == 0) {
111  nPartons -= 2;
112  np_lo -= 2;
113  np_nlo -= 2;
114  }
115 
116  // Set number of requested partons.
117  if (np_nlo > -1){
118  settingsPtr->mode("Merging:nRequested", np_nlo);
119  np_lo = -1;
120  } else if (np_lo > -1){
121  settingsPtr->mode("Merging:nRequested", np_lo);
122  np_nlo = -1;
123  } else {
124  settingsPtr->mode("Merging:nRequested", nPartons);
125  np_nlo = -1;
126  np_lo = nPartons;
127  }
128 
129  // Reset the event weight to incorporate corrective factor.
130  bool updateWgt = true;
131 
132  // Choose randomly if this event should be treated as subtraction or
133  // as regular event. Put the correct settings accordingly.
134  if (isunlops && np_nlo == 0 && np_lo == -1) {
135  settingsPtr->flag("Merging:doUNLOPSTree", false);
136  settingsPtr->flag("Merging:doUNLOPSSubt", false);
137  settingsPtr->flag("Merging:doUNLOPSLoop", true);
138  settingsPtr->flag("Merging:doUNLOPSSubtNLO", false);
139  settingsPtr->mode("Merging:nRecluster",0);
140  normFactor *= 1.;
141  } else if (isunlops && np_nlo > 0 && np_lo == -1) {
142  normFactor *= 2.;
143  if (rndmPtr->flat() < 0.5) {
144  normFactor *= -1.;
145  settingsPtr->flag("Merging:doUNLOPSTree", false);
146  settingsPtr->flag("Merging:doUNLOPSSubt", false);
147  settingsPtr->flag("Merging:doUNLOPSLoop", false);
148  settingsPtr->flag("Merging:doUNLOPSSubtNLO", true);
149  settingsPtr->mode("Merging:nRecluster",1);
150  } else {
151  settingsPtr->flag("Merging:doUNLOPSTree", false);
152  settingsPtr->flag("Merging:doUNLOPSSubt", false);
153  settingsPtr->flag("Merging:doUNLOPSLoop", true);
154  settingsPtr->flag("Merging:doUNLOPSSubtNLO", false);
155  settingsPtr->mode("Merging:nRecluster",0);
156  }
157  } else if (isunlops && np_nlo == -1 && np_lo > -1) {
158  normFactor *= 2.;
159  if (rndmPtr->flat() < 0.5) {
160  normFactor *= -1.;
161  settingsPtr->flag("Merging:doUNLOPSTree", false);
162  settingsPtr->flag("Merging:doUNLOPSSubt", true);
163  settingsPtr->flag("Merging:doUNLOPSLoop", false);
164  settingsPtr->flag("Merging:doUNLOPSSubtNLO", false);
165  settingsPtr->mode("Merging:nRecluster",1);
166 
167  // Double reclustering for exclusive NLO cross sections.
168  bool isnlotilde = settingsPtr->flag("Merging:doUNLOPSTilde");
169  int nmaxNLO = settingsPtr->mode("Merging:nJetMaxNLO");
170  if ( isnlotilde
171  && nmaxNLO > 0
172  && np_lo <= nmaxNLO+1
173  && np_lo > 1 ){
174  normFactor *= 2.;
175  if (rndmPtr->flat() < 0.5)
176  settingsPtr->mode("Merging:nRecluster",2);
177  }
178  } else {
179  settingsPtr->flag("Merging:doUNLOPSTree", true);
180  settingsPtr->flag("Merging:doUNLOPSSubt", false);
181  settingsPtr->flag("Merging:doUNLOPSLoop", false);
182  settingsPtr->flag("Merging:doUNLOPSSubtNLO", false);
183  settingsPtr->mode("Merging:nRecluster",0);
184  }
185  } else if (isumeps && np_lo == 0) {
186  settingsPtr->flag("Merging:doUMEPSTree", true);
187  settingsPtr->flag("Merging:doUMEPSSubt", false);
188  settingsPtr->mode("Merging:nRecluster",0);
189  } else if (isumeps && np_lo > 0) {
190  normFactor *= 2.;
191  if (rndmPtr->flat() < 0.5) {
192  normFactor *= -1.;
193  settingsPtr->flag("Merging:doUMEPSTree", false);
194  settingsPtr->flag("Merging:doUMEPSSubt", true);
195  settingsPtr->mode("Merging:nRecluster",1);
196  } else {
197  settingsPtr->flag("Merging:doUMEPSTree", true);
198  settingsPtr->flag("Merging:doUMEPSSubt", false);
199  settingsPtr->mode("Merging:nRecluster",0);
200  }
201  }
202  // Reset the event weight to incorporate corrective factor.
203  if ( updateWgt) {
204  infoPtr->weightContainerPtr->weightNominal *= normFactor;
205  normFactor = 1.;
206  }
207 
208  // Done
209  return false;
210  }
211 
212 private:
213 
214  int mergingScheme;
215  double normFactor;
216 };
217 
218 //==========================================================================
219 
220 } // end namespace Pythia8
221 
222 #endif // end Pythia8_aMCatNLOHooks_H
Event
Definition: AgUStep.h:26
Generated by   doxygen 1.8.5