Update to pythi8.170
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8170 / examples / main84.cc
1 // main84.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 <time.h>
11 #include "Pythia.h"
12
13 #include "HepMCInterface.h"
14 #include "HepMC/GenEvent.h"
15 #include "HepMC/IO_GenEvent.h"
16 // Following line to be used with HepMC 2.04 onwards.
17 #include "HepMC/Units.h"
18
19 using namespace Pythia8;
20
21 // Functions for histogramming
22 #include "fastjet/PseudoJet.hh"
23 #include "fastjet/ClusterSequence.hh"
24 #include "fastjet/CDFMidPointPlugin.hh"
25 #include "fastjet/CDFJetCluPlugin.hh"
26 #include "fastjet/D0RunIIConePlugin.hh"
27
28 //==========================================================================
29
30 // Find the Durham kT separation of the clustering from
31 // nJetMin --> nJetMin-1 jets in te input event  
32
33 double pTfirstJet( const Event& event, int nJetMin, double Rparam) {
34
35   double yPartonMax = 4.;
36
37   // Fastjet analysis - select algorithm and parameters
38   fastjet::Strategy               strategy = fastjet::Best;
39   fastjet::RecombinationScheme    recombScheme = fastjet::E_scheme;
40   fastjet::JetDefinition         *jetDef = NULL;
41   // For hadronic collision, use hadronic Durham kT measure
42   if(event[3].colType() != 0 || event[4].colType() != 0)
43     jetDef = new fastjet::JetDefinition(fastjet::kt_algorithm, Rparam,
44                                       recombScheme, strategy);
45   // For e+e- collision, use e+e- Durham kT measure
46   else
47     jetDef = new fastjet::JetDefinition(fastjet::ee_kt_algorithm,
48                                       recombScheme, strategy);
49   // Fastjet input
50   std::vector <fastjet::PseudoJet> fjInputs;
51   // Reset Fastjet input
52   fjInputs.resize(0);
53
54   // Loop over event record to decide what to pass to FastJet
55   for (int i = 0; i < event.size(); ++i) {
56     // (Final state && coloured+photons) only!
57     if ( !event[i].isFinal()
58       || event[i].isLepton()
59       || event[i].id() == 23
60       || abs(event[i].id()) == 24
61       || abs(event[i].y()) > yPartonMax)
62       continue;
63
64     // Store as input to Fastjet
65     fjInputs.push_back( fastjet::PseudoJet (event[i].px(),
66             event[i].py(), event[i].pz(),event[i].e() ) );
67   }
68
69   // Do nothing for empty input 
70   if (int(fjInputs.size()) == 0) {
71     delete jetDef;
72     return 0.0;
73   }
74
75   // Run Fastjet algorithm
76   fastjet::ClusterSequence clustSeq(fjInputs, *jetDef);
77   // Extract kT of first clustering
78   double pTFirst = sqrt(clustSeq.exclusive_dmerge_max(nJetMin-1));
79
80   delete jetDef;
81   // Return kT
82   return pTFirst;
83
84 }
85
86 //==========================================================================
87
88 int main( int argc, char* argv[] ){
89
90   // Check that correct number of command-line arguments
91   if (argc != 6) {
92     cerr << " Unexpected number of command-line arguments. \n You are"
93          << " expected to provide the arguments" << endl
94          << " 1. Input file for settings" << endl
95          << " 2. Name of output HepMC file" << endl
96          << " 3. Maximal number of additional jets"
97          << " (not used internally in Pythia, only used to construct the full"
98          << " name of lhe files with additional jets, and to label output"
99          << " histograms)" << endl
100          << " 4. Full name of the input LHE file (with path"
101          << " , without any _0.lhe suffix)" << endl
102          << " 5. Path for output histogram files" << endl
103          << " Program stopped. " << endl;
104     return 1;
105   }
106
107   Pythia pythia;
108
109   // First argument: Get input from an input file
110   pythia.readFile(argv[1]);
111   int nEvent = pythia.mode("Main:numberOfEvents");
112
113   // Interface for conversion from Pythia8::Event to HepMC one. 
114   HepMC::I_Pythia8 ToHepMC;
115   //  ToHepMC.set_crash_on_problem();
116   // Specify file where HepMC events will be stored.
117   HepMC::IO_GenEvent ascii_io(argv[2], std::ios::out);
118   // Following two lines are deprecated alternative.
119   // HepMC::IO_Ascii ascii_io(argv[2], std::ios::out);
120   // HepMC::IO_AsciiParticles ascii_io(argv[2], std::ios::out);
121
122   // Third argument: Maximal number of additional jets
123   int njet = atoi(argv[3]);
124
125   // Read input and output paths
126   string iPath = string(argv[4]);
127   string oPath = string(argv[5]);
128
129   // To write correctly normalized events to hepmc file, first get
130   // a reasonable accurate of the cross section
131   int njetCounterEstimate = njet;
132   vector<double> xsecEstimate;
133
134   vector<double> nTrialEstimate;
135   vector<double> nAcceptEstimate;
136
137   pythia.readString("Random:setSeed = on");
138   pythia.readString("Random:seed = 42390964");
139
140   while(njetCounterEstimate >= 0) {
141
142     // Number of runs
143     int nRun = 1;
144     double nTrial = 0.;
145     double nAccept = 0.;
146
147     int countEvents = 0;
148
149     // Run pythia nRun times with the same lhe file to get nRun times
150     // higher statistics in the histograms
151     for(int n = 1; n <= nRun ; ++n ) {
152
153       // Get process and events from LHE file, initialize only the
154       // first time
155       bool skipInit = false;
156       if(n > 1){
157         skipInit = true;
158         pythia.readString("Main:LHEFskipInit = on");
159       }
160
161       // From njet, choose LHE file
162       stringstream in;
163       in   << "_" << njetCounterEstimate << ".lhe";
164
165       string LHEfile = iPath + in.str();
166
167       pythia.readString("HadronLevel:all = off");
168
169       // Read in ME configurations
170       pythia.init(LHEfile,skipInit);
171
172       for( int iEvent=0; iEvent<nEvent; ++iEvent ){
173         countEvents++;
174
175         nTrial += 1.;
176         if(iEvent == 0) pythia.stat();
177
178         // Generate next event
179         if(pythia.next()) nAccept += 1.;
180
181         if(countEvents == nEvent*nRun-1){
182           xsecEstimate.push_back(pythia.info.sigmaGen());
183           nTrialEstimate.push_back(nTrial+1.);
184           nAcceptEstimate.push_back(nAccept+1.);
185         }
186
187
188       } // end loop over events to generate
189
190     } // end outer loop to rerun pythia with the same lhe file
191
192     // Restart with ME of a reduced the number of jets
193     if( njetCounterEstimate > 0 )
194       njetCounterEstimate--;
195     else
196       break;
197
198   } // end loop over different jet multiplicities
199
200   cout << endl << "Finished estimating cross section"
201     << endl;
202
203   for(int i=0; i < int(xsecEstimate.size()); ++i)
204     cout << "  Cross section estimate for " << njet-i << " jets :"
205       << scientific << setprecision(8) << xsecEstimate[i]
206       << endl;
207   for(int i=0; i < int(nTrialEstimate.size()); ++i)
208     cout << "  Trial events for " << njet-i << " jets :"
209       << scientific << setprecision(3) << nTrialEstimate[i]
210       << "  Accepted events for " << njet-i << " jets :"
211       << scientific << setprecision(3) << nAcceptEstimate[i]
212       << endl;
213   cout << endl;
214
215   // Now start merging procedure
216   int njetCounter = njet;
217
218   Hist histPTFirstSum("pT of first jet",100,0.,100.);
219   Hist histPTSecondSum("pT of second jet",100,0.,100.);
220
221   pythia.readString("Random:setSeed = on");
222   pythia.readString("Random:seed = 42390964");
223
224   // Sum of event weights
225   double sigma = 0.0;
226   double sigma2 = 0.0;
227
228   while(njetCounter >= 0) {
229
230     cout << "   Path to lhe files: " << iPath << "_*" << endl;
231     cout << "   Output written to: " << oPath << "'name'.dat" << endl;
232
233     // Set up histograms of pT of the first jet
234     Hist histPTFirst("pT of first jet",100,0.,200.);
235     Hist histPTSecond("pT of second jet",100,0.,200.);
236     Hist histPTThird("pT of third jet",100,0.,200.);
237     Hist histPTFourth("pT of fourth jet",50,0.,100.);
238     Hist histPTFifth("pT of fifth jet",30,0.,50.);
239     Hist histPTSixth("pT of sixth jet",30,0.,50.);
240
241     // Number of runs
242     int nRun = 1;
243     // Number of tried events
244     int nTriedEvents = 0;
245     // Number of accepted events
246     int nAccepEvents = 0;
247
248     // Run pythia nRun times with the same lhe file to get nRun times
249     // higher statistics in the histograms
250     for(int n = 1; n <= nRun ; ++n ) {
251
252       // Get process and events from LHE file, initialize only the
253       // first time
254       bool skipInit = false;
255       if(n > 1){
256         skipInit = true;
257         pythia.readString("Main:LHEFskipInit = on");
258       }
259
260       // From njet, choose LHE file
261       stringstream in;
262       in   << "_" << njetCounter << ".lhe";
263
264       string LHEfile = iPath + in.str();
265
266       cout << endl << endl
267         << "\t LHE FILE FOR + " << njetCounter
268         << " JET SAMPLE READ FROM " << LHEfile
269         << endl << endl;
270
271       cout << "Normalise with xsection " << xsecEstimate[njet-njetCounter]
272         << endl << endl;
273
274       pythia.readString("HadronLevel:all = on");
275
276       // Read in ME configurations
277       pythia.init(LHEfile,skipInit);
278
279       if(pythia.flag("Main:showChangedSettings")) {
280         pythia.settings.listChanged();
281       }
282
283       for( int iEvent=0; iEvent<nEvent; ++iEvent ){
284
285         nTriedEvents++;
286         if(iEvent == 0) pythia.stat();
287
288         // Generate next event
289         if( pythia.next()) {
290
291           double weight = pythia.info.mergingWeight();
292           nAccepEvents++;
293
294           // Jet pT's
295           double D = 0.4;
296           double pTfirst = pTfirstJet(pythia.event,1, D);
297           double pTsecnd = pTfirstJet(pythia.event,2, D);
298           double pTthird = pTfirstJet(pythia.event,3, D);
299           double pTfourt = pTfirstJet(pythia.event,4, D);
300           double pTfifth = pTfirstJet(pythia.event,5, D);
301           double pTsixth = pTfirstJet(pythia.event,6, D);
302           histPTFirst.fill( pTfirst, weight);
303           histPTSecond.fill( pTsecnd, weight);
304           histPTThird.fill( pTthird, weight);
305           histPTFourth.fill( pTfourt, weight);
306           histPTFifth.fill( pTfifth, weight);
307           histPTSixth.fill( pTsixth, weight);
308
309           if(weight > 0.){
310             // Construct new empty HepMC event. Form with arguments is only
311             // meaningful for HepMC 2.04 onwards, and even then unnecessary  
312             // if HepMC was built with GeV and mm as units from the onset. 
313             HepMC::GenEvent* hepmcevt = new HepMC::GenEvent();
314             //HepMC::GenEvent* hepmcevt = new HepMC::GenEvent(
315             //  HepMC::Units::GEV, HepMC::Units::MM); 
316
317             double normhepmc = 1.* xsecEstimate[njet-njetCounter]
318                 * nTrialEstimate[njet-njetCounter]
319                 / nAcceptEstimate[njet-njetCounter]
320                 * 1./ (double(nRun)*double(nEvent));
321
322             sigma += weight*normhepmc;
323             sigma2 += pow(weight*normhepmc,2);
324             // Set event weight
325             hepmcevt->weights().push_back(weight*normhepmc);
326
327             // Fill summed histograms
328             histPTFirstSum.fill( pTfirst, weight*normhepmc);
329             histPTSecondSum.fill( pTsecnd, weight*normhepmc);
330
331             // Fill HepMC event, including PDF info.
332             // ToHepMC.fill_next_event( pythia, hepmcevt );
333             // This alternative older method fills event, without PDF info.
334             ToHepMC.fill_next_event( pythia.event, hepmcevt );
335
336             // Report cross section to hepmc
337             HepMC::GenCrossSection xsec;
338             xsec.set_cross_section( sigma*1e9,
339               pythia.info.sigmaErr()*1e9 );
340             hepmcevt->set_cross_section( xsec );
341
342             // Write the HepMC event to file. Done with it.
343             ascii_io << hepmcevt;
344             delete hepmcevt;
345           }
346
347         } // if( pythia.next() )
348
349         if(nTriedEvents%10000 == 0)
350           cout << nTriedEvents << endl;
351
352       } // end loop over events to generate
353
354       // print cross section, errors
355       pythia.stat();
356
357     } // end outer loop to rerun pythia with the same lhe file
358
359     // Normalise histograms for this particular multiplicity
360     double norm = 1.
361                 * pythia.info.sigmaGen()
362                 * double(nTriedEvents)/double(nAccepEvents)
363                 * 1./ (double(nRun)*double(nEvent));
364
365     histPTFirst           *= norm;
366     histPTSecond          *= norm;
367     histPTThird           *= norm;
368     histPTFourth          *= norm;
369     histPTFifth           *= norm;
370     histPTSixth           *= norm;
371
372     // Write histograms for this particular multiplicity to file
373     ofstream write;
374     stringstream suffix;
375     suffix << njet << "_" << njetCounter;
376     suffix << "_wv.dat";
377
378     write.open( (char*)(oPath + "PTjet1_" + suffix.str()).c_str());
379     histPTFirst.table(write);
380     write.close();
381
382     write.open( (char*)(oPath + "PTjet2_" + suffix.str()).c_str());
383     histPTSecond.table(write);
384     write.close();
385
386     write.open( (char*)(oPath + "PTjet3_" + suffix.str()).c_str());
387     histPTThird.table(write);
388     write.close();
389
390     write.open( (char*)(oPath + "PTjet4_" + suffix.str()).c_str());
391     histPTFourth.table(write);
392     write.close();
393
394     write.open( (char*)(oPath + "PTjet5_" + suffix.str()).c_str());
395     histPTFifth.table(write);
396     write.close();
397
398     write.open( (char*)(oPath + "PTjet6_" + suffix.str()).c_str());
399     histPTSixth.table(write);
400     write.close();
401
402     histPTFirst.null();
403     histPTSecond.null();
404     histPTThird.null();
405     histPTFourth.null();
406     histPTFifth.null();
407     histPTSixth.null();
408
409     // Restart with ME of a reduced the number of jets
410     if( njetCounter > 0 )
411       njetCounter--;
412     else
413       break;
414
415   } // end loop over different jet multiplicities
416
417   // Since the histograms have been filled with the correct weight for
418   // each jet multiplicity, no normalisation is needed.
419   // Write summed histograms to file.
420   ofstream writeSum;
421   stringstream suffixSum;
422   suffixSum << njet << "_wv.dat";
423
424   writeSum.open( (char*)(oPath + "PTjet1Sum_" + suffixSum.str()).c_str());
425   histPTFirstSum.table(writeSum);
426   writeSum.close();
427
428   writeSum.open( (char*)(oPath + "PTjet2Sum_" + suffixSum.str()).c_str());
429   histPTSecondSum.table(writeSum);
430   writeSum.close();
431
432   for(int i=0; i < int(xsecEstimate.size()); ++i)
433     cout << "  Cross section estimate for " << njet-i << " jets :"
434       << scientific << setprecision(8) << xsecEstimate[i]
435       << endl;
436   for(int i=0; i < int(nTrialEstimate.size()); ++i)
437     cout << "  Trial events for " << njet-i << " jets :"
438       << scientific << setprecision(3) << nTrialEstimate[i]
439       << "  Accepted events for " << njet-i << " jets :"
440       << scientific << setprecision(3) << nAcceptEstimate[i]
441       << endl;
442   cout << endl;
443
444   cout << "Histogrammed cross section for "
445      << iPath << " with " << njet << " additional jets is " 
446      << scientific << setprecision(8) << sigma
447      << " error " << sqrt(sigma2) << endl;
448
449   return 0;
450 }