Fork me on GitHub

source: svn/trunk/readers/DelphesPythia8.cpp

Last change on this file was 1395, checked in by Pavel Demin, 10 years ago

synchronize with GitHub

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id Revision Date
File size: 7.3 KB
Line 
1#include <stdexcept>
2#include <iostream>
3#include <sstream>
4
5#include <signal.h>
6
7#include "Pythia8/Pythia.h"
8
9#include "TROOT.h"
10#include "TApplication.h"
11
12#include "TFile.h"
13#include "TObjArray.h"
14#include "TStopwatch.h"
15#include "TDatabasePDG.h"
16#include "TParticlePDG.h"
17#include "TLorentzVector.h"
18
19#include "modules/Delphes.h"
20#include "classes/DelphesClasses.h"
21#include "classes/DelphesFactory.h"
22
23#include "ExRootAnalysis/ExRootTreeWriter.h"
24#include "ExRootAnalysis/ExRootTreeBranch.h"
25#include "ExRootAnalysis/ExRootProgressBar.h"
26
27using namespace std;
28
29//---------------------------------------------------------------------------
30
31void ConvertInput(Long64_t eventCounter, Pythia8::Pythia *pythia,
32 ExRootTreeBranch *branch, DelphesFactory *factory,
33 TObjArray *allParticleOutputArray, TObjArray *stableParticleOutputArray, TObjArray *partonOutputArray,
34 TStopwatch *readStopWatch, TStopwatch *procStopWatch)
35{
36 int i;
37
38 HepMCEvent *element;
39 Candidate *candidate;
40 TDatabasePDG *pdg;
41 TParticlePDG *pdgParticle;
42 Int_t pdgCode;
43
44 Int_t pid, status;
45 Double_t px, py, pz, e, mass;
46 Double_t x, y, z, t;
47
48 // event information
49 element = static_cast<HepMCEvent *>(branch->NewEntry());
50
51 element->Number = eventCounter;
52
53 element->ProcessID = pythia->info.code();
54 element->MPI = 1;
55 element->Weight = pythia->info.weight();
56 element->Scale = pythia->info.QRen();
57 element->AlphaQED = pythia->info.alphaEM();
58 element->AlphaQCD = pythia->info.alphaS();
59
60 element->ID1 = pythia->info.id1();
61 element->ID2 = pythia->info.id2();
62 element->X1 = pythia->info.x1();
63 element->X2 = pythia->info.x2();
64 element->ScalePDF = pythia->info.QFac();
65 element->PDF1 = pythia->info.pdf1();
66 element->PDF2 = pythia->info.pdf2();
67
68 element->ReadTime = readStopWatch->RealTime();
69 element->ProcTime = procStopWatch->RealTime();
70
71 pdg = TDatabasePDG::Instance();
72
73 for(i = 0; i < pythia->event.size(); ++i)
74 {
75 Pythia8::Particle &particle = pythia->event[i];
76
77 pid = particle.id();
78 status = pythia->event.statusHepMC(i);
79 px = particle.px(); py = particle.py(); pz = particle.pz(); e = particle.e(); mass = particle.m();
80 x = particle.xProd(); y = particle.yProd(); z = particle.zProd(); t = particle.tProd();
81
82 candidate = factory->NewCandidate();
83
84 candidate->PID = pid;
85 pdgCode = TMath::Abs(candidate->PID);
86
87 candidate->Status = status;
88
89 candidate->M1 = particle.mother1() - 1;
90 candidate->M2 = particle.mother2() - 1;
91
92 candidate->D1 = particle.daughter1() - 1;
93 candidate->D2 = particle.daughter2() - 1;
94
95 pdgParticle = pdg->GetParticle(pid);
96 candidate->Charge = pdgParticle ? Int_t(pdgParticle->Charge()/3.0) : -999;
97 candidate->Mass = mass;
98
99 candidate->Momentum.SetPxPyPzE(px, py, pz, e);
100
101 candidate->Position.SetXYZT(x, y, z, t);
102
103 allParticleOutputArray->Add(candidate);
104
105 if(!pdgParticle) continue;
106
107 if(status == 1)
108 {
109 stableParticleOutputArray->Add(candidate);
110 }
111 else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
112 {
113 partonOutputArray->Add(candidate);
114 }
115 }
116}
117
118//---------------------------------------------------------------------------
119
120static bool interrupted = false;
121
122void SignalHandler(int sig)
123{
124 interrupted = true;
125}
126
127//---------------------------------------------------------------------------
128
129int main(int argc, char *argv[])
130{
131 char appName[] = "DelphesPythia8";
132 stringstream message;
133 TFile *outputFile = 0;
134 TStopwatch readStopWatch, procStopWatch;
135 ExRootTreeWriter *treeWriter = 0;
136 ExRootTreeBranch *branchEvent = 0;
137 ExRootConfReader *confReader = 0;
138 Delphes *modularDelphes = 0;
139 DelphesFactory *factory = 0;
140 TObjArray *stableParticleOutputArray = 0, *allParticleOutputArray = 0, *partonOutputArray = 0;
141 Long64_t eventCounter, errorCounter;
142 Long64_t numberOfEvents, timesAllowErrors;
143
144 Pythia8::Pythia *pythia = 0;
145
146 if(argc != 4)
147 {
148 cout << " Usage: " << appName << " config_file" << " pythia_card" << " output_file" << endl;
149 cout << " config_file - configuration file in Tcl format," << endl;
150 cout << " pythia_card - Pythia8 configuration file," << endl;
151 cout << " output_file - output file in ROOT format." << endl;
152 return 1;
153 }
154
155 signal(SIGINT, SignalHandler);
156
157 gROOT->SetBatch();
158
159 int appargc = 1;
160 char *appargv[] = {appName};
161 TApplication app(appName, &appargc, appargv);
162
163 try
164 {
165 outputFile = TFile::Open(argv[3], "CREATE");
166
167 if(outputFile == NULL)
168 {
169 message << "can't create output file " << argv[3];
170 throw runtime_error(message.str());
171 }
172
173 treeWriter = new ExRootTreeWriter(outputFile, "Delphes");
174
175 branchEvent = treeWriter->NewBranch("Event", HepMCEvent::Class());
176
177 confReader = new ExRootConfReader;
178 confReader->ReadFile(argv[1]);
179
180 modularDelphes = new Delphes("Delphes");
181 modularDelphes->SetConfReader(confReader);
182 modularDelphes->SetTreeWriter(treeWriter);
183
184 factory = modularDelphes->GetFactory();
185 allParticleOutputArray = modularDelphes->ExportArray("allParticles");
186 stableParticleOutputArray = modularDelphes->ExportArray("stableParticles");
187 partonOutputArray = modularDelphes->ExportArray("partons");
188
189 modularDelphes->InitTask();
190
191 // Initialize pythia
192 pythia = new Pythia8::Pythia;
193
194 if(pythia == NULL)
195 {
196 throw runtime_error("can't create Pythia instance");
197 }
198
199 // Read in commands from configuration file
200 pythia->readFile(argv[2]);
201
202 // Extract settings to be used in the main program
203 numberOfEvents = pythia->mode("Main:numberOfEvents");
204 timesAllowErrors = pythia->mode("Main:timesAllowErrors");
205
206 pythia->init();
207
208 // ExRootProgressBar progressBar(numberOfEvents - 1);
209 ExRootProgressBar progressBar(-1);
210
211 // Loop over all events
212 errorCounter = 0;
213 treeWriter->Clear();
214 modularDelphes->Clear();
215 readStopWatch.Start();
216 for(eventCounter = 0; eventCounter < numberOfEvents && !interrupted; ++eventCounter)
217 {
218 if(!pythia->next())
219 {
220 // If failure because reached end of file then exit event loop
221 if (pythia->info.atEndOfFile())
222 {
223 cerr << "Aborted since reached end of Les Houches Event File" << endl;
224 break;
225 }
226
227 // First few failures write off as "acceptable" errors, then quit
228 if (++errorCounter < timesAllowErrors) continue;
229 cerr << "Event generation aborted prematurely, owing to error!" << endl;
230 break;
231 }
232
233 readStopWatch.Stop();
234
235 procStopWatch.Start();
236 ConvertInput(eventCounter, pythia, branchEvent, factory,
237 allParticleOutputArray, stableParticleOutputArray, partonOutputArray,
238 &readStopWatch, &procStopWatch);
239 modularDelphes->ProcessTask();
240 procStopWatch.Stop();
241
242 treeWriter->Fill();
243
244 treeWriter->Clear();
245 modularDelphes->Clear();
246
247 readStopWatch.Start();
248 progressBar.Update(eventCounter, eventCounter);
249 }
250
251 progressBar.Update(eventCounter, eventCounter, kTRUE);
252 progressBar.Finish();
253
254 pythia->statistics();
255
256 modularDelphes->FinishTask();
257 treeWriter->Write();
258
259 cout << "** Exiting..." << endl;
260
261 delete pythia;
262 delete modularDelphes;
263 delete confReader;
264 delete treeWriter;
265 delete outputFile;
266
267 return 0;
268 }
269 catch(runtime_error &e)
270 {
271 if(treeWriter) delete treeWriter;
272 if(outputFile) delete outputFile;
273 cerr << "** ERROR: " << e.what() << endl;
274 return 1;
275 }
276}
Note: See TracBrowser for help on using the repository browser.