PYTHIA8/pythia8170/examples/main83.cc

   1 // main83.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2012 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 te 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 // Class for user interaction with the merging
  82
  83 class MyMergingHooks : public MergingHooks {
  84
  85 private:
  86
  87 public:
  88
  89   // Default constructor
  90   MyMergingHooks();
  91   // Destructor
  92   ~MyMergingHooks();
  93
  94   // Functional definition of the merging scale
  95   virtual double tmsDefinition( const Event& event);
  96
  97   // Function to dampen weights calculated from histories with lowest
  98   // multiplicity reclustered events that do not pass the ME cuts
  99   virtual double dampenIfFailCuts( const Event& inEvent );
 100
 101   // Helper function for tms definition
 102   double myKTdurham(const Particle& RadAfterBranch,
 103            const Particle& EmtAfterBranch, int Type, double D );
 104
 105 };
 106
 107 //--------------------------------------------------------------------------
 108
 109 // Constructor
 110 MyMergingHooks::MyMergingHooks() {}
 111
 112 // Desctructor
 113 MyMergingHooks::~MyMergingHooks() {}
 114
 115 //--------------------------------------------------------------------------
 116
 117 double MyMergingHooks::dampenIfFailCuts( const Event& inEvent ){
 118
 119   // Get pT for pure QCD 2->2 state
 120   double pT = 0.;
 121   for( int i=0; i < inEvent.size(); ++i)
 122     if(inEvent[i].isFinal() && inEvent[i].colType() != 0) {
 123       pT = sqrt(pow(inEvent[i].px(),2) + pow(inEvent[i].py(),2));
 124       break;
 125     }
 126
 127   // Veto history if lowest multiplicity event does not pass ME cuts
 128   if(pT < 10.) return 0.;
 129
 130   return 1.;
 131
 132 }
 133
 134 //--------------------------------------------------------------------------
 135
 136 // Definition of the merging scale
 137
 138 double MyMergingHooks::tmsDefinition( const Event& event){
 139
 140   // Cut only on QCD partons!
 141   // Count particle types
 142   int nFinalColoured = 0;
 143   int nFinalNow =0;
 144   for( int i=0; i < event.size(); ++i) {
 145     if(event[i].isFinal()){
 146       if(event[i].id() != 23 && abs(event[i].id()) != 24)
 147         nFinalNow++;
 148       if( event[i].colType() != 0)
 149         nFinalColoured++;
 150     }
 151   }
 152
 153   // Use MergingHooks in-built functions to get information on the hard process
 154   int nLeptons = nHardOutLeptons();
 155   int nQuarks  = nHardOutPartons();
 156   int nResNow  = nResInCurrent();
 157
 158   // Check if photons, electrons etc. have been produced. If so, do not veto
 159   if(nFinalNow - ( (nLeptons+nQuarks)/2 - nResNow)*2 != nFinalColoured){
 160     // Sometimes, Pythia detaches the decay products even though no
 161     // resonance was put into the LHE file, to catch this, add another
 162     // if statement
 163     if(nFinalNow != nFinalColoured) return 0.;
 164   }
 165
 166   // Check that one parton has been produced. If not (e.g. in MPI), do not veto
 167   int nMPI = infoPtr->nMPI();
 168   if(nMPI > 1) return 0.;
 169
 170   // Declare kT algorithm parameters
 171   double Dparam = 0.4;
 172   int kTtype = -1;
 173   // Declare final parton vector
 174   vector <int> FinalPartPos;
 175   FinalPartPos.clear();
 176   // Search event record for final state partons
 177   for (int i=0; i < event.size(); ++i)
 178     if(event[i].isFinal() && event[i].colType() != 0)
 179       FinalPartPos.push_back(i);
 180
 181   // Find minimal Durham kT in event, using own function: Check
 182   // definition of separation
 183   int type = (event[3].colType() == 0 && event[4].colType() == 0) ? 1 : kTtype;
 184   // Find minimal kT
 185   double ktmin = event[0].e();
 186   for(int i=0; i < int(FinalPartPos.size()); ++i){
 187     double kt12  = ktmin;
 188     // Compute separation to the beam axis for hadronic collisions
 189     if(type == -1 || type == -2) {
 190       double temp = event[FinalPartPos[i]].pT();
 191       kt12 = min(kt12, temp);
 192     }
 193     // Compute separation to other final state jets
 194     for(int j=i+1; j < int(FinalPartPos.size()); ++j) {
 195       double temp = kTdurham( event[FinalPartPos[i]], event[FinalPartPos[j]],
 196                       type, Dparam);
 197       kt12 = min(kt12, temp);
 198     }
 199     // Keep the minimal Durham separation
 200     ktmin = min(ktmin,kt12);
 201   }
 202
 203   // Return minimal Durham kT
 204   return ktmin;
 205
 206 }
 207
 208 //--------------------------------------------------------------------------
 209
 210 // Function to compute durham y separation from Particle input
 211
 212 double MyMergingHooks::myKTdurham(const Particle& RadAfterBranch,
 213   const Particle& EmtAfterBranch, int Type, double D ){
 214
 215   // Declare return variable
 216   double ktdur;
 217   // Save 4-momenta of final state particles
 218   Vec4 jet1 = RadAfterBranch.p();
 219   Vec4 jet2 = EmtAfterBranch.p();
 220
 221   if( Type == 1) {
 222     // Get angle between jets for e+e- collisions, make sure that
 223     // -1 <= cos(theta) <= 1
 224     double costh;
 225     if (jet1.pAbs()*jet2.pAbs() <=0.) costh = 1.;
 226     else {
 227       costh = costheta(jet1,jet2);
 228     }
 229     // Calculate kt durham separation between jets for e+e- collisions
 230     ktdur = 2.0*min( pow(jet1.e(),2) , (pow(jet2.e(),2)) )*(1.0 - costh);
 231   } else if( Type == -1 ){
 232     // Get delta_eta and cosh(Delta_eta) for hadronic collisions
 233     double eta1 = 0.5*log( (jet1.e() + jet1.pz()) / (jet1.e() - jet1.pz()) );
 234     double eta2 = 0.5*log( (jet2.e() + jet2.pz()) / (jet2.e() - jet2.pz()) );
 235     // Get delta_phi and cos(Delta_phi) for hadronic collisions
 236     double pt1 = sqrt( pow(jet1.px(),2) + pow(jet1.py(),2) );
 237     double pt2 = sqrt( pow(jet2.px(),2) + pow(jet2.py(),2) );
 238     double cosdPhi = ( jet1.px()*jet2.px() + jet1.py()*jet2.py() ) / (pt1*pt2);
 239     double dPhi = acos( cosdPhi );
 240     // Calculate kT durham like fastjet
 241      ktdur = min( pow(pt1,2),pow(pt2,2) )
 242            * ( pow(eta1-eta2,2) + pow(dPhi,2) ) / pow(D,2);
 243   } else if( Type == -2 ){
 244     // Get delta_eta and cosh(Delta_eta) for hadronic collisions
 245     double eta1 = 0.5*log( (jet1.e() + jet1.pz()) / (jet1.e() - jet1.pz()) );
 246     double eta2 = 0.5*log( (jet2.e() + jet2.pz()) / (jet2.e() - jet2.pz()) );
 247      double coshdEta = cosh( eta1 - eta2 );
 248     // Get delta_phi and cos(Delta_phi) for hadronic collisions
 249     double pt1 = sqrt( pow(jet1.px(),2) + pow(jet1.py(),2) );
 250     double pt2 = sqrt( pow(jet2.px(),2) + pow(jet2.py(),2) );
 251     double cosdPhi = ( jet1.px()*jet2.px() + jet1.py()*jet2.py() ) / (pt1*pt2);
 252     // Calculate kT durham separation "SHERPA-like"
 253      ktdur = 2.0*min( pow(pt1,2),pow(pt2,2) )
 254            * ( coshdEta - cosdPhi ) / pow(D,2);
 255   } else {
 256     ktdur = 0.0;
 257   }
 258   // Return kT
 259   return sqrt(ktdur);
 260 }
 261
 262 //==========================================================================
 263
 264 // Example main programm to illustrate merging
 265
 266 int main( int argc, char* argv[] ){
 267
 268   // Check that correct number of command-line arguments
 269   if (argc != 4) {
 270     cerr << " Unexpected number of command-line arguments. \n You are"
 271          << " expected to provide the arguments \n"
 272          << " 1. Input file for settings \n"
 273          << " 2. Full name of the input LHE file (with path) \n"
 274          << " 3. Path for output histogram files \n"
 275          << " Program stopped. " << endl;
 276     return 1;
 277   }
 278
 279   Pythia pythia;
 280
 281   // Input parameters:
 282   //  1. Input file for settings
 283   //  2. Path to input LHE file
 284   //  3. OUtput histogram path
 285   pythia.readFile(argv[1]);
 286   string iPath = string(argv[2]);
 287   string oPath = string(argv[3]);
 288
 289   // Number of events
 290   int nEvent = pythia.mode("Main:numberOfEvents");
 291
 292   // Construct user inut for merging
 293   MergingHooks* myMergingHooks = new MyMergingHooks();
 294   pythia.setMergingHooksPtr( myMergingHooks );
 295
 296   // For ISR regularisation off
 297   pythia.settings.forceParm("SpaceShower:pT0Ref",0.);
 298
 299   // Declare histograms
 300   Hist histPTFirst("pT of first jet",100,0.,100.);
 301   Hist histPTSecond("pT of second jet",100,0.,100.);
 302
 303   // Read in ME configurations
 304   pythia.init(iPath,false);
 305
 306   if(pythia.flag("Main:showChangedSettings")) {
 307     pythia.settings.listChanged();
 308   }
 309
 310   // Start generation loop
 311   for( int iEvent=0; iEvent<nEvent; ++iEvent ){
 312
 313     // Generate next event
 314     if( ! pythia.next()) continue;
 315
 316     // Get CKKWL weight of current event
 317     double weight = pythia.info.mergingWeight();
 318
 319     // Fill bins with CKKWL weight
 320     double pTfirst = pTfirstJet(pythia.event,1, 0.4);
 321     double pTsecnd = pTfirstJet(pythia.event,2, 0.4);
 322     histPTFirst.fill( pTfirst, weight);
 323     histPTSecond.fill( pTsecnd, weight);
 324
 325     if(iEvent%1000 == 0) cout << iEvent << endl;
 326
 327   } // end loop over events to generate
 328
 329   // print cross section, errors
 330   pythia.stat();
 331
 332   // Normalise histograms
 333   double norm = 1.
 334               * pythia.info.sigmaGen()
 335               * 1./ double(nEvent);
 336
 337   histPTFirst           *= norm;
 338   histPTSecond          *= norm;
 339
 340   // Get the number of jets in the LHE file from the file name
 341   string jetsInLHEF = iPath.substr(iPath.size()-5, iPath.size());
 342   jetsInLHEF = jetsInLHEF.substr(0, jetsInLHEF.size()-4);
 343
 344   // Write histograms to dat file. Use "jetsInLHEF" to label the files
 345   // Once all the samples up to the maximal desired jet multiplicity from the
 346   // matrix element are run, add all histograms to produce a
 347   // matrix-element-merged prediction
 348
 349   ofstream write;
 350   stringstream suffix;
 351   suffix << jetsInLHEF << "_wv.dat";
 352
 353   // Write histograms to file
 354   write.open( (char*)(oPath + "PTjet1_" + suffix.str()).c_str());
 355   histPTFirst.table(write);
 356   write.close();
 357
 358   write.open( (char*)(oPath + "PTjet2_" + suffix.str()).c_str());
 359   histPTSecond.table(write);
 360   write.close();
 361
 362
 363   delete myMergingHooks;
 364   return 0;
 365
 366   // Done
 367 }