StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
main83.cc
1 // main83.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2014 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 "Pythia8/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
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 }