10 #include "Pythia8/Pythia.h"
11 #include "Pythia8/Pythia8ToHepMC.h"
14 #include "HepMC/GenEvent.h"
15 #include "HepMC/IO_GenEvent.h"
17 #include "HepMC/Units.h"
19 using namespace Pythia8;
25 int main(
int argc,
char* argv[] ){
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. Path for output histogram files" << endl
35 <<
" Program stopped. " << endl;
45 pythia.readFile(argv[1]);
52 ToHepMC.set_print_inconsistency(
false);
53 ToHepMC.set_free_parton_warnings(
false);
55 ToHepMC.set_store_pdf(
false);
56 ToHepMC.set_store_proc(
false);
57 ToHepMC.set_store_xsec(
false);
61 string iPath = string(argv[2]);
64 int nEvent = pythia.mode(
"Main:numberOfEvents");
66 int nMaxLO = pythia.mode(
"Merging:nJetMax");
68 int nMaxNLO = pythia.mode(
"Merging:nJetMaxNLO");
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);
85 pythia.settings.flag(
"Merging:doXSectionEstimate",
true);
86 pythia.settings.flag(
"Merging:doUNLOPSTree",
true);
88 int njetcounterLO = nMaxLO;
89 string iPathTree = iPath +
"_tree";
92 vector<double> xsecLO;
93 vector<double> nSelectedLO;
94 vector<double> nAcceptLO;
95 vector<int> strategyLO;
97 cout << endl << endl << endl;
98 cout <<
"Start estimating unlops tree level cross section" << endl;
100 while(njetcounterLO >= 0){
104 in <<
"_" << njetcounterLO <<
".lhe";
106 if(access( (iPathTree+in.str()+
".gz").c_str(), F_OK) != -1) in <<
".gz";
108 string LHEfile = iPathTree + in.str();
109 LHAupLHEF lhareader((
char*)(LHEfile).c_str());
110 pythia.settings.mode(
"Merging:nRequested", njetcounterLO);
111 pythia.settings.word(
"Beams:LHEF", LHEfile);
112 pythia.init(&lhareader);
115 for(
int iEvent=0; iEvent<nEvent; ++iEvent ){
117 if( !pythia.next() ) {
118 if( pythia.info.atEndOfFile() ){
128 xsecLO.push_back(pythia.info.sigmaGen());
129 nSelectedLO.push_back(pythia.info.nSelected());
130 nAcceptLO.push_back(pythia.info.nAccepted());
131 strategyLO.push_back(pythia.info.lhaStrategy());
134 if( njetcounterLO > 0 )
144 cout << endl << endl << endl;
145 cout <<
"Start estimating unlops virtual corrections cross section" << endl;
147 pythia.settings.flag(
"Merging:doUNLOPSTree",
false);
148 pythia.settings.flag(
"Merging:doUNLOPSLoop",
true);
150 int njetcounterNLO = nMaxNLO;
151 string iPathLoop = iPath +
"_powheg";
154 vector<double> xsecNLO;
155 vector<double> nSelectedNLO;
156 vector<double> nAcceptNLO;
157 vector<int> strategyNLO;
159 while(njetcounterNLO >= 0){
163 in <<
"_" << njetcounterNLO <<
".lhe";
165 if(access( (iPathLoop+in.str()+
".gz").c_str(), F_OK) != -1) in <<
".gz";
167 string LHEfile = iPathLoop + in.str();
168 LHAupLHEF lhareader((
char*)(LHEfile).c_str());
169 pythia.settings.mode(
"Merging:nRequested", njetcounterNLO);
170 pythia.settings.word(
"Beams:LHEF", LHEfile);
171 pythia.init(&lhareader);
174 for(
int iEvent=0; iEvent<nEvent; ++iEvent ){
176 if( !pythia.next() ) {
177 if( pythia.info.atEndOfFile() ){
187 xsecNLO.push_back(pythia.info.sigmaGen());
188 nSelectedNLO.push_back(pythia.info.nSelected());
189 nAcceptNLO.push_back(pythia.info.nAccepted());
190 strategyNLO.push_back(pythia.info.lhaStrategy());
193 if( njetcounterNLO > 0 )
200 int sizeLO = int(xsecLO.size());
201 int sizeNLO = int(xsecNLO.size());
204 pythia.settings.flag(
"Merging:doXSectionEstimate",
false);
207 pythia.settings.flag(
"PartonLevel:FSR",fsr);
208 pythia.settings.flag(
"PartonLevel:ISR",isr);
209 pythia.settings.flag(
"HadronLevel:all",had);
210 pythia.settings.flag(
"PartonLevel:MPI",mpi);
216 double sigmaTemp = 0.;
217 vector<double> sampleXStree;
218 vector<double> sampleXSvirt;
219 vector<double> sampleXSsubtTree;
220 vector<double> sampleXSsubtVirt;
222 double sigmaTotal = 0.;
223 double errorTotal = 0.;
226 pythia.settings.flag(
"Merging:doUNLOPSTree",
true);
227 pythia.settings.flag(
"Merging:doUNLOPSLoop",
false);
228 pythia.settings.flag(
"Merging:doUNLOPSSubt",
false);
229 pythia.settings.flag(
"Merging:doUNLOPSSubtNLO",
false);
230 pythia.settings.mode(
"Merging:nRecluster",0);
233 njetcounterLO = nMaxLO;
234 iPathTree = iPath +
"_tree";
236 while(njetcounterLO >= 0){
240 in <<
"_" << njetcounterLO <<
".lhe";
242 if(access( (iPathTree+in.str()+
".gz").c_str(), F_OK) != -1) in <<
".gz";
244 string LHEfile = iPathTree + in.str();
245 LHAupLHEF lhareader((
char*)(LHEfile).c_str());
247 cout << endl << endl << endl
248 <<
"Start tree level treatment for " << njetcounterLO <<
" jets"
252 if ( njetcounterLO == 0 )
break;
253 pythia.settings.mode(
"Merging:nRequested", njetcounterLO);
254 pythia.settings.word(
"Beams:LHEF", LHEfile);
255 pythia.init(&lhareader);
258 int iNow = sizeLO-1-njetcounterLO;
261 for(
int iEvent=0; iEvent<nEvent; ++iEvent ){
264 if( !pythia.next() ) {
265 if( pythia.info.atEndOfFile() )
break;
270 double weightNLO = pythia.info.mergingWeightNLO();
271 double evtweight = pythia.info.weight();
272 weightNLO *= evtweight;
274 if ( weightNLO == 0. )
continue;
279 double normhepmc = xsecLO[iNow] / nAcceptLO[iNow];
281 hepmcevt->
weights().push_back(weightNLO*normhepmc);
283 ToHepMC.fill_next_event( pythia, hepmcevt );
285 sigmaTotal += weightNLO*normhepmc;
286 sigmaTemp += weightNLO*normhepmc;
287 errorTotal += pow2(weightNLO*normhepmc);
293 ascii_io << hepmcevt;
301 if( njetcounterLO > 0 )
307 sampleXStree.push_back(sigmaTemp);
315 cout << endl << endl << endl;
316 cout <<
"Start unlops virtual corrections part" << endl;
319 pythia.settings.flag(
"Merging:doUNLOPSTree",
false);
320 pythia.settings.flag(
"Merging:doUNLOPSLoop",
true);
321 pythia.settings.flag(
"Merging:doUNLOPSSubt",
false);
322 pythia.settings.flag(
"Merging:doUNLOPSSubtNLO",
false);
323 pythia.settings.mode(
"Merging:nRecluster",0);
325 njetcounterNLO = nMaxNLO;
326 iPathLoop= iPath +
"_powheg";
328 while(njetcounterNLO >= 0){
332 in <<
"_" << njetcounterNLO <<
".lhe";
334 if(access( (iPathLoop+in.str()+
".gz").c_str(), F_OK) != -1) in <<
".gz";
336 string LHEfile = iPathLoop + in.str();
337 LHAupLHEF lhareader((
char*)(LHEfile).c_str());
339 cout << endl << endl << endl
340 <<
"Start loop level treatment for " << njetcounterNLO <<
" jets"
343 pythia.settings.mode(
"Merging:nRequested", njetcounterNLO);
344 pythia.settings.word(
"Beams:LHEF", LHEfile);
345 pythia.init(&lhareader);
348 int iNow = sizeNLO-1-njetcounterNLO;
351 for(
int iEvent=0; iEvent<nEvent; ++iEvent ){
354 if( !pythia.next() ) {
355 if( pythia.info.atEndOfFile() )
break;
360 double weightNLO = pythia.info.mergingWeightNLO();
361 double evtweight = pythia.info.weight();
362 weightNLO *= evtweight;
364 if ( weightNLO == 0. )
continue;
369 double normhepmc = xsecNLO[iNow] / nAcceptNLO[iNow];
371 if( abs(strategyNLO[iNow]) == 4)
372 normhepmc = 1. / (1e9*nSelectedNLO[iNow]);
374 hepmcevt->
weights().push_back(weightNLO*normhepmc);
376 ToHepMC.fill_next_event( pythia, hepmcevt );
378 sigmaTotal += weightNLO*normhepmc;
379 sigmaTemp += weightNLO*normhepmc;
380 errorTotal += pow2(weightNLO*normhepmc);
386 ascii_io << hepmcevt;
394 sampleXSvirt.push_back(sigmaTemp);
398 if( njetcounterNLO > 0)
408 cout << endl << endl << endl;
409 cout <<
"Shower subtractive events" << endl;
412 pythia.settings.flag(
"Merging:doUNLOPSTree",
false);
413 pythia.settings.flag(
"Merging:doUNLOPSLoop",
false);
414 pythia.settings.flag(
"Merging:doUNLOPSSubt",
true);
415 pythia.settings.flag(
"Merging:doUNLOPSSubtNLO",
false);
416 pythia.settings.mode(
"Merging:nRecluster",1);
419 int njetcounterCT = nMaxCT;
420 string iPathSubt= iPath +
"_tree";
422 while(njetcounterCT >= 1){
426 in <<
"_" << njetcounterCT <<
".lhe";
428 if(access( (iPathSubt+in.str()+
".gz").c_str(), F_OK) != -1) in <<
".gz";
430 string LHEfile = iPathSubt + in.str();
431 LHAupLHEF lhareader((
char*)(LHEfile).c_str());
433 cout << endl << endl << endl
434 <<
"Start subtractive treatment for " << njetcounterCT <<
" jets"
437 pythia.settings.mode(
"Merging:nRequested", njetcounterCT);
438 pythia.settings.word(
"Beams:LHEF", LHEfile);
439 pythia.init(&lhareader);
442 int iNow = sizeLO-1-njetcounterCT;
445 for(
int iEvent=0; iEvent<nEvent; ++iEvent ){
448 if( !pythia.next() ) {
449 if( pythia.info.atEndOfFile() )
break;
454 double weightNLO = pythia.info.mergingWeightNLO();
455 double evtweight = pythia.info.weight();
456 weightNLO *= evtweight;
458 if ( weightNLO == 0. )
continue;
463 double normhepmc = -1*xsecLO[iNow] / nAcceptLO[iNow];
465 hepmcevt->
weights().push_back(weightNLO*normhepmc);
467 ToHepMC.fill_next_event( pythia, hepmcevt );
469 sigmaTotal += weightNLO*normhepmc;
470 sigmaTemp += weightNLO*normhepmc;
471 errorTotal += pow2(weightNLO*normhepmc);
477 ascii_io << hepmcevt;
485 sampleXSsubtTree.push_back(sigmaTemp);
489 if( njetcounterCT > 1 )
499 cout << endl << endl << endl;
500 cout <<
"Shower subtractive events" << endl;
502 pythia.settings.flag(
"Merging:doUNLOPSTree",
false);
503 pythia.settings.flag(
"Merging:doUNLOPSLoop",
false);
504 pythia.settings.flag(
"Merging:doUNLOPSSubt",
false);
505 pythia.settings.flag(
"Merging:doUNLOPSSubtNLO",
true);
506 pythia.settings.mode(
"Merging:nRecluster",1);
509 njetcounterCT = nMaxCT;
510 iPathSubt= iPath +
"_powheg";
512 while(njetcounterCT >= 1){
516 in <<
"_" << njetcounterCT <<
".lhe";
518 if(access( (iPathSubt+in.str()+
".gz").c_str(), F_OK) != -1) in <<
".gz";
520 string LHEfile = iPathSubt + in.str();
521 LHAupLHEF lhareader((
char*)(LHEfile).c_str());
523 cout << endl << endl << endl
524 <<
"Start subtractive treatment for " << njetcounterCT <<
" nlo jets"
527 pythia.settings.mode(
"Merging:nRequested", njetcounterCT);
528 pythia.settings.word(
"Beams:LHEF", LHEfile);
529 pythia.init(&lhareader);
532 int iNow = sizeNLO-1-njetcounterCT;
535 for(
int iEvent=0; iEvent<nEvent; ++iEvent ){
538 if( !pythia.next() ) {
539 if( pythia.info.atEndOfFile() )
break;
544 double weightNLO = pythia.info.mergingWeightNLO();
545 double evtweight = pythia.info.weight();
546 weightNLO *= evtweight;
548 if ( weightNLO == 0. )
continue;
553 double normhepmc = -1*xsecNLO[iNow] / nAcceptNLO[iNow];
555 if( abs(strategyNLO[iNow]) == 4)
556 normhepmc = -1. / (1e9*nSelectedNLO[iNow]);
558 hepmcevt->
weights().push_back(weightNLO*normhepmc);
560 ToHepMC.fill_next_event( pythia, hepmcevt );
562 sigmaTotal += weightNLO*normhepmc;
563 sigmaTemp += weightNLO*normhepmc;
564 errorTotal += pow2(weightNLO*normhepmc);
570 ascii_io << hepmcevt;
578 sampleXSsubtVirt.push_back(sigmaTemp);
582 if( njetcounterCT > 1 )
590 cout << endl << endl;
591 cout <<
" *---------------------------------------------------*" << endl;
592 cout <<
" | |" << endl;
593 cout <<
" | Sample cross sections after UNLOPS merging |" << endl;
594 cout <<
" | |" << endl;
595 cout <<
" | Leading order cross sections (mb): |" << endl;
596 for (
int i = 0; i < int(sampleXStree.size()); ++i)
597 cout <<
" | " << sampleXStree.size()-1-i+1 <<
"-jet: "
598 << setw(17) << scientific << setprecision(6)
599 << sampleXStree[i] <<
" |" << endl;
600 cout <<
" | (No 0-jet tree-level sample in UNLOPS) |" << endl;
601 cout <<
" | |" << endl;
602 cout <<
" | NLO order cross sections (mb): |" << endl;
603 for (
int i = 0; i < int(sampleXSvirt.size()); ++i)
604 cout <<
" | " << sampleXSvirt.size()-1-i <<
"-jet: "
605 << setw(17) << scientific << setprecision(6)
606 << sampleXSvirt[i] <<
" |" << endl;
607 cout <<
" | |" << endl;
608 cout <<
" | Leading-order subtractive cross sections (mb): |" << endl;
609 for (
int i = 0; i < int(sampleXSsubtTree.size()); ++i)
610 cout <<
" | " << sampleXSsubtTree.size()-1-i+1 <<
"-jet: "
611 << setw(17) << scientific << setprecision(6)
612 << sampleXSsubtTree[i] <<
" |" << endl;
613 cout <<
" | |" << endl;
614 if ( sampleXSsubtVirt.size() > 0) {
615 cout <<
" | NLO subtractive cross sections (mb): |" << endl;
616 for (
int i = 0; i < int(sampleXSsubtVirt.size()); ++i)
617 cout <<
" | " << sampleXSsubtVirt.size()-1-i+1 <<
"-jet: "
618 << setw(17) << scientific << setprecision(6)
619 << sampleXSsubtVirt[i] <<
" |" << endl;
620 cout <<
" | |" << endl;
622 cout <<
" |---------------------------------------------------|" << endl;
623 cout <<
" |---------------------------------------------------|" << endl;
624 cout <<
" | Inclusive cross sections: |" << endl;
625 cout <<
" | |" << endl;
626 cout <<
" | UNLOPS merged inclusive cross section: |" << endl;
627 cout <<
" | " << setw(17) << scientific << setprecision(6)
628 << sigmaTotal <<
" +- " << setw(17) << sqrt(errorTotal) <<
" mb "
630 cout <<
" | |" << endl;
631 cout <<
" | NLO inclusive cross section: |" << endl;
632 cout <<
" | " << setw(17) << scientific << setprecision(6)
633 << xsecNLO.back() <<
" mb |" << endl;
634 cout <<
" | |" << endl;
635 cout <<
" | LO inclusive cross section: |" << endl;
636 cout <<
" | " << setw(17) << scientific << setprecision(6)
637 << xsecLO.back() <<
" mb |" << endl;
638 cout <<
" | |" << endl;
639 cout <<
" *---------------------------------------------------*" << endl;
640 cout << endl << endl;
void set_cross_section(double xs, double xs_err)
Set cross section and error in pb.
The GenCrossSection class stores the generated cross section.
void set_cross_section(const GenCrossSection &)
provide a pointer to the GenCrossSection container
The GenEvent class is the core of HepMC.
IO_GenEvent also deals with HeavyIon and PdfInfo.
WeightContainer & weights()
direct access to WeightContainer