PYTHIA8/pythia8175/examples/main81.cc

   1 // main81.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2013 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 "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 }