Fork me on GitHub

source: git/readers/DelphesProMC.cpp@ 7c0fcd5

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 7c0fcd5 was 7c0fcd5, checked in by Pavel Demin <demin@…>, 10 years ago

delete duplicate license file and prepend GPLv3 header to all source code files

  • Property mode set to 100644
File size: 8.0 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#include <memory>
23
24#include <map>
25
26#include <stdlib.h>
27#include <signal.h>
28#include <stdio.h>
29
30#include "TROOT.h"
31#include "TApplication.h"
32
33#include "TFile.h"
34#include "TObjArray.h"
35#include "TStopwatch.h"
36#include "TDatabasePDG.h"
37#include "TParticlePDG.h"
38#include "TLorentzVector.h"
39
40#include "modules/Delphes.h"
41#include "classes/DelphesStream.h"
42#include "classes/DelphesClasses.h"
43#include "classes/DelphesFactory.h"
44
45#include "ExRootAnalysis/ExRootTreeWriter.h"
46#include "ExRootAnalysis/ExRootTreeBranch.h"
47#include "ExRootAnalysis/ExRootProgressBar.h"
48
49#include "ProMC/ProMC.pb.h"
50#include "ProMC/ProMCBook.h"
51#include "ProMC/ProMCHeader.pb.h"
52
53using namespace std;
54
55//---------------------------------------------------------------------------
56
57void ConvertInput(ProMCEvent &event, ExRootTreeBranch *branch, DelphesFactory *factory,
58 TObjArray *allParticleOutputArray, TObjArray *stableParticleOutputArray, TObjArray *partonOutputArray,
59 TStopwatch *readStopWatch, TStopwatch *procStopWatch)
60{
61 Int_t i;
62
63 ProMCEvent_Event *mutableEvent;
64 ProMCEvent_Particles *mutableParticles;
65
66 HepMCEvent *element;
67 Candidate *candidate;
68 TDatabasePDG *pdg;
69 TParticlePDG *pdgParticle;
70 Int_t pdgCode;
71
72 Int_t pid, status;
73 Double_t px, py, pz, mass;
74 Double_t x, y, z, t;
75
76 pdg = TDatabasePDG::Instance();
77
78 // event information
79 mutableEvent = event.mutable_event();
80
81 element = static_cast<HepMCEvent *>(branch->NewEntry());
82
83 element->Number = mutableEvent->number();
84
85 element->ProcessID = mutableEvent->process_id();
86 element->MPI = mutableEvent->mpi();
87 element->Weight = mutableEvent->weight();
88 element->Scale = mutableEvent->scale();
89 element->AlphaQED = mutableEvent->alpha_qed();
90 element->AlphaQCD = mutableEvent->alpha_qcd();
91
92 element->ID1 = mutableEvent->id1();
93 element->ID2 = mutableEvent->id2();
94 element->X1 = mutableEvent->x1();
95 element->X2 = mutableEvent->x2();
96 element->ScalePDF = mutableEvent->scale_pdf();
97 element->PDF1 = mutableEvent->pdf1();
98 element->PDF2 = mutableEvent->pdf2();
99
100 element->ReadTime = readStopWatch->RealTime();
101 element->ProcTime = procStopWatch->RealTime();
102
103 mutableParticles = event.mutable_particles();
104
105 for(i = 0; i < mutableParticles->pdg_id_size(); ++i)
106 {
107 pid = mutableParticles->pdg_id(i);
108 status = mutableParticles->status(i);
109 px = mutableParticles->px(i); py = mutableParticles->py(i); pz = mutableParticles->pz(i); mass = mutableParticles->mass(i);
110 x = mutableParticles->x(i); y = mutableParticles->y(i); z = mutableParticles->z(i); t = mutableParticles->t(i);
111
112 candidate = factory->NewCandidate();
113
114 candidate->PID = pid;
115 pdgCode = TMath::Abs(candidate->PID);
116
117 candidate->Status = status;
118
119 candidate->M1 = mutableParticles->mother1(i);
120 candidate->M2 = mutableParticles->mother2(i);
121
122 candidate->D1 = mutableParticles->daughter1(i);
123 candidate->D2 = mutableParticles->daughter2(i);
124
125 pdgParticle = pdg->GetParticle(pid);
126 candidate->Charge = pdgParticle ? Int_t(pdgParticle->Charge()/3.0) : -999;
127 candidate->Mass = mass;
128
129 candidate->Momentum.SetXYZM(px, py, pz, mass);
130
131 candidate->Position.SetXYZT(x, y, z, t);
132
133 allParticleOutputArray->Add(candidate);
134
135 if(!pdgParticle) continue;
136
137 if(status == 1)
138 {
139 stableParticleOutputArray->Add(candidate);
140 }
141 else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
142 {
143 partonOutputArray->Add(candidate);
144 }
145 }
146}
147
148//---------------------------------------------------------------------------
149
150static bool interrupted = false;
151
152void SignalHandler(int sig)
153{
154 interrupted = true;
155}
156
157//---------------------------------------------------------------------------
158
159int main(int argc, char *argv[])
160{
161 char appName[] = "DelphesProMC";
162 stringstream message;
163 ProMCBook *inputFile = 0;
164 TFile *outputFile = 0;
165 TStopwatch readStopWatch, procStopWatch;
166 ExRootTreeWriter *treeWriter = 0;
167 ExRootTreeBranch *branchEvent = 0;
168 ExRootConfReader *confReader = 0;
169 Delphes *modularDelphes = 0;
170 DelphesFactory *factory = 0;
171 TObjArray *allParticleOutputArray = 0, *stableParticleOutputArray = 0, *partonOutputArray = 0;
172 Int_t i;
173 Long64_t eventCounter, numberOfEvents;
174
175 if(argc < 4)
176 {
177 cout << " Usage: " << appName << " config_file" << " output_file" << " input_file(s)" << endl;
178 cout << " config_file - configuration file in Tcl format," << endl;
179 cout << " output_file - output file in ROOT format," << endl;
180 cout << " input_file(s) - input file(s) in ProMC format." << endl;
181 return 1;
182 }
183
184 signal(SIGINT, SignalHandler);
185
186 gROOT->SetBatch();
187
188 int appargc = 1;
189 char *appargv[] = {appName};
190 TApplication app(appName, &appargc, appargv);
191
192 try
193 {
194 outputFile = TFile::Open(argv[2], "CREATE");
195
196 if(outputFile == NULL)
197 {
198 message << "can't open " << argv[2] << endl;
199 throw runtime_error(message.str());
200 }
201
202 treeWriter = new ExRootTreeWriter(outputFile, "Delphes");
203
204 branchEvent = treeWriter->NewBranch("Event", HepMCEvent::Class());
205
206 confReader = new ExRootConfReader;
207 confReader->ReadFile(argv[1]);
208
209 modularDelphes = new Delphes("Delphes");
210 modularDelphes->SetConfReader(confReader);
211 modularDelphes->SetTreeWriter(treeWriter);
212
213 factory = modularDelphes->GetFactory();
214 allParticleOutputArray = modularDelphes->ExportArray("allParticles");
215 stableParticleOutputArray = modularDelphes->ExportArray("stableParticles");
216 partonOutputArray = modularDelphes->ExportArray("partons");
217
218 modularDelphes->InitTask();
219
220 for(i = 3; i < argc && !interrupted; ++i)
221 {
222 cout << "** Reading " << argv[i] << endl;
223
224 inputFile = new ProMCBook(argv[i], "r");
225
226 if(inputFile == NULL)
227 {
228 message << "can't open " << argv[i] << endl;
229 throw runtime_error(message.str());
230 }
231
232 numberOfEvents = inputFile->getEvents();
233
234 if(numberOfEvents <= 0) continue;
235
236 ExRootProgressBar progressBar(numberOfEvents - 1);
237
238 // Loop over all objects
239 modularDelphes->Clear();
240 treeWriter->Clear();
241 readStopWatch.Start();
242 for(eventCounter = 0; eventCounter < numberOfEvents && !interrupted; ++eventCounter)
243 {
244 if(inputFile->next() != 0) continue;
245 ProMCEvent event = inputFile->get();
246
247 readStopWatch.Stop();
248
249 procStopWatch.Start();
250 ConvertInput(event, branchEvent, factory,
251 allParticleOutputArray, stableParticleOutputArray, partonOutputArray,
252 &readStopWatch, &procStopWatch);
253 modularDelphes->ProcessTask();
254 procStopWatch.Stop();
255
256 treeWriter->Fill();
257
258 modularDelphes->Clear();
259 treeWriter->Clear();
260
261 readStopWatch.Start();
262 progressBar.Update(eventCounter);
263 }
264
265 progressBar.Update(eventCounter, eventCounter, kTRUE);
266 progressBar.Finish();
267
268 inputFile->close();
269 delete inputFile;
270 }
271
272 modularDelphes->FinishTask();
273 treeWriter->Write();
274
275 cout << "** Exiting..." << endl;
276
277 delete modularDelphes;
278 delete confReader;
279 delete treeWriter;
280 delete outputFile;
281
282 return 0;
283 }
284 catch(runtime_error &e)
285 {
286 if(treeWriter) delete treeWriter;
287 if(outputFile) delete outputFile;
288 cerr << "** ERROR: " << e.what() << endl;
289 return 1;
290 }
291}
Note: See TracBrowser for help on using the repository browser.