StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
main86.cc
1 // main86.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2011 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 UMEPS merging,
8 // see the Matrix Element Merging page in the online manual.
9 
10 #include "Pythia8/Pythia.h"
11 #include "Pythia8/Pythia8ToHepMC.h"
12 #include <unistd.h>
13 
14 #include "HepMC/GenEvent.h"
15 #include "HepMC/IO_GenEvent.h"
16 // Following line to be used with HepMC 2.04 onwards.
17 #include "HepMC/Units.h"
18 
19 using namespace Pythia8;
20 
21 //==========================================================================
22 
23 // Example main programm to illustrate merging
24 
25 int main( int argc, char* argv[] ){
26 
27  // Check that correct number of command-line arguments
28  if (argc != 4) {
29  cerr << " Unexpected number of command-line arguments ("<<argc<<"). \n"
30  << " You are expected to provide the arguments" << endl
31  << " 1. Input file for settings" << endl
32  << " 2. Name of the input LHE file (with path), up to the '_tree'"
33  << " identifier" << endl
34  << " 3. Path for output histogram files" << endl
35  << " Program stopped. " << endl;
36  return 1;
37  }
38 
39  Pythia pythia;
40 
41  // Input parameters:
42  pythia.readFile(argv[1]);
43  // Interface for conversion from Pythia8::Event to HepMC one.
44  HepMC::Pythia8ToHepMC ToHepMC;
45  // Specify file where HepMC events will be stored.
46  HepMC::IO_GenEvent ascii_io(argv[3], std::ios::out);
47  // Switch off warnings for parton-level events.
48  ToHepMC.set_print_inconsistency(false);
49  ToHepMC.set_free_parton_warnings(false);
50  // Do not store cross section information, as this will be done manually.
51  ToHepMC.set_store_pdf(false);
52  ToHepMC.set_store_proc(false);
53  ToHepMC.set_store_xsec(false);
54 
55  // Path to input events, with name up to the "_tree" identifier included.
56  string iPath = string(argv[2]);
57 
58  // Number of events
59  int nEvent = pythia.mode("Main:numberOfEvents");
60  // Maximal number of additional LO jets.
61  int nMaxLO = pythia.mode("Merging:nJetMax");
62 
63 //-----------------------------------------------------------------------------
64 //-----------------------------------------------------------------------------
65 
66  // Switch off all showering and MPI when extimating the cross section after
67  // the merging scale cut.
68  bool fsr = pythia.flag("PartonLevel:FSR");
69  bool isr = pythia.flag("PartonLevel:ISR");
70  bool mpi = pythia.flag("PartonLevel:MPI");
71  bool had = pythia.flag("HadronLevel:all");
72  pythia.settings.flag("PartonLevel:FSR",false);
73  pythia.settings.flag("PartonLevel:ISR",false);
74  pythia.settings.flag("HadronLevel:all",false);
75  pythia.settings.flag("PartonLevel:MPI",false);
76 
77  // Switch on cross section estimation procedure.
78  pythia.settings.flag("Merging:doXSectionEstimate", true);
79  pythia.settings.flag("Merging:doUMEPSTree",true);
80 
81  int njetcounterLO = nMaxLO;
82  string iPathTree = iPath + "_tree";
83 
84  // Save estimates in vectors.
85  vector<double> xsecLO;
86  vector<double> nAcceptLO;
87 
88  cout << endl << endl << endl;
89  cout << "Start estimating umeps tree level cross section" << endl;
90 
91  while(njetcounterLO >= 0) {
92 
93  // From njet, choose LHE file
94  stringstream in;
95  in << "_" << njetcounterLO << ".lhe";
96 #ifdef GZIPSUPPORT
97  if(access( (iPathTree+in.str()+".gz").c_str(), F_OK) != -1) in << ".gz";
98 #endif
99  string LHEfile = iPathTree + in.str();
100 
101  LHAupLHEF lhareader((char*)(LHEfile).c_str());
102  pythia.settings.mode("Merging:nRequested", njetcounterLO);
103  pythia.settings.word("Beams:LHEF", LHEfile);
104  pythia.init(&lhareader);
105 
106  // Start generation loop
107  for( int iEvent=0; iEvent<nEvent; ++iEvent ){
108  // Generate next event
109  if( !pythia.next() ) {
110  if( pythia.info.atEndOfFile() ){
111  break;
112  }
113  else continue;
114  }
115  } // end loop over events to generate
116 
117  // print cross section, errors
118  pythia.stat();
119 
120  xsecLO.push_back(pythia.info.sigmaGen());
121  nAcceptLO.push_back(pythia.info.nAccepted());
122 
123  // Restart with ME of a reduced the number of jets
124  if( njetcounterLO > 0 )
125  njetcounterLO--;
126  else
127  break;
128 
129  } // end loop over different jet multiplicities
130 
131 
132  // Switch off cross section estimation.
133  pythia.settings.flag("Merging:doXSectionEstimate", false);
134 
135  // Switch showering and multiple interaction back on.
136  pythia.settings.flag("PartonLevel:FSR",fsr);
137  pythia.settings.flag("PartonLevel:ISR",isr);
138  pythia.settings.flag("HadronLevel:all",had);
139  pythia.settings.flag("PartonLevel:MPI",mpi);
140 
141  // Declare sample cross section for output.
142  double sigmaTemp = 0.;
143  vector<double> sampleXStree;
144  vector<double> sampleXSsubtTree;
145  // Cross section an error.
146  double sigmaTotal = 0.;
147  double errorTotal = 0.;
148 
149  int sizeLO = int(xsecLO.size());
150  njetcounterLO = nMaxLO;
151  iPathTree = iPath + "_tree";
152 
153  while(njetcounterLO >= 0){
154 
155  // From njet, choose LHE file
156  stringstream in;
157  in << "_" << njetcounterLO << ".lhe";
158 #ifdef GZIPSUPPORT
159  if(access( (iPathTree+in.str()+".gz").c_str(), F_OK) != -1) in << ".gz";
160 #endif
161  string LHEfile = iPathTree + in.str();
162 
163  pythia.settings.flag("Merging:doUMEPSTree",true);
164  pythia.settings.flag("Merging:doUMEPSSubt",false);
165  pythia.settings.mode("Merging:nRecluster",0);
166  LHAupLHEF lhareader((char*)(LHEfile).c_str());
167 
168  cout << endl << endl << endl
169  << "Start tree level treatment for " << njetcounterLO << " jets"
170  << endl;
171 
172  pythia.settings.mode("Merging:nRequested", njetcounterLO);
173  pythia.settings.word("Beams:LHEF", LHEfile);
174  pythia.init(&lhareader);
175 
176  // Remember position in vector of cross section estimates.
177  int iNow = sizeLO-1-njetcounterLO;
178 
179  // Start generation loop
180  for( int iEvent=0; iEvent<nEvent; ++iEvent ){
181 
182  // Generate next event
183  if( !pythia.next() ) {
184  if( pythia.info.atEndOfFile() ) break;
185  else continue;
186  }
187 
188  // Get event weight(s).
189  double weight = pythia.info.mergingWeight();
190  double evtweight = pythia.info.weight();
191  weight *= evtweight;
192  // Do not print zero-weight events.
193  if ( weight == 0. ) continue;
194 
195  // Construct new empty HepMC event.
196  HepMC::GenEvent* hepmcevt = new HepMC::GenEvent();
197  // Get correct cross section from previous estimate.
198  double normhepmc = xsecLO[iNow] / nAcceptLO[iNow];
199  // Set event weight
200  hepmcevt->weights().push_back(weight*normhepmc);
201  // Fill HepMC event.
202  ToHepMC.fill_next_event( pythia, hepmcevt );
203  // Add the weight of the current event to the cross section.
204  sigmaTotal += weight*normhepmc;
205  sigmaTemp += weight*normhepmc;
206  errorTotal += pow2(weight*normhepmc);
207  // Report cross section to hepmc
208  HepMC::GenCrossSection xsec;
209  xsec.set_cross_section( sigmaTotal*1e9, pythia.info.sigmaErr()*1e9 );
210  hepmcevt->set_cross_section( xsec );
211  // Write the HepMC event to file. Done with it.
212  ascii_io << hepmcevt;
213  delete hepmcevt;
214 
215  } // end loop over events to generate
216 
217  // print cross section, errors
218  pythia.stat();
219  // Save sample cross section for output.
220  sampleXStree.push_back(sigmaTemp);
221  sigmaTemp = 0.;
222 
223  // Restart with ME of a reduced the number of jets
224  if( njetcounterLO > 0 )
225  njetcounterLO--;
226  else
227  break;
228 
229  }
230 
231  cout << endl << endl << endl;
232  cout << "Do UMEPS subtraction" << endl;
233 
234  int njetcounterLS = nMaxLO;
235  string iPathSubt = iPath + "_tree";
236 
237  while(njetcounterLS >= 1){
238 
239  // From njet, choose LHE file
240  stringstream in;
241  in << "_" << njetcounterLS << ".lhe";
242 #ifdef GZIPSUPPORT
243  if(access( (iPathSubt+in.str()+".gz").c_str(), F_OK) != -1) in << ".gz";
244 #endif
245  string LHEfile = iPathSubt + in.str();
246 
247  pythia.settings.flag("Merging:doUMEPSTree",false);
248  pythia.settings.flag("Merging:doUMEPSSubt",true);
249  pythia.settings.mode("Merging:nRecluster",1);
250  LHAupLHEF lhareader((char*)(LHEfile).c_str());
251 
252  cout << endl << endl << endl
253  << "Start subtractive treatment for " << njetcounterLS << " jets"
254  << endl;
255 
256  pythia.settings.mode("Merging:nRequested", njetcounterLS);
257  pythia.settings.word("Beams:LHEF", LHEfile);
258  pythia.init(&lhareader);
259 
260  // Remember position in vector of cross section estimates.
261  int iNow = sizeLO-1-njetcounterLS;
262 
263  // Start generation loop
264  for( int iEvent=0; iEvent<nEvent; ++iEvent ){
265 
266  // Generate next event
267  if( !pythia.next() ) {
268  if( pythia.info.atEndOfFile() ) break;
269  else continue;
270  }
271 
272  // Get event weight(s).
273  double weight = pythia.info.mergingWeight();
274  double evtweight = pythia.info.weight();
275  weight *= evtweight;
276  // Do not print zero-weight events.
277  if ( weight == 0. ) continue;
278 
279  // Construct new empty HepMC event.
280  HepMC::GenEvent* hepmcevt = new HepMC::GenEvent();
281  // Get correct cross section from previous estimate.
282  double normhepmc = -1*xsecLO[iNow] / nAcceptLO[iNow];
283  // Set event weight
284  hepmcevt->weights().push_back(weight*normhepmc);
285  // Fill HepMC event.
286  ToHepMC.fill_next_event( pythia, hepmcevt );
287  // Add the weight of the current event to the cross section.
288  sigmaTotal += weight*normhepmc;
289  sigmaTemp += weight*normhepmc;
290  errorTotal += pow2(weight*normhepmc);
291  // Report cross section to hepmc.
292  HepMC::GenCrossSection xsec;
293  xsec.set_cross_section( sigmaTotal*1e9, pythia.info.sigmaErr()*1e9 );
294  hepmcevt->set_cross_section( xsec );
295  // Write the HepMC event to file. Done with it.
296  ascii_io << hepmcevt;
297  delete hepmcevt;
298 
299  } // end loop over events to generate
300 
301  // print cross section, errors
302  pythia.stat();
303  // Save sample cross section for output.
304  sampleXSsubtTree.push_back(sigmaTemp);
305  sigmaTemp = 0.;
306 
307  // Restart with ME of a reduced the number of jets
308  if( njetcounterLS > 1 )
309  njetcounterLS--;
310  else
311  break;
312  }
313 
314  // Print cross section information.
315  cout << endl << endl;
316  cout << " *---------------------------------------------------*" << endl;
317  cout << " | |" << endl;
318  cout << " | Sample cross sections after UMEPS merging |" << endl;
319  cout << " | |" << endl;
320  cout << " | Leading order cross sections (mb): |" << endl;
321  for (int i = 0; i < int(sampleXStree.size()); ++i)
322  cout << " | " << sampleXStree.size()-1-i << "-jet: "
323  << setw(17) << scientific << setprecision(6)
324  << sampleXStree[i] << " |" << endl;
325  cout << " | |" << endl;
326  cout << " | Leading-order subtractive cross sections (mb): |" << endl;
327  for (int i = 0; i < int(sampleXSsubtTree.size()); ++i)
328  cout << " | " << sampleXSsubtTree.size()-1-i+1 << "-jet: "
329  << setw(17) << scientific << setprecision(6)
330  << sampleXSsubtTree[i] << " |" << endl;
331  cout << " | |" << endl;
332  cout << " |---------------------------------------------------|" << endl;
333  cout << " |---------------------------------------------------|" << endl;
334  cout << " | Inclusive cross sections: |" << endl;
335  cout << " | |" << endl;
336  cout << " | UMEPS merged inclusive cross section: |" << endl;
337  cout << " | " << setw(17) << scientific << setprecision(6)
338  << sigmaTotal << " +- " << setw(17) << sqrt(errorTotal) << " mb "
339  << " |" << endl;
340  cout << " | |" << endl;
341  cout << " | LO inclusive cross section: |" << endl;
342  cout << " | " << setw(17) << scientific << setprecision(6)
343  << xsecLO.back() << " mb |" << endl;
344  cout << " | |" << endl;
345  cout << " *---------------------------------------------------*" << endl;
346  cout << endl << endl;
347 
348  // Done
349  return 0;
350 
351 }
HepMC::GenCrossSection::set_cross_section
void set_cross_section(double xs, double xs_err)
Set cross section and error in pb.
Definition: GenCrossSection.h:98
HepMC::GenCrossSection
The GenCrossSection class stores the generated cross section.
Definition: GenCrossSection.h:32
HepMC::GenEvent::set_cross_section
void set_cross_section(const GenCrossSection &)
provide a pointer to the GenCrossSection container
Definition: GenEvent.h:752
HepMC::GenEvent
The GenEvent class is the core of HepMC.
Definition: GenEvent.h:155
HepMC::IO_GenEvent
IO_GenEvent also deals with HeavyIon and PdfInfo.
Definition: IO_GenEvent.h:63
HepMC::GenEvent::weights
WeightContainer & weights()
direct access to WeightContainer
Definition: GenEvent.h:699
Generated by   doxygen 1.8.5