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) 2018 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 
111  // Set number of requested partons.
112  if (np_nlo > -1){
113  settingsPtr->mode("Merging:nRequested", np_nlo);
114  np_lo = -1;
115  } else if (np_lo > -1){
116  settingsPtr->mode("Merging:nRequested", np_lo);
117  np_nlo = -1;
118  } else {
119  settingsPtr->mode("Merging:nRequested", nPartons);
120  np_nlo = -1;
121  np_lo = nPartons;
122  }
123 
124  // Reset the event weight to incorporate corrective factor.
125  bool updateWgt = settingsPtr->flag("Merging:includeWeightInXsection");
126  double norm = (abs(infoPtr->lhaStrategy()) == 4) ? 1./1e9 : 1.;
127 
128  // Choose randomly if this event should be treated as subtraction or
129  // as regular event. Put the correct settings accordingly.
130  if (isunlops && np_nlo == 0 && np_lo == -1) {
131  settingsPtr->flag("Merging:doUNLOPSTree", false);
132  settingsPtr->flag("Merging:doUNLOPSSubt", false);
133  settingsPtr->flag("Merging:doUNLOPSLoop", true);
134  settingsPtr->flag("Merging:doUNLOPSSubtNLO", false);
135  settingsPtr->mode("Merging:nRecluster",0);
136  normFactor *= 1.;
137  } else if (isunlops && np_nlo > 0 && np_lo == -1) {
138  normFactor *= 2.;
139  if (rndmPtr->flat() < 0.5) {
140  normFactor *= -1.;
141  settingsPtr->flag("Merging:doUNLOPSTree", false);
142  settingsPtr->flag("Merging:doUNLOPSSubt", false);
143  settingsPtr->flag("Merging:doUNLOPSLoop", false);
144  settingsPtr->flag("Merging:doUNLOPSSubtNLO", true);
145  settingsPtr->mode("Merging:nRecluster",1);
146  } else {
147  settingsPtr->flag("Merging:doUNLOPSTree", false);
148  settingsPtr->flag("Merging:doUNLOPSSubt", false);
149  settingsPtr->flag("Merging:doUNLOPSLoop", true);
150  settingsPtr->flag("Merging:doUNLOPSSubtNLO", false);
151  settingsPtr->mode("Merging:nRecluster",0);
152  }
153  } else if (isunlops && np_nlo == -1 && np_lo > -1) {
154  normFactor *= 2.;
155  if (rndmPtr->flat() < 0.5) {
156  normFactor *= -1.;
157  settingsPtr->flag("Merging:doUNLOPSTree", false);
158  settingsPtr->flag("Merging:doUNLOPSSubt", true);
159  settingsPtr->flag("Merging:doUNLOPSLoop", false);
160  settingsPtr->flag("Merging:doUNLOPSSubtNLO", false);
161  settingsPtr->mode("Merging:nRecluster",1);
162 
163  // Double reclustering for exclusive NLO cross sections.
164  bool isnlotilde = settingsPtr->flag("Merging:doUNLOPSTilde");
165  int nmaxNLO = settingsPtr->mode("Merging:nJetMaxNLO");
166  if ( isnlotilde
167  && nmaxNLO > 0
168  && np_lo <= nmaxNLO+1
169  && np_lo > 1 ){
170  normFactor *= 2.;
171  if (rndmPtr->flat() < 0.5)
172  settingsPtr->mode("Merging:nRecluster",2);
173  }
174  } else {
175  settingsPtr->flag("Merging:doUNLOPSTree", true);
176  settingsPtr->flag("Merging:doUNLOPSSubt", false);
177  settingsPtr->flag("Merging:doUNLOPSLoop", false);
178  settingsPtr->flag("Merging:doUNLOPSSubtNLO", false);
179  settingsPtr->mode("Merging:nRecluster",0);
180  }
181  } else if (isumeps && np_lo == 0) {
182  settingsPtr->flag("Merging:doUMEPSTree", true);
183  settingsPtr->flag("Merging:doUMEPSSubt", false);
184  settingsPtr->mode("Merging:nRecluster",0);
185  } else if (isumeps && np_lo > 0) {
186  normFactor *= 2.;
187  if (rndmPtr->flat() < 0.5) {
188  normFactor *= -1.;
189  settingsPtr->flag("Merging:doUMEPSTree", false);
190  settingsPtr->flag("Merging:doUMEPSSubt", true);
191  settingsPtr->mode("Merging:nRecluster",1);
192  } else {
193  settingsPtr->flag("Merging:doUMEPSTree", true);
194  settingsPtr->flag("Merging:doUMEPSSubt", false);
195  settingsPtr->mode("Merging:nRecluster",0);
196  }
197  // Reset the event weight to incorporate corrective factor.
198  if ( updateWgt) {
199  infoPtr->updateWeight(infoPtr->weight()*norm*normFactor);
200  normFactor = 1.;
201  }
202  }
203 
204  // Done
205  return false;
206  }
207 
208 private:
209 
210  int mergingScheme;
211  double normFactor;
212 };
213 
214 //==========================================================================
215 
216 } // end namespace Pythia8
217 
218 #endif // end Pythia8_aMCatNLOHooks_H
Pythia8::amcnlo_unitarised_interface
Definition: aMCatNLOHooks.h:28
Generated by   doxygen 1.8.5