1 | // main81.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 the 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 | // Example main programm to illustrate merging |
---|
82 | |
---|
83 | int main( int argc, char* argv[] ){ |
---|
84 | |
---|
85 | // Check that correct number of command-line arguments |
---|
86 | if (argc != 4) { |
---|
87 | cerr << " Unexpected number of command-line arguments ("<<argc<<"). \n" |
---|
88 | << " You are expected to provide the arguments" << endl |
---|
89 | << " 1. Input file for settings" << endl |
---|
90 | << " 2. Full name of the input LHE file (with path)" << endl |
---|
91 | << " 3. Path for output histogram files" << endl |
---|
92 | << " Program stopped. " << endl; |
---|
93 | return 1; |
---|
94 | } |
---|
95 | |
---|
96 | Pythia pythia; |
---|
97 | |
---|
98 | // Input parameters: |
---|
99 | // 1. Input file for settings |
---|
100 | // 2. Path to input LHE file |
---|
101 | // 3. OUtput histogram path |
---|
102 | pythia.readFile(argv[1]); |
---|
103 | string iPath = string(argv[2]); |
---|
104 | string oPath = string(argv[3]); |
---|
105 | // Number of events |
---|
106 | int nEvent = pythia.mode("Main:numberOfEvents"); |
---|
107 | |
---|
108 | // For ISR regularisation off |
---|
109 | pythia.settings.forceParm("SpaceShower:pT0Ref",0.); |
---|
110 | |
---|
111 | // Declare histograms |
---|
112 | Hist histPTFirst("pT of first jet",100,0.,100.); |
---|
113 | Hist histPTSecond("pT of second jet",100,0.,100.); |
---|
114 | Hist histPTThird("pT of third jet",100,0.,100.); |
---|
115 | |
---|
116 | // Read in ME configurations |
---|
117 | pythia.init(iPath,false); |
---|
118 | |
---|
119 | if(pythia.flag("Main:showChangedSettings")) { |
---|
120 | pythia.settings.listChanged(); |
---|
121 | } |
---|
122 | |
---|
123 | // Start generation loop |
---|
124 | for( int iEvent=0; iEvent<nEvent; ++iEvent ){ |
---|
125 | |
---|
126 | // Generate next event |
---|
127 | if( ! pythia.next()) continue; |
---|
128 | |
---|
129 | // Get CKKWL weight of current event |
---|
130 | double weight = pythia.info.mergingWeight(); |
---|
131 | |
---|
132 | // Fill bins with CKKWL weight |
---|
133 | // Functions use fastjet to get first / second jet |
---|
134 | double pTfirst = pTfirstJet(pythia.event,1, 0.4); |
---|
135 | double pTsecnd = pTfirstJet(pythia.event,2, 0.4); |
---|
136 | double pTthird = pTfirstJet(pythia.event,3, 0.4); |
---|
137 | histPTFirst.fill( pTfirst, weight); |
---|
138 | histPTSecond.fill( pTsecnd, weight); |
---|
139 | histPTThird.fill( pTthird, weight); |
---|
140 | |
---|
141 | if(iEvent%1000 == 0) cout << iEvent << endl; |
---|
142 | |
---|
143 | } // end loop over events to generate |
---|
144 | |
---|
145 | // print cross section, errors |
---|
146 | pythia.stat(); |
---|
147 | |
---|
148 | // Normalise histograms |
---|
149 | double norm = 1. |
---|
150 | * pythia.info.sigmaGen() |
---|
151 | * 1./ double(nEvent); |
---|
152 | |
---|
153 | histPTFirst *= norm; |
---|
154 | histPTSecond *= norm; |
---|
155 | histPTThird *= norm; |
---|
156 | |
---|
157 | // Get the number of jets in the LHE file from the file name |
---|
158 | string jetsInLHEF = iPath.substr(iPath.size()-5, iPath.size()); |
---|
159 | jetsInLHEF = jetsInLHEF.substr(0, jetsInLHEF.size()-4); |
---|
160 | |
---|
161 | // Write histograms to dat file. Use "jetsInLHEF" to label the files |
---|
162 | // Once all the samples up to the maximal desired jet multiplicity from the |
---|
163 | // matrix element are run, add all histograms to produce a |
---|
164 | // matrix-element-merged prediction |
---|
165 | |
---|
166 | ofstream write; |
---|
167 | stringstream suffix; |
---|
168 | suffix << jetsInLHEF << "_wv.dat"; |
---|
169 | |
---|
170 | // Write histograms to file |
---|
171 | write.open( (char*)(oPath + "PTjet1_" + suffix.str()).c_str()); |
---|
172 | histPTFirst.table(write); |
---|
173 | write.close(); |
---|
174 | |
---|
175 | write.open( (char*)(oPath + "PTjet2_" + suffix.str()).c_str()); |
---|
176 | histPTSecond.table(write); |
---|
177 | write.close(); |
---|
178 | |
---|
179 | write.open( (char*)(oPath + "PTjet3_" + suffix.str()).c_str()); |
---|
180 | histPTThird.table(write); |
---|
181 | write.close(); |
---|
182 | |
---|
183 | // Done |
---|
184 | return 0; |
---|
185 | |
---|
186 | } |
---|