StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
main87.cc
1 // main87.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 NL3 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  << " or '_powheg' identifiers" << endl
34  << " 3. Output hepmc file name" << endl
35  << " Program stopped. " << endl;
36  return 1;
37  }
38 
39  Pythia pythia;
40 
41  // Input parameters:
42  // 1. Input file for settings
43  // 2. Path to input LHE file
44  // 3. Output histogram path
45  pythia.readFile(argv[1]);
46 
47  // Interface for conversion from Pythia8::Event to HepMC one.
48  HepMC::Pythia8ToHepMC ToHepMC;
49  // Specify file where HepMC events will be stored.
50  HepMC::IO_GenEvent ascii_io(argv[3], std::ios::out);
51  // Switch off warnings for parton-level events.
52  ToHepMC.set_print_inconsistency(false);
53  ToHepMC.set_free_parton_warnings(false);
54  // Do not store cross section information, as this will be done manually.
55  ToHepMC.set_store_pdf(false);
56  ToHepMC.set_store_proc(false);
57  ToHepMC.set_store_xsec(false);
58 
59  // Path to input events, with name up to the "_tree", "_powheg" identifier
60  // included.
61  string iPath = string(argv[2]);
62 
63  // Number of events
64  int nEvent = pythia.mode("Main:numberOfEvents");
65  // Maximal number of additional LO jets.
66  int nMaxLO = pythia.mode("Merging:nJetMax");
67  // maximal number of additional NLO jets.
68  int nMaxNLO = pythia.mode("Merging:nJetMaxNLO");
69 
70 //-----------------------------------------------------------------------------
71 //-----------------------------------------------------------------------------
72 
73  // Switch off all showering and MPI when extimating the cross section after
74  // the merging scale cut.
75  bool fsr = pythia.flag("PartonLevel:FSR");
76  bool isr = pythia.flag("PartonLevel:ISR");
77  bool mpi = pythia.flag("PartonLevel:MPI");
78  bool had = pythia.flag("HadronLevel:all");
79  pythia.settings.flag("PartonLevel:FSR",false);
80  pythia.settings.flag("PartonLevel:ISR",false);
81  pythia.settings.flag("HadronLevel:all",false);
82  pythia.settings.flag("PartonLevel:MPI",false);
83 
84  // Switch on cross section estimation procedure.
85  pythia.settings.flag("Merging:doXSectionEstimate", true);
86  pythia.settings.flag("Merging:doNL3Tree", true);
87 
88  int njetcounterLO = nMaxLO;
89  string iPathTree = iPath + "_tree";
90 
91  vector<double> xsecLO;
92  vector<double> nSelectedLO;
93  vector<double> nAcceptLO;
94 
95  cout << endl << endl << endl;
96  cout << "Start estimating nl3 tree level cross section" << endl;
97 
98  while(njetcounterLO >= 0){
99 
100  // From njetcounter, choose LHE file
101  stringstream in;
102  in << "_" << njetcounterLO << ".lhe";
103 #ifdef GZIPSUPPORT
104  if(access( (iPathTree+in.str()+".gz").c_str(), F_OK) != -1) in << ".gz";
105 #endif
106  string LHEfile = iPathTree + in.str();
107 
108  LHAupLHEF lhareader((char*)(LHEfile).c_str());
109  pythia.settings.mode("Merging:nRequested", njetcounterLO);
110  pythia.settings.word("Beams:LHEF", LHEfile);
111  pythia.init(&lhareader);
112 
113  // Start generation loop
114  for( int iEvent=0; iEvent<nEvent; ++iEvent ){
115  // Generate next event
116  if( !pythia.next() ) {
117  if( pythia.info.atEndOfFile() ){
118  break;
119  }
120  else continue;
121  }
122  } // end loop over events to generate
123 
124  // print cross section, errors
125  pythia.stat();
126 
127  xsecLO.push_back(pythia.info.sigmaGen());
128  nSelectedLO.push_back(pythia.info.nSelected());
129  nAcceptLO.push_back(pythia.info.nAccepted());
130 
131  // Restart with ME of a reduced the number of jets
132  if( njetcounterLO > 0 )
133  njetcounterLO--;
134  else
135  break;
136 
137  }
138 
139 //-----------------------------------------------------------------------------
140 //-----------------------------------------------------------------------------
141 
142  cout << endl << endl << endl;
143  cout << "Start estimating nl3 virtual corrections cross section" << endl;
144 
145  pythia.settings.flag("Merging:doNL3Tree",false);
146  pythia.settings.flag("Merging:doNL3Loop",true);
147 
148  int njetcounterNLO = nMaxNLO;
149  string iPathLoop= iPath + "_powheg";
150 
151  vector<double> xsecNLO;
152  vector<double> nSelectedNLO;
153  vector<double> nAcceptNLO;
154  vector<int> strategyNLO;
155 
156  while(njetcounterNLO >= 0){
157 
158  // From njetcounter, choose LHE file
159  stringstream in;
160  in << "_" << njetcounterNLO << ".lhe";
161 #ifdef GZIPSUPPORT
162  if(access( (iPathLoop+in.str()+".gz").c_str(), F_OK) != -1) in << ".gz";
163 #endif
164  string LHEfile = iPathLoop + in.str();
165 
166  LHAupLHEF lhareader((char*)(LHEfile).c_str());
167  pythia.settings.mode("Merging:nRequested", njetcounterNLO);
168  pythia.settings.word("Beams:LHEF", LHEfile);
169  pythia.init(&lhareader);
170 
171  // Start generation loop
172  for( int iEvent=0; iEvent<nEvent; ++iEvent ){
173  // Generate next event
174  if( !pythia.next() ) {
175  if( pythia.info.atEndOfFile() ){
176  break;
177  }
178  else continue;
179  }
180  } // end loop over events to generate
181 
182  // print cross section, errors
183  pythia.stat();
184 
185  xsecNLO.push_back(pythia.info.sigmaGen());
186  nSelectedNLO.push_back(pythia.info.nSelected());
187  nAcceptNLO.push_back(pythia.info.nAccepted());
188  strategyNLO.push_back(pythia.info.lhaStrategy());
189 
190  // Restart with ME of a reduced the number of jets
191  if( njetcounterNLO > 0 )
192  njetcounterNLO--;
193  else
194  break;
195 
196  }
197 
198  // Set k-factors
199  int sizeLO = int(xsecLO.size());
200  int sizeNLO = int(xsecNLO.size());
201  double k0 = 1.;
202  double k1 = 1.;
203  double k2 = 1.;
204  // Lowest order k-factor only
205  if ( false ) k1 = k2 = k0 = xsecNLO.back() / xsecLO.back();
206  // No k-factors
207  if ( true ) k0 = k1 = k2 = 1.;
208 
209  cout << " K-Factors :" << endl;
210  cout << "k0 = " << k0 << endl;
211  cout << "k1 = " << k1 << endl;
212  cout << "k2 = " << k2 << endl;
213 
214  // Switch off cross section estimation.
215  pythia.settings.flag("Merging:doXSectionEstimate", false);
216 
217  // Switch showering and multiple interaction back on.
218  pythia.settings.flag("PartonLevel:FSR",fsr);
219  pythia.settings.flag("PartonLevel:ISR",isr);
220  pythia.settings.flag("HadronLevel:all",had);
221  pythia.settings.flag("PartonLevel:MPI",mpi);
222 
223 //-----------------------------------------------------------------------------
224 //-----------------------------------------------------------------------------
225 
226  // Declare sample cross section for output.
227  double sigmaTemp = 0.;
228  vector<double> sampleXStree;
229  vector<double> sampleXSvirt;
230  vector<double> sampleXSsubtTree;
231  // Cross section an error.
232  double sigmaTotal = 0.;
233  double errorTotal = 0.;
234 
235  // Switch on tree-level processing.
236  pythia.settings.flag("Merging:doNL3Tree",true);
237  pythia.settings.flag("Merging:doNL3Loop",false);
238  pythia.settings.flag("Merging:doNL3Subt",false);
239  pythia.settings.mode("Merging:nRecluster",0);
240  pythia.settings.mode("Merging:nRequested", -1);
241 
242  njetcounterLO = nMaxLO;
243  iPathTree = iPath + "_tree";
244 
245  while(njetcounterLO >= 0){
246 
247  // Set k factors
248  pythia.settings.parm("Merging:kFactor0j", k0);
249  pythia.settings.parm("Merging:kFactor1j", k1);
250  pythia.settings.parm("Merging:kFactor2j", k2);
251 
252  // From njetcounter, choose LHE file
253  stringstream in;
254  in << "_" << njetcounterLO << ".lhe";
255 #ifdef GZIPSUPPORT
256  if(access( (iPathTree+in.str()+".gz").c_str(), F_OK) != -1) in << ".gz";
257 #endif
258  string LHEfile = iPathTree + in.str();
259 
260  cout << endl << endl << endl
261  << "Start tree level treatment for " << njetcounterLO << " jets"
262  << endl;
263 
264  LHAupLHEF lhareader((char*)(LHEfile).c_str());
265  pythia.settings.mode("Merging:nRequested", njetcounterLO);
266  pythia.settings.word("Beams:LHEF", LHEfile);
267  pythia.init(&lhareader);
268 
269  // Remember position in vector of cross section estimates.
270  int iNow = sizeLO-1-njetcounterLO;
271 
272  // Start generation loop
273  for( int iEvent=0; iEvent<nEvent; ++iEvent ){
274 
275  // Generate next event
276  if( !pythia.next() ) {
277  if( pythia.info.atEndOfFile() ) break;
278  else continue;
279  }
280 
281  // Get event weight(s).
282  double weightNLO = pythia.info.mergingWeightNLO();
283  double evtweight = pythia.info.weight();
284  weightNLO *= evtweight;
285  // Do not print zero-weight events.
286  if ( weightNLO == 0. ) continue;
287 
288  // Construct new empty HepMC event.
289  HepMC::GenEvent* hepmcevt = new HepMC::GenEvent();
290  // Get correct cross section from previous estimate.
291  double normhepmc = xsecLO[iNow] / nAcceptLO[iNow];
292  // powheg box weighted events
293  if( abs(pythia.info.lhaStrategy()) == 4 )
294  normhepmc = 1. / (1e9*nSelectedLO[iNow]);
295  // Set event weight.
296  hepmcevt->weights().push_back(weightNLO*normhepmc);
297  // Fill HepMC event.
298  ToHepMC.fill_next_event( pythia, hepmcevt );
299  // Add the weight of the current event to the cross section.
300  sigmaTotal += weightNLO*normhepmc;
301  sigmaTemp += weightNLO*normhepmc;
302  errorTotal += pow2(weightNLO*normhepmc);
303  // Report cross section to hepmc
304  HepMC::GenCrossSection xsec;
305  xsec.set_cross_section( sigmaTotal*1e9, pythia.info.sigmaErr()*1e9 );
306  hepmcevt->set_cross_section( xsec );
307  // Write the HepMC event to file. Done with it.
308  ascii_io << hepmcevt;
309  delete hepmcevt;
310 
311  } // end loop over events to generate
312 
313  // print cross section, errors
314  pythia.stat();
315  // Save sample cross section for output.
316  sampleXStree.push_back(sigmaTemp);
317  sigmaTemp = 0.;
318 
319  // Restart with ME of a reduced the number of jets
320  if( njetcounterLO > 0 )
321  njetcounterLO--;
322  else
323  break;
324 
325  }
326 
327 //-----------------------------------------------------------------------------
328 //-----------------------------------------------------------------------------
329 
330  cout << endl << endl << endl;
331  cout << "Start nl3 virtual corrections part" << endl;
332 
333  // Switch on loop-level processing.
334  pythia.settings.flag("Merging:doNL3Tree",false);
335  pythia.settings.flag("Merging:doNL3Loop",true);
336  pythia.settings.flag("Merging:doNL3Subt",false);
337  pythia.settings.mode("Merging:nRecluster",0);
338 
339  njetcounterNLO = nMaxNLO;
340  iPathLoop = iPath + "_powheg";
341 
342  while(njetcounterNLO >= 0){
343 
344  // From njetcounter, choose LHE file
345  stringstream in;
346  in << "_" << njetcounterNLO << ".lhe";
347 #ifdef GZIPSUPPORT
348  if(access( (iPathLoop+in.str()+".gz").c_str(), F_OK) != -1) in << ".gz";
349 #endif
350  string LHEfile = iPathLoop + in.str();
351 
352  cout << endl << endl << endl
353  << "Start loop level treatment for " << njetcounterNLO << " jets"
354  << endl;
355 
356  LHAupLHEF lhareader((char*)(LHEfile).c_str());
357  pythia.settings.mode("Merging:nRequested", njetcounterNLO);
358  pythia.settings.word("Beams:LHEF", LHEfile);
359  pythia.init(&lhareader);
360 
361  // Remember position in vector of cross section estimates.
362  int iNow = sizeNLO-1-njetcounterNLO;
363 
364  // Start generation loop
365  for( int iEvent=0; iEvent<nEvent; ++iEvent ){
366 
367  // Generate next event
368  if( !pythia.next() ) {
369  if( pythia.info.atEndOfFile() ) break;
370  else continue;
371  }
372 
373  // Get event weight(s).
374  double weightNLO = pythia.info.mergingWeightNLO();
375  double evtweight = pythia.info.weight();
376  weightNLO *= evtweight;
377  // Do not print zero-weight events.
378  if ( weightNLO == 0. ) continue;
379 
380  // Construct new empty HepMC event.
381  HepMC::GenEvent* hepmcevt = new HepMC::GenEvent();
382  // Get correct cross section from previous estimate.
383  double normhepmc = xsecNLO[iNow] / nAcceptNLO[iNow];
384  // powheg box weighted events
385  if( abs(pythia.info.lhaStrategy()) == 4 )
386  normhepmc = 1. / (1e9*nSelectedNLO[iNow]);
387  // Set event weight.
388  hepmcevt->weights().push_back(weightNLO*normhepmc);
389  // Fill HepMC event.
390  ToHepMC.fill_next_event( pythia, hepmcevt );
391  // Add the weight of the current event to the cross section.
392  sigmaTotal += weightNLO*normhepmc;
393  sigmaTemp += weightNLO*normhepmc;
394  errorTotal += pow2(weightNLO*normhepmc);
395  // Report cross section to hepmc
396  HepMC::GenCrossSection xsec;
397  xsec.set_cross_section( sigmaTotal*1e9, pythia.info.sigmaErr()*1e9 );
398  hepmcevt->set_cross_section( xsec );
399  // Write the HepMC event to file. Done with it.
400  ascii_io << hepmcevt;
401  delete hepmcevt;
402 
403  } // end loop over events to generate
404 
405  // print cross section, errors
406  pythia.stat();
407  // Save sample cross section for output.
408  sampleXSvirt.push_back(sigmaTemp);
409  sigmaTemp = 0.;
410 
411  // Restart with ME of a reduced the number of jets
412  if( njetcounterNLO > 0)
413  njetcounterNLO--;
414  else
415  break;
416 
417  }
418 
419 //-----------------------------------------------------------------------------
420 //-----------------------------------------------------------------------------
421 
422  cout << endl << endl << endl;
423  cout << "Shower subtractive events" << endl;
424 
425  // Switch on processing of counter-events.
426  pythia.settings.flag("Merging:doNL3Tree",false);
427  pythia.settings.flag("Merging:doNL3Loop",false);
428  pythia.settings.flag("Merging:doNL3Subt",true);
429  pythia.settings.mode("Merging:nRecluster",1);
430  pythia.settings.mode("Merging:nRequested", -1);
431 
432  int nMaxCT = nMaxNLO + 1;
433  int njetcounterCT = nMaxCT;
434  string iPathSubt = iPath + "_tree";
435 
436  while(njetcounterCT >= 1){
437 
438  // From njetcounter, choose LHE file
439  stringstream in;
440  in << "_" << njetcounterCT << ".lhe";
441 #ifdef GZIPSUPPORT
442  if(access( (iPathSubt+in.str()+".gz").c_str(), F_OK) != -1) in << ".gz";
443 #endif
444  string LHEfile = iPathSubt + in.str();
445 
446  cout << endl << endl << endl
447  << "Start subtractive treatment for " << njetcounterCT << " jets"
448  << endl;
449 
450  LHAupLHEF lhareader((char*)(LHEfile).c_str());
451  pythia.settings.mode("Merging:nRequested", njetcounterCT);
452  pythia.settings.word("Beams:LHEF", LHEfile);
453  pythia.init(&lhareader);
454 
455  // Remember position in vector of cross section estimates.
456  int iNow = sizeLO-1-njetcounterCT;
457 
458  // Start generation loop
459  for( int iEvent=0; iEvent<nEvent; ++iEvent ){
460 
461  // Generate next event
462  if( !pythia.next() ) {
463  if( pythia.info.atEndOfFile() ) break;
464  else continue;
465  }
466 
467  // Get event weight(s).
468  double weightNLO = pythia.info.mergingWeightNLO();
469  double evtweight = pythia.info.weight();
470  weightNLO *= evtweight;
471  // Do not print zero-weight events.
472  if ( weightNLO == 0. ) continue;
473 
474  // Construct new empty HepMC event.
475  HepMC::GenEvent* hepmcevt = new HepMC::GenEvent();
476  // Get correct cross section from previous estimate.
477  double normhepmc = -1.*xsecLO[iNow] / nAcceptLO[iNow];
478  // powheg box weighted events
479  if( abs(pythia.info.lhaStrategy()) == 4 )
480  normhepmc = -1. / (1e9*nSelectedLO[iNow]);
481  // Set event weight.
482  hepmcevt->weights().push_back( weightNLO*normhepmc);
483  // Fill HepMC event.
484  ToHepMC.fill_next_event( pythia, hepmcevt );
485  // Add the weight of the current event to the cross section.
486  sigmaTotal += weightNLO*normhepmc;
487  sigmaTemp += weightNLO*normhepmc;
488  errorTotal += pow2(weightNLO*normhepmc);
489  // Report cross section to hepmc
490  HepMC::GenCrossSection xsec;
491  xsec.set_cross_section( sigmaTotal*1e9, pythia.info.sigmaErr()*1e9 );
492  hepmcevt->set_cross_section( xsec );
493 
494  // Write the HepMC event to file. Done with it.
495  ascii_io << hepmcevt;
496  delete hepmcevt;
497 
498  } // end loop over events to generate
499 
500  // print cross section, errors
501  pythia.stat();
502  // Save sample cross section for output.
503  sampleXSsubtTree.push_back(sigmaTemp);
504  sigmaTemp = 0.;
505 
506  // Restart with ME of a reduced the number of jets
507  if( njetcounterCT > 1 )
508  njetcounterCT--;
509  else
510  break;
511 
512  }
513 
514  // Print cross section information.
515  cout << endl << endl;
516  cout << " *---------------------------------------------------*" << endl;
517  cout << " | |" << endl;
518  cout << " | Sample cross sections after NL3 merging |" << endl;
519  cout << " | |" << endl;
520  cout << " | Leading order cross sections (mb): |" << endl;
521  for (int i = 0; i < int(sampleXStree.size()); ++i)
522  cout << " | " << sampleXStree.size()-1-i << "-jet: "
523  << setw(17) << scientific << setprecision(6)
524  << sampleXStree[i] << " |" << endl;
525  cout << " | |" << endl;
526  cout << " | NLO order cross sections (mb): |" << endl;
527  for (int i = 0; i < int(sampleXSvirt.size()); ++i)
528  cout << " | " << sampleXSvirt.size()-1-i << "-jet: "
529  << setw(17) << scientific << setprecision(6)
530  << sampleXSvirt[i] << " |" << endl;
531  cout << " | |" << endl;
532  cout << " | Leading-order subtractive cross sections (mb): |" << endl;
533  for (int i = 0; i < int(sampleXSsubtTree.size()); ++i)
534  cout << " | " << sampleXSsubtTree.size()-1-i+1 << "-jet: "
535  << setw(17) << scientific << setprecision(6)
536  << sampleXSsubtTree[i] << " |" << endl;
537  cout << " | |" << endl;
538  cout << " |---------------------------------------------------|" << endl;
539  cout << " |---------------------------------------------------|" << endl;
540  cout << " | Inclusive cross sections: |" << endl;
541  cout << " | |" << endl;
542  cout << " | NL3 merged inclusive cross section: |" << endl;
543  cout << " | " << setw(17) << scientific << setprecision(6)
544  << sigmaTotal << " +- " << setw(17) << sqrt(errorTotal) << " mb "
545  << " |" << endl;
546  cout << " | |" << endl;
547  cout << " | NLO inclusive cross section: |" << endl;
548  cout << " | " << setw(17) << scientific << setprecision(6)
549  << xsecNLO.back() << " mb |" << endl;
550  cout << " | |" << endl;
551  cout << " | LO inclusive cross section: |" << endl;
552  cout << " | " << setw(17) << scientific << setprecision(6)
553  << xsecLO.back() << " mb |" << endl;
554  cout << " | |" << endl;
555  cout << " *---------------------------------------------------*" << endl;
556  cout << endl << endl;
557 
558  // Done
559  return 0;
560 
561 }
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