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