Fork me on GitHub

source: git/readers/DelphesProMC.cpp

Last change on this file was 341014c, checked in by Pavel Demin <pavel-demin@…>, 6 years ago

apply .clang-format to all .h, .cc and .cpp files

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