Fork me on GitHub

source: git/readers/DelphesPythia8.cpp@ 0f8390e

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 0f8390e was 1fa50c2, checked in by Pavel Demin <pavel.demin@…>, 10 years ago

fix GPLv3 header

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