Fork me on GitHub

source: git/readers/DelphesCMSFWLite.cpp @ 341014c

ImprovedOutputFileTimingllp
Last change on this file since 341014c was 341014c, checked in by Pavel Demin <pavel-demin@…>, 18 months ago

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

  • Property mode set to 100644
File size: 13.2 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 <algorithm>
20#include <iostream>
21#include <memory>
22#include <sstream>
23#include <stdexcept>
24
25#include <map>
26#include <vector>
27
28#include <signal.h>
29#include <stdio.h>
30#include <stdlib.h>
31
32#include "TApplication.h"
33#include "TROOT.h"
34
35#include "TDatabasePDG.h"
36#include "TFile.h"
37#include "TLorentzVector.h"
38#include "TObjArray.h"
39#include "TParticlePDG.h"
40#include "TStopwatch.h"
41
42#include "classes/DelphesClasses.h"
43#include "classes/DelphesFactory.h"
44#include "classes/DelphesStream.h"
45#include "modules/Delphes.h"
46
47#include "ExRootAnalysis/ExRootProgressBar.h"
48#include "ExRootAnalysis/ExRootTreeBranch.h"
49#include "ExRootAnalysis/ExRootTreeWriter.h"
50
51#include "DataFormats/FWLite/interface/Event.h"
52#include "DataFormats/FWLite/interface/Handle.h"
53#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
54#include "DataFormats/PatCandidates/interface/PackedGenParticle.h"
55#include "FWCore/FWLite/interface/FWLiteEnabler.h"
56#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
57#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
58#include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h"
59#include "SimDataFormats/GeneratorProducts/interface/WeightsInfo.h"
60
61using namespace std;
62
63//---------------------------------------------------------------------------
64
65void ConvertInput(fwlite::Event &event, Long64_t eventCounter,
66  ExRootTreeBranch *branchEvent, ExRootTreeBranch *branchWeight,
67  DelphesFactory *factory, TObjArray *allParticleOutputArray,
68  TObjArray *stableParticleOutputArray, TObjArray *partonOutputArray, Bool_t firstEvent)
69{
70
71  fwlite::Handle<GenEventInfoProduct> handleGenEventInfo;
72  fwlite::Handle<LHEEventProduct> handleLHEEvent;
73  fwlite::Handle<vector<reco::GenParticle>> handleParticle;
74  fwlite::Handle<vector<pat::PackedGenParticle>> handlePackedParticle;
75
76  vector<reco::GenParticle>::const_iterator itParticle;
77  vector<pat::PackedGenParticle>::const_iterator itPackedParticle;
78
79  vector<const reco::Candidate *> vectorCandidate;
80  vector<const reco::Candidate *>::iterator itCandidate;
81
82  handleGenEventInfo.getByLabel(event, "generator");
83
84  if(!((handleLHEEvent.getBranchNameFor(event, "source")).empty()))
85  {
86    handleLHEEvent.getByLabel(event, "source");
87  }
88  else if(!((handleLHEEvent.getBranchNameFor(event, "externalLHEProducer")).empty()))
89  {
90    handleLHEEvent.getByLabel(event, "externalLHEProducer");
91  }
92  else if(firstEvent)
93  {
94    cout << "Wrong LHEEvent Label! Please, check the input file." << endl;
95  }
96
97  if(!((handleParticle.getBranchNameFor(event, "genParticles")).empty()))
98  {
99    handleParticle.getByLabel(event, "genParticles");
100  }
101  else if(!((handlePackedParticle.getBranchNameFor(event, "packedGenParticles")).empty()) && !((handleParticle.getBranchNameFor(event, "prunedGenParticles")).empty()))
102  {
103    handleParticle.getByLabel(event, "prunedGenParticles");
104    handlePackedParticle.getByLabel(event, "packedGenParticles");
105  }
106  else
107  {
108    std::cout << "Wrong GenParticle Label! Please, check the input file." << std::endl;
109    exit(-1);
110  }
111
112  Bool_t foundLHE = !((handleLHEEvent.getBranchNameFor(event, "source")).empty()) || !((handleLHEEvent.getBranchNameFor(event, "externalLHEProducer")).empty());
113  Bool_t isMiniAOD = !((handlePackedParticle.getBranchNameFor(event, "packedGenParticles")).empty()) && ((handleParticle.getBranchNameFor(event, "genParticles")).empty());
114
115  HepMCEvent *element;
116  Weight *weight;
117  Candidate *candidate;
118  TDatabasePDG *pdg;
119  TParticlePDG *pdgParticle;
120  Int_t pdgCode;
121
122  Int_t pid, status;
123  Double_t px, py, pz, e, mass;
124  Double_t x, y, z;
125
126  element = static_cast<HepMCEvent *>(branchEvent->NewEntry());
127
128  element->Number = eventCounter;
129
130  element->ProcessID = handleGenEventInfo->signalProcessID();
131  element->MPI = 1;
132  element->Weight = handleGenEventInfo->weight();
133  element->Scale = handleGenEventInfo->qScale();
134  element->AlphaQED = handleGenEventInfo->alphaQED();
135  element->AlphaQCD = handleGenEventInfo->alphaQCD();
136
137  element->ID1 = 0;
138  element->ID2 = 0;
139  element->X1 = 0.0;
140  element->X2 = 0.0;
141  element->ScalePDF = 0.0;
142  element->PDF1 = 0.0;
143  element->PDF2 = 0.0;
144
145  element->ReadTime = 0.0;
146  element->ProcTime = 0.0;
147
148  if(foundLHE)
149  {
150    const vector<gen::WeightsInfo> &vectorWeightsInfo = handleLHEEvent->weights();
151    vector<gen::WeightsInfo>::const_iterator itWeightsInfo;
152
153    for(itWeightsInfo = vectorWeightsInfo.begin(); itWeightsInfo != vectorWeightsInfo.end(); ++itWeightsInfo)
154    {
155      weight = static_cast<Weight *>(branchWeight->NewEntry());
156      weight->Weight = itWeightsInfo->wgt;
157    }
158  }
159
160  pdg = TDatabasePDG::Instance();
161
162  for(itParticle = handleParticle->begin(); itParticle != handleParticle->end(); ++itParticle)
163  {
164    const reco::GenParticle &particle = *itParticle;
165    if(!isMiniAOD || particle.status() != 1) vectorCandidate.push_back(&*itParticle);
166  }
167
168  for(itParticle = handleParticle->begin(); itParticle != handleParticle->end(); ++itParticle)
169  {
170    const reco::GenParticle &particle = *itParticle;
171
172    pid = particle.pdgId();
173    status = particle.status();
174    if(isMiniAOD && particle.status() == 1) continue;
175    px = particle.px();
176    py = particle.py();
177    pz = particle.pz();
178    e = particle.energy();
179    mass = particle.mass();
180    x = particle.vx();
181    y = particle.vy();
182    z = particle.vz();
183
184    candidate = factory->NewCandidate();
185
186    candidate->PID = pid;
187    pdgCode = TMath::Abs(candidate->PID);
188
189    candidate->Status = status;
190
191    if(particle.mother())
192    {
193      itCandidate = find(vectorCandidate.begin(), vectorCandidate.end(), particle.mother());
194      if(itCandidate != vectorCandidate.end()) candidate->M1 = distance(vectorCandidate.begin(), itCandidate);
195    }
196
197    itCandidate = find(vectorCandidate.begin(), vectorCandidate.end(), particle.daughter(0));
198    if(itCandidate != vectorCandidate.end()) candidate->D1 = distance(vectorCandidate.begin(), itCandidate);
199
200    itCandidate = find(vectorCandidate.begin(), vectorCandidate.end(), particle.daughter(particle.numberOfDaughters() - 1));
201    if(itCandidate != vectorCandidate.end()) candidate->D2 = distance(vectorCandidate.begin(), itCandidate);
202
203    pdgParticle = pdg->GetParticle(pid);
204    candidate->Charge = pdgParticle ? Int_t(pdgParticle->Charge() / 3.0) : -999;
205    candidate->Mass = mass;
206
207    candidate->Momentum.SetPxPyPzE(px, py, pz, e);
208
209    candidate->Position.SetXYZT(x * 10.0, y * 10.0, z * 10.0, 0.0);
210
211    allParticleOutputArray->Add(candidate);
212
213    if(!pdgParticle) continue;
214
215    if(status == 1)
216    {
217      // Prevent duplicated particle.
218      if(!isMiniAOD) stableParticleOutputArray->Add(candidate);
219    }
220    else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
221    {
222      partonOutputArray->Add(candidate);
223    }
224  }
225
226  if(!isMiniAOD) return;
227  // For MiniAOD sample,
228  // Only status==1 particles are saved to packedGenParticles.
229  for(itPackedParticle = handlePackedParticle->begin(); itPackedParticle != handlePackedParticle->end(); ++itPackedParticle)
230  {
231    vectorCandidate.push_back(&*itPackedParticle);
232  }
233
234  for(itPackedParticle = handlePackedParticle->begin(); itPackedParticle != handlePackedParticle->end(); ++itPackedParticle)
235  {
236    const pat::PackedGenParticle &particle = *itPackedParticle;
237
238    pid = particle.pdgId();
239    status = particle.status();
240    px = particle.px();
241    py = particle.py();
242    pz = particle.pz();
243    e = particle.energy();
244    mass = particle.mass();
245    x = particle.vx();
246    y = particle.vy();
247    z = particle.vz();
248
249    candidate = factory->NewCandidate();
250
251    candidate->PID = pid;
252    pdgCode = TMath::Abs(candidate->PID);
253
254    candidate->Status = status;
255
256    if(particle.mother(0))
257    {
258      itCandidate = find(vectorCandidate.begin(), vectorCandidate.end(), particle.mother(0));
259      if(itCandidate != vectorCandidate.end()) candidate->M1 = distance(vectorCandidate.begin(), itCandidate);
260    }
261
262    itCandidate = find(vectorCandidate.begin(), vectorCandidate.end(), particle.daughter(0));
263    if(itCandidate != vectorCandidate.end()) candidate->D1 = distance(vectorCandidate.begin(), itCandidate);
264
265    itCandidate = find(vectorCandidate.begin(), vectorCandidate.end(), particle.daughter(particle.numberOfDaughters() - 1));
266    if(itCandidate != vectorCandidate.end()) candidate->D2 = distance(vectorCandidate.begin(), itCandidate);
267
268    pdgParticle = pdg->GetParticle(pid);
269    candidate->Charge = pdgParticle ? Int_t(pdgParticle->Charge() / 3.0) : -999;
270    candidate->Mass = mass;
271
272    candidate->Momentum.SetPxPyPzE(px, py, pz, e);
273
274    candidate->Position.SetXYZT(x * 10.0, y * 10.0, z * 10.0, 0.0);
275
276    allParticleOutputArray->Add(candidate);
277
278    if(!pdgParticle) continue;
279
280    if(status == 1)
281    {
282      stableParticleOutputArray->Add(candidate);
283    }
284  }
285}
286
287//---------------------------------------------------------------------------
288
289static bool interrupted = false;
290
291void SignalHandler(int sig)
292{
293  interrupted = true;
294}
295
296//---------------------------------------------------------------------------
297
298int main(int argc, char *argv[])
299{
300  char appName[] = "DelphesCMSFWLite";
301  stringstream message;
302  TFile *inputFile = 0;
303  TFile *outputFile = 0;
304  TStopwatch eventStopWatch;
305  ExRootTreeWriter *treeWriter = 0;
306  ExRootTreeBranch *branchEvent = 0, *branchWeight = 0;
307  ExRootConfReader *confReader = 0;
308  Delphes *modularDelphes = 0;
309  DelphesFactory *factory = 0;
310  TObjArray *allParticleOutputArray = 0, *stableParticleOutputArray = 0, *partonOutputArray = 0;
311  Int_t i;
312  Long64_t eventCounter, numberOfEvents;
313  Bool_t firstEvent = kTRUE;
314
315  if(argc < 4)
316  {
317    cout << " Usage: " << appName << " config_file"
318         << " output_file"
319         << " input_file(s)" << endl;
320    cout << " config_file - configuration file in Tcl format," << endl;
321    cout << " output_file - output file in ROOT format," << endl;
322    cout << " input_file(s) - input file(s) in ROOT format." << endl;
323    return 1;
324  }
325
326  signal(SIGINT, SignalHandler);
327
328  gROOT->SetBatch();
329
330  int appargc = 1;
331  char *appargv[] = {appName};
332  TApplication app(appName, &appargc, appargv);
333
334  FWLiteEnabler::enable();
335
336  try
337  {
338    outputFile = TFile::Open(argv[2], "CREATE");
339
340    if(outputFile == NULL)
341    {
342      message << "can't open " << argv[2] << endl;
343      throw runtime_error(message.str());
344    }
345
346    treeWriter = new ExRootTreeWriter(outputFile, "Delphes");
347
348    branchEvent = treeWriter->NewBranch("Event", HepMCEvent::Class());
349    branchWeight = treeWriter->NewBranch("Weight", Weight::Class());
350
351    confReader = new ExRootConfReader;
352    confReader->ReadFile(argv[1]);
353
354    modularDelphes = new Delphes("Delphes");
355    modularDelphes->SetConfReader(confReader);
356    modularDelphes->SetTreeWriter(treeWriter);
357
358    factory = modularDelphes->GetFactory();
359    allParticleOutputArray = modularDelphes->ExportArray("allParticles");
360    stableParticleOutputArray = modularDelphes->ExportArray("stableParticles");
361    partonOutputArray = modularDelphes->ExportArray("partons");
362
363    modularDelphes->InitTask();
364
365    for(i = 3; i < argc && !interrupted; ++i)
366    {
367      cout << "** Reading " << argv[i] << endl;
368
369      inputFile = TFile::Open(argv[i]);
370
371      if(inputFile == NULL)
372      {
373        message << "can't open " << argv[i] << endl;
374        throw runtime_error(message.str());
375      }
376
377      fwlite::Event event(inputFile);
378
379      numberOfEvents = event.size();
380
381      if(numberOfEvents <= 0) continue;
382
383      // ExRootProgressBar progressBar(numberOfEvents - 1);
384      ExRootProgressBar progressBar(-1);
385
386      // Loop over all objects
387      eventCounter = 0;
388      modularDelphes->Clear();
389      treeWriter->Clear();
390
391      for(event.toBegin(); !event.atEnd() && !interrupted; ++event)
392      {
393        ConvertInput(event, eventCounter, branchEvent, branchWeight, factory,
394          allParticleOutputArray, stableParticleOutputArray, partonOutputArray, firstEvent);
395        modularDelphes->ProcessTask();
396
397        firstEvent = kFALSE;
398
399        treeWriter->Fill();
400
401        modularDelphes->Clear();
402        treeWriter->Clear();
403
404        progressBar.Update(eventCounter, eventCounter);
405        ++eventCounter;
406      }
407
408      progressBar.Update(eventCounter, eventCounter, kTRUE);
409      progressBar.Finish();
410
411      inputFile->Close();
412    }
413
414    modularDelphes->FinishTask();
415    treeWriter->Write();
416
417    cout << "** Exiting..." << endl;
418
419    delete modularDelphes;
420    delete confReader;
421    delete treeWriter;
422    delete outputFile;
423
424    return 0;
425  }
426  catch(runtime_error &e)
427  {
428    if(treeWriter) delete treeWriter;
429    if(outputFile) delete outputFile;
430    cerr << "** ERROR: " << e.what() << endl;
431    return 1;
432  }
433}
Note: See TracBrowser for help on using the repository browser.