source: HiSusy/trunk/Pythia8/pythia8170/examples/main82.cc @ 1

Last change on this file since 1 was 1, checked in by zerwas, 11 years ago

first import of structure, PYTHIA8 and DELPHES

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