1 | // main83.cc is a part of the PYTHIA event generator. |
---|
2 | // Copyright (C) 2012 Torbjorn Sjostrand. |
---|
3 | // PYTHIA is licenced under the GNU GPL version 2, 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 CKKW-L merging, |
---|
8 | // see the Matrix Element Merging page in the online manual. |
---|
9 | |
---|
10 | #include "Pythia.h" |
---|
11 | |
---|
12 | using namespace Pythia8; |
---|
13 | |
---|
14 | // Functions for histogramming |
---|
15 | #include "fastjet/PseudoJet.hh" |
---|
16 | #include "fastjet/ClusterSequence.hh" |
---|
17 | #include "fastjet/CDFMidPointPlugin.hh" |
---|
18 | #include "fastjet/CDFJetCluPlugin.hh" |
---|
19 | #include "fastjet/D0RunIIConePlugin.hh" |
---|
20 | |
---|
21 | //========================================================================== |
---|
22 | |
---|
23 | // Find the Durham kT separation of the clustering from |
---|
24 | // nJetMin --> nJetMin-1 jets in te input event |
---|
25 | |
---|
26 | double pTfirstJet( const Event& event, int nJetMin, double Rparam) { |
---|
27 | |
---|
28 | double yPartonMax = 4.; |
---|
29 | |
---|
30 | // Fastjet analysis - select algorithm and parameters |
---|
31 | fastjet::Strategy strategy = fastjet::Best; |
---|
32 | fastjet::RecombinationScheme recombScheme = fastjet::E_scheme; |
---|
33 | fastjet::JetDefinition *jetDef = NULL; |
---|
34 | // For hadronic collision, use hadronic Durham kT measure |
---|
35 | if(event[3].colType() != 0 || event[4].colType() != 0) |
---|
36 | jetDef = new fastjet::JetDefinition(fastjet::kt_algorithm, Rparam, |
---|
37 | recombScheme, strategy); |
---|
38 | // For e+e- collision, use e+e- Durham kT measure |
---|
39 | else |
---|
40 | jetDef = new fastjet::JetDefinition(fastjet::ee_kt_algorithm, |
---|
41 | recombScheme, strategy); |
---|
42 | // Fastjet input |
---|
43 | std::vector <fastjet::PseudoJet> fjInputs; |
---|
44 | // Reset Fastjet input |
---|
45 | fjInputs.resize(0); |
---|
46 | |
---|
47 | // Loop over event record to decide what to pass to FastJet |
---|
48 | for (int i = 0; i < event.size(); ++i) { |
---|
49 | // (Final state && coloured+photons) only! |
---|
50 | if ( !event[i].isFinal() |
---|
51 | || event[i].isLepton() |
---|
52 | || event[i].id() == 23 |
---|
53 | || abs(event[i].id()) == 24 |
---|
54 | || abs(event[i].y()) > yPartonMax) |
---|
55 | continue; |
---|
56 | |
---|
57 | // Store as input to Fastjet |
---|
58 | fjInputs.push_back( fastjet::PseudoJet (event[i].px(), |
---|
59 | event[i].py(), event[i].pz(),event[i].e() ) ); |
---|
60 | } |
---|
61 | |
---|
62 | // Do nothing for empty input |
---|
63 | if (int(fjInputs.size()) == 0) { |
---|
64 | delete jetDef; |
---|
65 | return 0.0; |
---|
66 | } |
---|
67 | |
---|
68 | // Run Fastjet algorithm |
---|
69 | fastjet::ClusterSequence clustSeq(fjInputs, *jetDef); |
---|
70 | // Extract kT of first clustering |
---|
71 | double pTFirst = sqrt(clustSeq.exclusive_dmerge_max(nJetMin-1)); |
---|
72 | |
---|
73 | delete jetDef; |
---|
74 | // Return kT |
---|
75 | return pTFirst; |
---|
76 | |
---|
77 | } |
---|
78 | |
---|
79 | //========================================================================== |
---|
80 | |
---|
81 | // Class for user interaction with the merging |
---|
82 | |
---|
83 | class MyMergingHooks : public MergingHooks { |
---|
84 | |
---|
85 | private: |
---|
86 | |
---|
87 | public: |
---|
88 | |
---|
89 | // Default constructor |
---|
90 | MyMergingHooks(); |
---|
91 | // Destructor |
---|
92 | ~MyMergingHooks(); |
---|
93 | |
---|
94 | // Functional definition of the merging scale |
---|
95 | virtual double tmsDefinition( const Event& event); |
---|
96 | |
---|
97 | // Function to dampen weights calculated from histories with lowest |
---|
98 | // multiplicity reclustered events that do not pass the ME cuts |
---|
99 | virtual double dampenIfFailCuts( const Event& inEvent ); |
---|
100 | |
---|
101 | // Helper function for tms definition |
---|
102 | double myKTdurham(const Particle& RadAfterBranch, |
---|
103 | const Particle& EmtAfterBranch, int Type, double D ); |
---|
104 | |
---|
105 | }; |
---|
106 | |
---|
107 | //-------------------------------------------------------------------------- |
---|
108 | |
---|
109 | // Constructor |
---|
110 | MyMergingHooks::MyMergingHooks() {} |
---|
111 | |
---|
112 | // Desctructor |
---|
113 | MyMergingHooks::~MyMergingHooks() {} |
---|
114 | |
---|
115 | //-------------------------------------------------------------------------- |
---|
116 | |
---|
117 | double MyMergingHooks::dampenIfFailCuts( const Event& inEvent ){ |
---|
118 | |
---|
119 | // Get pT for pure QCD 2->2 state |
---|
120 | double pT = 0.; |
---|
121 | for( int i=0; i < inEvent.size(); ++i) |
---|
122 | if(inEvent[i].isFinal() && inEvent[i].colType() != 0) { |
---|
123 | pT = sqrt(pow(inEvent[i].px(),2) + pow(inEvent[i].py(),2)); |
---|
124 | break; |
---|
125 | } |
---|
126 | |
---|
127 | // Veto history if lowest multiplicity event does not pass ME cuts |
---|
128 | if(pT < 10.) return 0.; |
---|
129 | |
---|
130 | return 1.; |
---|
131 | |
---|
132 | } |
---|
133 | |
---|
134 | //-------------------------------------------------------------------------- |
---|
135 | |
---|
136 | // Definition of the merging scale |
---|
137 | |
---|
138 | double MyMergingHooks::tmsDefinition( const Event& event){ |
---|
139 | |
---|
140 | // Cut only on QCD partons! |
---|
141 | // Count particle types |
---|
142 | int nFinalColoured = 0; |
---|
143 | int nFinalNow =0; |
---|
144 | for( int i=0; i < event.size(); ++i) { |
---|
145 | if(event[i].isFinal()){ |
---|
146 | if(event[i].id() != 23 && abs(event[i].id()) != 24) |
---|
147 | nFinalNow++; |
---|
148 | if( event[i].colType() != 0) |
---|
149 | nFinalColoured++; |
---|
150 | } |
---|
151 | } |
---|
152 | |
---|
153 | // Use MergingHooks in-built functions to get information on the hard process |
---|
154 | int nLeptons = nHardOutLeptons(); |
---|
155 | int nQuarks = nHardOutPartons(); |
---|
156 | int nResNow = nResInCurrent(); |
---|
157 | |
---|
158 | // Check if photons, electrons etc. have been produced. If so, do not veto |
---|
159 | if(nFinalNow - ( (nLeptons+nQuarks)/2 - nResNow)*2 != nFinalColoured){ |
---|
160 | // Sometimes, Pythia detaches the decay products even though no |
---|
161 | // resonance was put into the LHE file, to catch this, add another |
---|
162 | // if statement |
---|
163 | if(nFinalNow != nFinalColoured) return 0.; |
---|
164 | } |
---|
165 | |
---|
166 | // Check that one parton has been produced. If not (e.g. in MPI), do not veto |
---|
167 | int nMPI = infoPtr->nMPI(); |
---|
168 | if(nMPI > 1) return 0.; |
---|
169 | |
---|
170 | // Declare kT algorithm parameters |
---|
171 | double Dparam = 0.4; |
---|
172 | int kTtype = -1; |
---|
173 | // Declare final parton vector |
---|
174 | vector <int> FinalPartPos; |
---|
175 | FinalPartPos.clear(); |
---|
176 | // Search event record for final state partons |
---|
177 | for (int i=0; i < event.size(); ++i) |
---|
178 | if(event[i].isFinal() && event[i].colType() != 0) |
---|
179 | FinalPartPos.push_back(i); |
---|
180 | |
---|
181 | // Find minimal Durham kT in event, using own function: Check |
---|
182 | // definition of separation |
---|
183 | int type = (event[3].colType() == 0 && event[4].colType() == 0) ? 1 : kTtype; |
---|
184 | // Find minimal kT |
---|
185 | double ktmin = event[0].e(); |
---|
186 | for(int i=0; i < int(FinalPartPos.size()); ++i){ |
---|
187 | double kt12 = ktmin; |
---|
188 | // Compute separation to the beam axis for hadronic collisions |
---|
189 | if(type == -1 || type == -2) { |
---|
190 | double temp = event[FinalPartPos[i]].pT(); |
---|
191 | kt12 = min(kt12, temp); |
---|
192 | } |
---|
193 | // Compute separation to other final state jets |
---|
194 | for(int j=i+1; j < int(FinalPartPos.size()); ++j) { |
---|
195 | double temp = kTdurham( event[FinalPartPos[i]], event[FinalPartPos[j]], |
---|
196 | type, Dparam); |
---|
197 | kt12 = min(kt12, temp); |
---|
198 | } |
---|
199 | // Keep the minimal Durham separation |
---|
200 | ktmin = min(ktmin,kt12); |
---|
201 | } |
---|
202 | |
---|
203 | // Return minimal Durham kT |
---|
204 | return ktmin; |
---|
205 | |
---|
206 | } |
---|
207 | |
---|
208 | //-------------------------------------------------------------------------- |
---|
209 | |
---|
210 | // Function to compute durham y separation from Particle input |
---|
211 | |
---|
212 | double MyMergingHooks::myKTdurham(const Particle& RadAfterBranch, |
---|
213 | const Particle& EmtAfterBranch, int Type, double D ){ |
---|
214 | |
---|
215 | // Declare return variable |
---|
216 | double ktdur; |
---|
217 | // Save 4-momenta of final state particles |
---|
218 | Vec4 jet1 = RadAfterBranch.p(); |
---|
219 | Vec4 jet2 = EmtAfterBranch.p(); |
---|
220 | |
---|
221 | if( Type == 1) { |
---|
222 | // Get angle between jets for e+e- collisions, make sure that |
---|
223 | // -1 <= cos(theta) <= 1 |
---|
224 | double costh; |
---|
225 | if (jet1.pAbs()*jet2.pAbs() <=0.) costh = 1.; |
---|
226 | else { |
---|
227 | costh = costheta(jet1,jet2); |
---|
228 | } |
---|
229 | // Calculate kt durham separation between jets for e+e- collisions |
---|
230 | ktdur = 2.0*min( pow(jet1.e(),2) , (pow(jet2.e(),2)) )*(1.0 - costh); |
---|
231 | } else if( Type == -1 ){ |
---|
232 | // Get delta_eta and cosh(Delta_eta) for hadronic collisions |
---|
233 | double eta1 = 0.5*log( (jet1.e() + jet1.pz()) / (jet1.e() - jet1.pz()) ); |
---|
234 | double eta2 = 0.5*log( (jet2.e() + jet2.pz()) / (jet2.e() - jet2.pz()) ); |
---|
235 | // Get delta_phi and cos(Delta_phi) for hadronic collisions |
---|
236 | double pt1 = sqrt( pow(jet1.px(),2) + pow(jet1.py(),2) ); |
---|
237 | double pt2 = sqrt( pow(jet2.px(),2) + pow(jet2.py(),2) ); |
---|
238 | double cosdPhi = ( jet1.px()*jet2.px() + jet1.py()*jet2.py() ) / (pt1*pt2); |
---|
239 | double dPhi = acos( cosdPhi ); |
---|
240 | // Calculate kT durham like fastjet |
---|
241 | ktdur = min( pow(pt1,2),pow(pt2,2) ) |
---|
242 | * ( pow(eta1-eta2,2) + pow(dPhi,2) ) / pow(D,2); |
---|
243 | } else if( Type == -2 ){ |
---|
244 | // Get delta_eta and cosh(Delta_eta) for hadronic collisions |
---|
245 | double eta1 = 0.5*log( (jet1.e() + jet1.pz()) / (jet1.e() - jet1.pz()) ); |
---|
246 | double eta2 = 0.5*log( (jet2.e() + jet2.pz()) / (jet2.e() - jet2.pz()) ); |
---|
247 | double coshdEta = cosh( eta1 - eta2 ); |
---|
248 | // Get delta_phi and cos(Delta_phi) for hadronic collisions |
---|
249 | double pt1 = sqrt( pow(jet1.px(),2) + pow(jet1.py(),2) ); |
---|
250 | double pt2 = sqrt( pow(jet2.px(),2) + pow(jet2.py(),2) ); |
---|
251 | double cosdPhi = ( jet1.px()*jet2.px() + jet1.py()*jet2.py() ) / (pt1*pt2); |
---|
252 | // Calculate kT durham separation "SHERPA-like" |
---|
253 | ktdur = 2.0*min( pow(pt1,2),pow(pt2,2) ) |
---|
254 | * ( coshdEta - cosdPhi ) / pow(D,2); |
---|
255 | } else { |
---|
256 | ktdur = 0.0; |
---|
257 | } |
---|
258 | // Return kT |
---|
259 | return sqrt(ktdur); |
---|
260 | } |
---|
261 | |
---|
262 | //========================================================================== |
---|
263 | |
---|
264 | // Example main programm to illustrate merging |
---|
265 | |
---|
266 | int main( int argc, char* argv[] ){ |
---|
267 | |
---|
268 | // Check that correct number of command-line arguments |
---|
269 | if (argc != 4) { |
---|
270 | cerr << " Unexpected number of command-line arguments. \n You are" |
---|
271 | << " expected to provide the arguments \n" |
---|
272 | << " 1. Input file for settings \n" |
---|
273 | << " 2. Full name of the input LHE file (with path) \n" |
---|
274 | << " 3. Path for output histogram files \n" |
---|
275 | << " Program stopped. " << endl; |
---|
276 | return 1; |
---|
277 | } |
---|
278 | |
---|
279 | Pythia pythia; |
---|
280 | |
---|
281 | // Input parameters: |
---|
282 | // 1. Input file for settings |
---|
283 | // 2. Path to input LHE file |
---|
284 | // 3. OUtput histogram path |
---|
285 | pythia.readFile(argv[1]); |
---|
286 | string iPath = string(argv[2]); |
---|
287 | string oPath = string(argv[3]); |
---|
288 | |
---|
289 | // Number of events |
---|
290 | int nEvent = pythia.mode("Main:numberOfEvents"); |
---|
291 | |
---|
292 | // Construct user inut for merging |
---|
293 | MergingHooks* myMergingHooks = new MyMergingHooks(); |
---|
294 | pythia.setMergingHooksPtr( myMergingHooks ); |
---|
295 | |
---|
296 | // For ISR regularisation off |
---|
297 | pythia.settings.forceParm("SpaceShower:pT0Ref",0.); |
---|
298 | |
---|
299 | // Declare histograms |
---|
300 | Hist histPTFirst("pT of first jet",100,0.,100.); |
---|
301 | Hist histPTSecond("pT of second jet",100,0.,100.); |
---|
302 | |
---|
303 | // Read in ME configurations |
---|
304 | pythia.init(iPath,false); |
---|
305 | |
---|
306 | if(pythia.flag("Main:showChangedSettings")) { |
---|
307 | pythia.settings.listChanged(); |
---|
308 | } |
---|
309 | |
---|
310 | // Start generation loop |
---|
311 | for( int iEvent=0; iEvent<nEvent; ++iEvent ){ |
---|
312 | |
---|
313 | // Generate next event |
---|
314 | if( ! pythia.next()) continue; |
---|
315 | |
---|
316 | // Get CKKWL weight of current event |
---|
317 | double weight = pythia.info.mergingWeight(); |
---|
318 | |
---|
319 | // Fill bins with CKKWL weight |
---|
320 | double pTfirst = pTfirstJet(pythia.event,1, 0.4); |
---|
321 | double pTsecnd = pTfirstJet(pythia.event,2, 0.4); |
---|
322 | histPTFirst.fill( pTfirst, weight); |
---|
323 | histPTSecond.fill( pTsecnd, weight); |
---|
324 | |
---|
325 | if(iEvent%1000 == 0) cout << iEvent << endl; |
---|
326 | |
---|
327 | } // end loop over events to generate |
---|
328 | |
---|
329 | // print cross section, errors |
---|
330 | pythia.stat(); |
---|
331 | |
---|
332 | // Normalise histograms |
---|
333 | double norm = 1. |
---|
334 | * pythia.info.sigmaGen() |
---|
335 | * 1./ double(nEvent); |
---|
336 | |
---|
337 | histPTFirst *= norm; |
---|
338 | histPTSecond *= norm; |
---|
339 | |
---|
340 | // Get the number of jets in the LHE file from the file name |
---|
341 | string jetsInLHEF = iPath.substr(iPath.size()-5, iPath.size()); |
---|
342 | jetsInLHEF = jetsInLHEF.substr(0, jetsInLHEF.size()-4); |
---|
343 | |
---|
344 | // Write histograms to dat file. Use "jetsInLHEF" to label the files |
---|
345 | // Once all the samples up to the maximal desired jet multiplicity from the |
---|
346 | // matrix element are run, add all histograms to produce a |
---|
347 | // matrix-element-merged prediction |
---|
348 | |
---|
349 | ofstream write; |
---|
350 | stringstream suffix; |
---|
351 | suffix << jetsInLHEF << "_wv.dat"; |
---|
352 | |
---|
353 | // Write histograms to file |
---|
354 | write.open( (char*)(oPath + "PTjet1_" + suffix.str()).c_str()); |
---|
355 | histPTFirst.table(write); |
---|
356 | write.close(); |
---|
357 | |
---|
358 | write.open( (char*)(oPath + "PTjet2_" + suffix.str()).c_str()); |
---|
359 | histPTSecond.table(write); |
---|
360 | write.close(); |
---|
361 | |
---|
362 | |
---|
363 | delete myMergingHooks; |
---|
364 | return 0; |
---|
365 | |
---|
366 | // Done |
---|
367 | } |
---|