Fork me on GitHub

source: git/readers/DelphesProIO.cpp@ e1e72fd

Last change on this file since e1e72fd 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: 12.4 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 * @author S.Chekanov (ANL)
18 */
19
20#include <iostream>
21#include <memory>
22#include <sstream>
23#include <stdexcept>
24
25#include <map>
26
27#include <signal.h>
28#include <stdio.h>
29#include <stdlib.h>
30
31#include "TApplication.h"
32#include "TROOT.h"
33
34#include "TDatabasePDG.h"
35#include "TFile.h"
36#include "TLorentzVector.h"
37#include "TObjArray.h"
38#include "TParticlePDG.h"
39#include "TStopwatch.h"
40
41#include "classes/DelphesClasses.h"
42#include "classes/DelphesFactory.h"
43#include "classes/DelphesStream.h"
44#include "modules/Delphes.h"
45
46#include "ExRootAnalysis/ExRootProgressBar.h"
47#include "ExRootAnalysis/ExRootTreeBranch.h"
48#include "ExRootAnalysis/ExRootTreeWriter.h"
49
50#include <proio/event.h>
51#include <proio/model/mc.pb.h>
52#include <proio/reader.h>
53namespace model = proio::model::mc;
54
55using namespace std;
56
57//---------------------------------------------------------------------------
58// This method dynamically checks the message type (varint or not) depending on
59// non-zero value of units momentumUnit and positionUnit.
60void ConvertInput(proio::Event *event, double momentumUnit, double positionUnit,
61 ExRootTreeBranch *branch, DelphesFactory *factory,
62 TObjArray *allParticleOutputArray, TObjArray *stableParticleOutputArray,
63 TObjArray *partonOutputArray, TStopwatch *readStopWatch, TStopwatch *procStopWatch)
64{
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 element = static_cast<HepMCEvent *>(branch->NewEntry());
80
81 int nID = 0;
82 double weight = 0;
83 int process_id = 0;
84 auto mcentries = event->TaggedEntries("MCParameters");
85 for(uint64_t mcentryID : mcentries)
86 {
87 auto mcpar = dynamic_cast<proio::model::mc::MCParameters *>(event->GetEntry(mcentryID));
88 nID = mcpar->number();
89 weight = mcpar->weight();
90 process_id = mcpar->processid();
91 break; // consider only most generic from 1st entry
92 };
93
94 element->Number = nID;
95 element->ProcessID = process_id;
96 element->Weight = weight;
97
98 /*
99 // Pythia8 specific
100 element->MPI = mutableEvent->mpi();
101 element->Scale = mutableEvent->scale();
102 element->AlphaQED = mutableEvent->alpha_qed();
103 element->AlphaQCD = mutableEvent->alpha_qcd();
104 element->ID1 = mutableEvent->id1();
105 element->ID2 = mutableEvent->id2();
106 element->X1 = mutableEvent->x1();
107 element->X2 = mutableEvent->x2();
108 element->ScalePDF = mutableEvent->scale_pdf();
109 element->PDF1 = mutableEvent->pdf1();
110 element->PDF2 = mutableEvent->pdf2();
111*/
112
113 element->ReadTime = readStopWatch->RealTime();
114 element->ProcTime = procStopWatch->RealTime();
115
116 if(momentumUnit > 0 && positionUnit > 0)
117 {
118
119 auto entries = event->TaggedEntries("VarintPackedParticles");
120
121 for(uint64_t entryID : entries)
122 {
123
124 auto mutableParticles = dynamic_cast<model::VarintPackedParticles *>(event->GetEntry(entryID));
125
126 for(int i = 0; i < mutableParticles->pdg_size(); ++i)
127 {
128 pid = mutableParticles->pdg(i);
129 status = mutableParticles->status(i);
130 px = mutableParticles->px(i) / momentumUnit;
131 py = mutableParticles->py(i) / momentumUnit;
132 pz = mutableParticles->pz(i) / momentumUnit;
133 mass = mutableParticles->mass(i) / momentumUnit;
134 x = mutableParticles->x(i) / positionUnit;
135 y = mutableParticles->y(i) / positionUnit;
136 z = mutableParticles->z(i) / positionUnit;
137 t = mutableParticles->t(i) / positionUnit;
138
139 candidate = factory->NewCandidate();
140 candidate->PID = pid;
141 pdgCode = TMath::Abs(candidate->PID);
142 candidate->Status = status;
143 candidate->M1 = mutableParticles->parent1(i);
144 candidate->M2 = mutableParticles->parent2(i);
145 candidate->D1 = mutableParticles->child1(i);
146 candidate->D2 = mutableParticles->child2(i);
147 pdgParticle = pdg->GetParticle(pid);
148 candidate->Charge = mutableParticles->charge(i) / 3.0;
149 //candidate->Charge = pdgParticle ? Int_t(pdgParticle->Charge()/3.0) : -999;
150 candidate->Mass = mass;
151 candidate->Momentum.SetXYZM(px, py, pz, mass);
152 candidate->Position.SetXYZT(x, y, z, t);
153 allParticleOutputArray->Add(candidate);
154 if(!pdgParticle) continue;
155
156 if(status == 1)
157 {
158 stableParticleOutputArray->Add(candidate);
159 }
160 else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
161 {
162 partonOutputArray->Add(candidate);
163 }
164 }
165 }
166 }
167 else
168 {
169
170 auto entries = event->TaggedEntries("Particle");
171
172 for(uint64_t entryID : entries)
173 {
174 auto part = dynamic_cast<model::Particle *>(event->GetEntry(entryID));
175 pid = part->pdg();
176 status = part->status();
177 model::XYZF pvector = part->p();
178 px = pvector.x();
179 py = pvector.y();
180 pz = pvector.z();
181 mass = part->mass();
182 auto v = part->vertex();
183 x = v.x();
184 y = v.y();
185 z = v.z();
186 t = v.t();
187
188 candidate = factory->NewCandidate();
189 candidate->PID = pid;
190 pdgCode = TMath::Abs(candidate->PID);
191 candidate->Status = status;
192
193 int M1 = 0;
194 int M2 = 0;
195 for(int k1 = 0; k1 < part->parent_size(); k1++)
196 {
197 if(k1 == 0)
198 {
199 auto mother = dynamic_cast<model::Particle *>(event->GetEntry(part->parent(0)));
200 M1 = mother->barcode();
201 }
202 if(k1 == 1)
203 {
204 auto mother = dynamic_cast<model::Particle *>(event->GetEntry(part->parent(1)));
205 M2 = mother->barcode();
206 }
207 }
208
209 int D1 = 0;
210 int D2 = 0;
211 for(int k1 = 0; k1 < part->child_size(); k1++)
212 {
213 if(k1 == 0)
214 {
215 auto child = dynamic_cast<model::Particle *>(event->GetEntry(part->child(0)));
216 D1 = child->barcode();
217 }
218
219 if(k1 == 1)
220 {
221 auto child = dynamic_cast<model::Particle *>(event->GetEntry(part->child(1)));
222 D2 = child->barcode();
223 };
224 }
225
226 candidate->M1 = M1;
227 candidate->M2 = M2;
228 candidate->D1 = D1;
229 candidate->D2 = D2;
230
231 pdgParticle = pdg->GetParticle(pid);
232 candidate->Charge = pdgParticle ? Int_t(pdgParticle->Charge() / 3.0) : -999;
233 candidate->Mass = mass;
234
235 candidate->Momentum.SetXYZM(px, py, pz, mass);
236
237 candidate->Position.SetXYZT(x, y, z, t);
238
239 allParticleOutputArray->Add(candidate);
240
241 if(!pdgParticle) continue;
242
243 if(status == 1)
244 {
245 stableParticleOutputArray->Add(candidate);
246 }
247 else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
248 {
249 partonOutputArray->Add(candidate);
250 }
251 }
252
253 } // end particle type
254}
255
256//---------------------------------------------------------------------------
257
258static bool interrupted = false;
259
260void SignalHandler(int sig)
261{
262 interrupted = true;
263}
264
265//---------------------------------------------------------------------------
266
267int main(int argc, char *argv[])
268{
269 char appName[] = "DelphesProIO";
270 stringstream message;
271 proio::Reader *inputFile = 0;
272 TFile *outputFile = 0;
273 TStopwatch readStopWatch, procStopWatch;
274 ExRootTreeWriter *treeWriter = 0;
275 ExRootTreeBranch *branchEvent = 0;
276 ExRootConfReader *confReader = 0;
277 Delphes *modularDelphes = 0;
278 DelphesFactory *factory = 0;
279 TObjArray *allParticleOutputArray = 0, *stableParticleOutputArray = 0, *partonOutputArray = 0;
280 Int_t i;
281 Long64_t eventCounter, numberOfEvents;
282
283 if(argc < 4)
284 {
285 cout << " Usage: " << appName << " config_file"
286 << " output_file"
287 << " input_file(s)" << endl;
288 cout << " config_file - configuration file in Tcl format," << endl;
289 cout << " output_file - output file in ROOT format," << endl;
290 cout << " input_file(s) - input file(s) in ProIO format." << endl;
291 return 1;
292 }
293
294 signal(SIGINT, SignalHandler);
295
296 gROOT->SetBatch();
297
298 int appargc = 1;
299 char *appargv[] = {appName};
300 TApplication app(appName, &appargc, appargv);
301
302 try
303 {
304 outputFile = TFile::Open(argv[2], "CREATE");
305
306 if(outputFile == NULL)
307 {
308 message << "can't open " << argv[2] << endl;
309 throw runtime_error(message.str());
310 }
311
312 treeWriter = new ExRootTreeWriter(outputFile, "Delphes");
313
314 branchEvent = treeWriter->NewBranch("Event", HepMCEvent::Class());
315
316 confReader = new ExRootConfReader;
317 confReader->ReadFile(argv[1]);
318
319 modularDelphes = new Delphes("Delphes");
320 modularDelphes->SetConfReader(confReader);
321 modularDelphes->SetTreeWriter(treeWriter);
322
323 factory = modularDelphes->GetFactory();
324 allParticleOutputArray = modularDelphes->ExportArray("allParticles");
325 stableParticleOutputArray = modularDelphes->ExportArray("stableParticles");
326 partonOutputArray = modularDelphes->ExportArray("partons");
327
328 modularDelphes->InitTask();
329
330 for(i = 3; i < argc && !interrupted; ++i)
331 {
332 cout << "** INFO: Reading " << argv[i] << endl;
333
334 inputFile = new proio::Reader(argv[i]);
335
336 if(inputFile == NULL)
337 {
338 message << "can't open " << argv[i] << endl;
339 throw runtime_error(message.str());
340 }
341
342 /*
343// this is slow method, but general
344 inputFile->SeekToStart();
345 int nn=0;
346 auto event = new proio::Event();
347 while(inputFile->Next(event)){
348 auto entries = event->TaggedEntries("Particle");
349 if (entries.size()>0) nn++;
350 }
351*/
352
353 auto event = new proio::Event();
354
355 double varint_energy = 0;
356 double varint_length = 0;
357
358 auto max_n_events = std::numeric_limits<uint64_t>::max();
359 auto nn = inputFile->Skip(max_n_events);
360 cout << "** INFO: " << nn - 1 << " events found in ProIO file" << endl;
361 inputFile->SeekToStart();
362 numberOfEvents = nn - 1; // last event has only metadata (no particle record)
363
364 if(numberOfEvents <= 0) continue;
365
366 ExRootProgressBar progressBar(numberOfEvents - 1);
367
368 // Loop over all objects
369 modularDelphes->Clear();
370 treeWriter->Clear();
371 readStopWatch.Start();
372
373 for(eventCounter = 0; eventCounter < numberOfEvents && !interrupted; ++eventCounter)
374 {
375 inputFile->Next(event);
376 if(event == 0) continue;
377
378 // get metadata
379 if(eventCounter == 0)
380 {
381 auto metadata = event->Metadata();
382 std::cout << "** INFO: ProIO file metadata:" << std::endl;
383 for(auto element : metadata)
384 {
385 string key = (string)element.first;
386 string value = (string)(*element.second);
387 std::cout << "** INFO: " << key << " = " << value << std::endl;
388 if(key == "info:varint_energy") varint_energy = std::stod(value);
389 if(key == "info:varint_length") varint_length = std::stod(value);
390 }
391 }
392
393 readStopWatch.Stop();
394
395 procStopWatch.Start();
396
397 ConvertInput(event, varint_energy, varint_length,
398 branchEvent, factory,
399 allParticleOutputArray, stableParticleOutputArray,
400 partonOutputArray, &readStopWatch, &procStopWatch);
401
402 modularDelphes->ProcessTask();
403 procStopWatch.Stop();
404
405 treeWriter->Fill();
406
407 modularDelphes->Clear();
408 treeWriter->Clear();
409
410 readStopWatch.Start();
411 progressBar.Update(eventCounter);
412 }
413
414 progressBar.Update(eventCounter, eventCounter, kTRUE);
415 progressBar.Finish();
416
417 delete inputFile;
418 }
419
420 modularDelphes->FinishTask();
421 treeWriter->Write();
422
423 cout << "** Exiting..." << endl;
424
425 delete modularDelphes;
426 delete confReader;
427 delete treeWriter;
428 return 0;
429 }
430 catch(runtime_error &e)
431 {
432 if(treeWriter) delete treeWriter;
433 if(outputFile) delete outputFile;
434 cerr << "** ERROR: " << e.what() << endl;
435 return 1;
436 }
437}
Note: See TracBrowser for help on using the repository browser.