Fork me on GitHub

source: git/readers/DelphesHepMC3.cpp@ 95a917c

Last change on this file since 95a917c was 95a917c, checked in by Pavel Demin <pavel.demin@…>, 3 years ago

add HepMC3 library

  • Property mode set to 100644
File size: 14.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 <iostream>
20#include <fstream>
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 "modules/Delphes.h"
44
45#include "ExRootAnalysis/ExRootProgressBar.h"
46#include "ExRootAnalysis/ExRootTreeBranch.h"
47#include "ExRootAnalysis/ExRootTreeWriter.h"
48
49#include "HepMC3/ReaderAscii.h"
50#include "HepMC3/GenEvent.h"
51#include "HepMC3/GenCrossSection.h"
52#include "HepMC3/GenPdfInfo.h"
53#include "HepMC3/GenParticle.h"
54#include "HepMC3/GenVertex.h"
55#include "HepMC3/Units.h"
56
57using namespace std;
58using namespace HepMC3;
59
60map<Int_t, pair<Int_t, Int_t> > gMotherMap;
61map<Int_t, pair<Int_t, Int_t> > gDaughterMap;
62
63//---------------------------------------------------------------------------
64
65void AnalyzeParticle(Bool_t in, Int_t counter,
66 Double_t momentumCoefficient,
67 Double_t positionCoefficient,
68 shared_ptr<HepMC3::GenVertex> vertex,
69 shared_ptr<HepMC3::GenParticle> particle,
70 DelphesFactory *factory,
71 TObjArray *allParticleOutputArray,
72 TObjArray *stableParticleOutputArray,
73 TObjArray *partonOutputArray)
74{
75 map<Int_t, pair<Int_t, Int_t> >::iterator itMotherMap;
76 map<Int_t, pair<Int_t, Int_t> >::iterator itDaughterMap;
77
78 Candidate *candidate;
79 TDatabasePDG *pdg;
80 TParticlePDG *pdgParticle;
81 Int_t pdgCode;
82
83 Int_t pid, status, inVertexCode, outVertexCode;
84 Double_t px, py, pz, e, mass;
85 Double_t x, y, z, t;
86
87 pdg = TDatabasePDG::Instance();
88
89 candidate = factory->NewCandidate();
90
91 pid = particle->pid();
92 px = particle->momentum().px();
93 py = particle->momentum().py();
94 pz = particle->momentum().pz();
95 e = particle->momentum().e();
96 mass = particle->generated_mass();
97 x = vertex->position().x();
98 y = vertex->position().y();
99 z = vertex->position().z();
100 t = vertex->position().t();
101 status = particle->status();
102
103 outVertexCode = vertex->id();
104 inVertexCode = particle->end_vertex() ? particle->end_vertex()->id() : 0;
105
106 candidate->PID = pid;
107 pdgCode = TMath::Abs(pid);
108
109 candidate->Status = status;
110
111 pdgParticle = pdg->GetParticle(pid);
112 candidate->Charge = pdgParticle ? int(pdgParticle->Charge() / 3.0) : -999;
113 candidate->Mass = mass;
114
115 candidate->Momentum.SetPxPyPzE(px, py, pz, e);
116 if(momentumCoefficient != 1.0)
117 {
118 candidate->Momentum *= momentumCoefficient;
119 }
120
121 candidate->M2 = 1;
122 candidate->D2 = 1;
123
124 if(in)
125 {
126 candidate->M1 = 1;
127 candidate->Position.SetXYZT(0.0, 0.0, 0.0, 0.0);
128 }
129 else
130 {
131 candidate->M1 = outVertexCode;
132 candidate->Position.SetXYZT(x, y, z, t);
133 if(positionCoefficient != 1.0)
134 {
135 candidate->Position *= positionCoefficient;
136 }
137
138 itDaughterMap = gDaughterMap.find(outVertexCode);
139 if(itDaughterMap == gDaughterMap.end())
140 {
141 gDaughterMap[outVertexCode] = make_pair(counter, counter);
142 }
143 else
144 {
145 itDaughterMap->second.second = counter;
146 }
147 }
148
149 if(inVertexCode < 0)
150 {
151 candidate->D1 = inVertexCode;
152
153 itMotherMap = gMotherMap.find(inVertexCode);
154 if(itMotherMap == gMotherMap.end())
155 {
156 gMotherMap[inVertexCode] = make_pair(counter, -1);
157 }
158 else
159 {
160 itMotherMap->second.second = counter;
161 }
162 }
163 else
164 {
165 candidate->D1 = 1;
166 }
167
168 allParticleOutputArray->Add(candidate);
169
170 if(!pdgParticle) return;
171
172 if(status == 1)
173 {
174 stableParticleOutputArray->Add(candidate);
175 }
176 else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
177 {
178 partonOutputArray->Add(candidate);
179 }
180}
181
182//---------------------------------------------------------------------------
183
184void AnalyzeEvent(GenEvent &event, ExRootTreeBranch *branchEvent, DelphesFactory *factory,
185 TObjArray *allParticleOutputArray, TObjArray *stableParticleOutputArray,
186 TObjArray *partonOutputArray, TStopwatch *readStopWatch, TStopwatch *procStopWatch)
187{
188 Int_t i, counter;
189 map<Int_t, pair<Int_t, Int_t> >::iterator itMotherMap;
190 map<Int_t, pair<Int_t, Int_t> >::iterator itDaughterMap;
191
192 HepMCEvent *element;
193 Candidate *candidate;
194 Double_t momentumCoefficient, positionCoefficient;
195
196 shared_ptr<IntAttribute> processID = event.attribute<IntAttribute>("signal_process_id");
197 shared_ptr<IntAttribute> mpi = event.attribute<IntAttribute>("mpi");
198 shared_ptr<DoubleAttribute> scale = event.attribute<DoubleAttribute>("event_scale");
199 shared_ptr<DoubleAttribute> alphaQED = event.attribute<DoubleAttribute>("alphaQED");
200 shared_ptr<DoubleAttribute> alphaQCD = event.attribute<DoubleAttribute>("alphaQCD");
201
202 shared_ptr<GenCrossSection> cs = event.attribute<GenCrossSection>("GenCrossSection");
203 shared_ptr<GenPdfInfo> pdf = event.attribute<GenPdfInfo>("GenPdfInfo");
204
205 element = static_cast<HepMCEvent *>(branchEvent->NewEntry());
206
207 element->Number = event.event_number();
208
209 element->ProcessID = processID ? processID->value() : 0;
210 element->MPI = mpi ? mpi->value() : 0;
211 element->Weight = event.weights().size() > 0 ? event.weights()[0] : 1.0;
212 element->Scale = scale ? scale->value() : 0.0;
213 element->AlphaQED = alphaQED ? alphaQED->value() : 0.0;
214 element->AlphaQCD = alphaQCD ? alphaQCD->value() : 0.0;
215
216 if(cs)
217 {
218 element->CrossSection = cs->xsec();
219 element->CrossSectionError = cs->xsec_err();
220 }
221 else
222 {
223 element->CrossSection = 0.0;
224 element->CrossSectionError = 0.0;;
225 }
226
227 if(pdf)
228 {
229 element->ID1 = pdf->parton_id[0];
230 element->ID2 = pdf->parton_id[1];
231 element->X1 = pdf->x[0];
232 element->X2 = pdf->x[1];
233 element->ScalePDF = pdf->scale;
234 element->PDF1 = pdf->xf[0];
235 element->PDF2 = pdf->xf[1];
236 }
237 else
238 {
239 element->ID1 = 0;
240 element->ID2 = 0;
241 element->X1 = 0.0;
242 element->X2 = 0.0;
243 element->ScalePDF = 0.0;
244 element->PDF1 = 0.0;
245 element->PDF2 = 0.0;
246 }
247
248 element->ReadTime = readStopWatch->RealTime();
249 element->ProcTime = procStopWatch->RealTime();
250
251 momentumCoefficient = (event.momentum_unit() == Units::MEV) ? 0.001 : 1.0;
252 positionCoefficient = (event.length_unit() == Units::MM) ? 1.0 : 10.0;
253
254 counter = 0;
255 for(auto vertex: event.vertices())
256 {
257 for(auto particle: vertex->particles_in())
258 {
259 if(!particle->production_vertex() || particle->production_vertex()->id() == 0)
260 {
261 AnalyzeParticle(kTRUE, counter, momentumCoefficient, positionCoefficient, vertex, particle,
262 factory, allParticleOutputArray, stableParticleOutputArray, partonOutputArray);
263 ++counter;
264 }
265 }
266 for(auto particle: vertex->particles_out())
267 {
268 AnalyzeParticle(kFALSE, counter, momentumCoefficient, positionCoefficient, vertex, particle,
269 factory, allParticleOutputArray, stableParticleOutputArray, partonOutputArray);
270 ++counter;
271 }
272 }
273
274 for(i = 0; i < allParticleOutputArray->GetEntriesFast(); ++i)
275 {
276 candidate = static_cast<Candidate *>(allParticleOutputArray->At(i));
277
278 if(candidate->M1 > 0)
279 {
280 candidate->M1 = -1;
281 candidate->M2 = -1;
282 }
283 else
284 {
285 itMotherMap = gMotherMap.find(candidate->M1);
286 if(itMotherMap == gMotherMap.end())
287 {
288 candidate->M1 = -1;
289 candidate->M2 = -1;
290 }
291 else
292 {
293 candidate->M1 = itMotherMap->second.first;
294 candidate->M2 = itMotherMap->second.second;
295 }
296 }
297 if(candidate->D1 > 0)
298 {
299 candidate->D1 = -1;
300 candidate->D2 = -1;
301 }
302 else
303 {
304 itDaughterMap = gDaughterMap.find(candidate->D1);
305 if(itDaughterMap == gDaughterMap.end())
306 {
307 candidate->D1 = -1;
308 candidate->D2 = -1;
309 }
310 else
311 {
312 candidate->D1 = itDaughterMap->second.first;
313 candidate->D2 = itDaughterMap->second.second;
314 }
315 }
316 }
317}
318
319//---------------------------------------------------------------------------
320
321void AnalyzeWeight(GenEvent &event, ExRootTreeBranch *branchWeight)
322{
323 Weight *element;
324
325 for(auto weight: event.weights())
326 {
327 element = static_cast<Weight *>(branchWeight->NewEntry());
328
329 element->Weight = weight;
330 }
331}
332
333//---------------------------------------------------------------------------
334
335static bool interrupted = false;
336
337void SignalHandler(int sig)
338{
339 interrupted = true;
340}
341
342//---------------------------------------------------------------------------
343
344int main(int argc, char *argv[])
345{
346 char appName[] = "DelphesHepMC3";
347 stringstream message;
348 ifstream inputFile;
349 filebuf *inputBuffer;
350 TFile *outputFile = 0;
351 TStopwatch readStopWatch, procStopWatch;
352 ExRootTreeWriter *treeWriter = 0;
353 ExRootTreeBranch *branchEvent = 0, *branchWeight = 0;
354 ExRootConfReader *confReader = 0;
355 Delphes *modularDelphes = 0;
356 DelphesFactory *factory = 0;
357 TObjArray *allParticleOutputArray = 0, *stableParticleOutputArray = 0, *partonOutputArray = 0;
358 ReaderAscii *reader = 0;
359 GenEvent event;
360 Int_t i, maxEvents, skipEvents;
361 Long64_t length, eventCounter;
362
363 if(argc < 3)
364 {
365 cout << " Usage: " << appName << " config_file"
366 << " output_file"
367 << " [input_file(s)]" << endl;
368 cout << " config_file - configuration file in Tcl format," << endl;
369 cout << " output_file - output file in ROOT format," << endl;
370 cout << " input_file(s) - input file(s) in HepMC format," << endl;
371 cout << " with no input_file, or when input_file is -, read standard input." << endl;
372 return 1;
373 }
374
375 signal(SIGINT, SignalHandler);
376
377 gROOT->SetBatch();
378
379 int appargc = 1;
380 char *appargv[] = {appName};
381 TApplication app(appName, &appargc, appargv);
382
383 try
384 {
385 outputFile = TFile::Open(argv[2], "CREATE");
386
387 if(outputFile == NULL)
388 {
389 message << "can't create output file " << argv[2];
390 throw runtime_error(message.str());
391 }
392
393 treeWriter = new ExRootTreeWriter(outputFile, "Delphes");
394
395 branchEvent = treeWriter->NewBranch("Event", HepMCEvent::Class());
396 branchWeight = treeWriter->NewBranch("Weight", Weight::Class());
397
398 confReader = new ExRootConfReader;
399 confReader->ReadFile(argv[1]);
400
401 maxEvents = confReader->GetInt("::MaxEvents", 0);
402 skipEvents = confReader->GetInt("::SkipEvents", 0);
403
404 if(maxEvents < 0)
405 {
406 throw runtime_error("MaxEvents must be zero or positive");
407 }
408
409 if(skipEvents < 0)
410 {
411 throw runtime_error("SkipEvents must be zero or positive");
412 }
413
414 modularDelphes = new Delphes("Delphes");
415 modularDelphes->SetConfReader(confReader);
416 modularDelphes->SetTreeWriter(treeWriter);
417
418 factory = modularDelphes->GetFactory();
419 allParticleOutputArray = modularDelphes->ExportArray("allParticles");
420 stableParticleOutputArray = modularDelphes->ExportArray("stableParticles");
421 partonOutputArray = modularDelphes->ExportArray("partons");
422
423 inputBuffer = inputFile.rdbuf();
424 reader = new ReaderAscii(inputFile);
425
426 modularDelphes->InitTask();
427
428 i = 3;
429 do
430 {
431 if(interrupted) break;
432
433 if(i == argc || strncmp(argv[i], "-", 2) == 0)
434 {
435 cout << "** Reading standard input" << endl;
436 inputFile.ios::rdbuf(cin.rdbuf());
437 length = -1;
438 }
439 else
440 {
441 cout << "** Reading " << argv[i] << endl;
442 inputFile.ios::rdbuf(inputBuffer);
443 inputFile.open(argv[i]);
444
445 if(inputFile.fail())
446 {
447 message << "can't open " << argv[i];
448 throw runtime_error(message.str());
449 }
450
451 inputFile.seekg(0, ios::end);
452 length = inputFile.tellg();
453 inputFile.seekg(0, ios::beg);
454
455 if(length <= 0)
456 {
457 inputFile.close();
458 inputFile.clear();
459 ++i;
460 continue;
461 }
462 }
463
464 ExRootProgressBar progressBar(length);
465
466 reader->skip(skipEvents);
467 progressBar.Update(inputFile.tellg(), skipEvents);
468
469 eventCounter = skipEvents;
470 modularDelphes->Clear();
471 treeWriter->Clear();
472 event.clear();
473 gMotherMap.clear();
474 gDaughterMap.clear();
475 readStopWatch.Start();
476 reader->read_event(event);
477 while((maxEvents <= 0 || eventCounter - skipEvents < maxEvents) && !reader->failed() && !interrupted)
478 {
479 readStopWatch.Stop();
480
481 ++eventCounter;
482
483 procStopWatch.Start();
484 AnalyzeEvent(event, branchEvent, factory,
485 allParticleOutputArray, stableParticleOutputArray,
486 partonOutputArray, &readStopWatch, &procStopWatch);
487 modularDelphes->ProcessTask();
488 AnalyzeWeight(event, branchWeight);
489 procStopWatch.Stop();
490
491 treeWriter->Fill();
492
493 modularDelphes->Clear();
494 treeWriter->Clear();
495 event.clear();
496 gMotherMap.clear();
497 gDaughterMap.clear();
498
499 progressBar.Update(inputFile.tellg(), eventCounter);
500
501 readStopWatch.Start();
502 reader->read_event(event);
503 }
504
505 progressBar.Update(length, eventCounter, kTRUE);
506 progressBar.Finish();
507
508 if(length > 0) inputFile.close();
509
510 ++i;
511 } while(i < argc);
512
513 modularDelphes->FinishTask();
514 treeWriter->Write();
515
516 cout << "** Exiting..." << endl;
517
518 delete reader;
519 delete modularDelphes;
520 delete confReader;
521 delete treeWriter;
522 delete outputFile;
523
524 return 0;
525 }
526 catch(runtime_error &e)
527 {
528 if(treeWriter) delete treeWriter;
529 if(outputFile) delete outputFile;
530 cerr << "** ERROR: " << e.what() << endl;
531 return 1;
532 }
533}
Note: See TracBrowser for help on using the repository browser.