StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
main81.cc
1 // main81.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2014 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 "Pythia8/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 }