Fork me on GitHub

source: git/readers/DelphesCMSFWLite.cpp

Last change on this file was 0c0c9af, checked in by GitHub <noreply@…>, 3 years ago

Fix eventCounter/skipEvents bug (#101)

  • Property mode set to 100644
File size: 13.8 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 if (pdgCode == 11 || pdgCode == 13) partonOutputArray->Add(candidate);
220 }
221 //else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
222 else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 11 || pdgCode == 13 || pdgCode == 15)
223 {
224 partonOutputArray->Add(candidate);
225 }
226 }
227
228 if(!isMiniAOD) return;
229 // For MiniAOD sample,
230 // Only status==1 particles are saved to packedGenParticles.
231 for(itPackedParticle = handlePackedParticle->begin(); itPackedParticle != handlePackedParticle->end(); ++itPackedParticle)
232 {
233 vectorCandidate.push_back(&*itPackedParticle);
234 }
235
236 for(itPackedParticle = handlePackedParticle->begin(); itPackedParticle != handlePackedParticle->end(); ++itPackedParticle)
237 {
238 const pat::PackedGenParticle &particle = *itPackedParticle;
239
240 pid = particle.pdgId();
241 status = particle.status();
242 px = particle.px();
243 py = particle.py();
244 pz = particle.pz();
245 e = particle.energy();
246 mass = particle.mass();
247 x = particle.vx();
248 y = particle.vy();
249 z = particle.vz();
250
251 candidate = factory->NewCandidate();
252
253 candidate->PID = pid;
254 pdgCode = TMath::Abs(candidate->PID);
255
256 candidate->Status = status;
257
258 if(particle.mother(0))
259 {
260 itCandidate = find(vectorCandidate.begin(), vectorCandidate.end(), particle.mother(0));
261 if(itCandidate != vectorCandidate.end()) candidate->M1 = distance(vectorCandidate.begin(), itCandidate);
262 }
263
264 itCandidate = find(vectorCandidate.begin(), vectorCandidate.end(), particle.daughter(0));
265 if(itCandidate != vectorCandidate.end()) candidate->D1 = distance(vectorCandidate.begin(), itCandidate);
266
267 itCandidate = find(vectorCandidate.begin(), vectorCandidate.end(), particle.daughter(particle.numberOfDaughters() - 1));
268 if(itCandidate != vectorCandidate.end()) candidate->D2 = distance(vectorCandidate.begin(), itCandidate);
269
270 pdgParticle = pdg->GetParticle(pid);
271 candidate->Charge = pdgParticle ? Int_t(pdgParticle->Charge() / 3.0) : -999;
272 candidate->Mass = mass;
273
274 candidate->Momentum.SetPxPyPzE(px, py, pz, e);
275
276 candidate->Position.SetXYZT(x * 10.0, y * 10.0, z * 10.0, 0.0);
277
278 allParticleOutputArray->Add(candidate);
279
280 if(!pdgParticle) continue;
281
282 if(status == 1)
283 {
284 stableParticleOutputArray->Add(candidate);
285 }
286 }
287}
288
289//---------------------------------------------------------------------------
290
291static bool interrupted = false;
292
293void SignalHandler(int sig)
294{
295 interrupted = true;
296}
297
298//---------------------------------------------------------------------------
299
300int main(int argc, char *argv[])
301{
302 char appName[] = "DelphesCMSFWLite";
303 stringstream message;
304 TFile *inputFile = 0;
305 TFile *outputFile = 0;
306 TStopwatch eventStopWatch;
307 ExRootTreeWriter *treeWriter = 0;
308 ExRootTreeBranch *branchEvent = 0, *branchWeight = 0;
309 ExRootConfReader *confReader = 0;
310 Delphes *modularDelphes = 0;
311 DelphesFactory *factory = 0;
312 TObjArray *allParticleOutputArray = 0, *stableParticleOutputArray = 0, *partonOutputArray = 0;
313 Int_t i, maxEvents, skipEvents;
314 Long64_t eventCounter, numberOfEvents;
315 Bool_t firstEvent = kTRUE;
316
317 if(argc < 4)
318 {
319 cout << " Usage: " << appName << " config_file"
320 << " output_file"
321 << " input_file(s)" << endl;
322 cout << " config_file - configuration file in Tcl format," << endl;
323 cout << " output_file - output file in ROOT format," << endl;
324 cout << " input_file(s) - input file(s) in ROOT format." << endl;
325 return 1;
326 }
327
328 signal(SIGINT, SignalHandler);
329
330 gROOT->SetBatch();
331
332 int appargc = 1;
333 char *appargv[] = {appName};
334 TApplication app(appName, &appargc, appargv);
335
336 FWLiteEnabler::enable();
337
338 try
339 {
340 outputFile = TFile::Open(argv[2], "CREATE");
341
342 if(outputFile == NULL)
343 {
344 message << "can't open " << argv[2] << endl;
345 throw runtime_error(message.str());
346 }
347
348 treeWriter = new ExRootTreeWriter(outputFile, "Delphes");
349
350 branchEvent = treeWriter->NewBranch("Event", HepMCEvent::Class());
351 branchWeight = treeWriter->NewBranch("Weight", Weight::Class());
352
353 confReader = new ExRootConfReader;
354 confReader->ReadFile(argv[1]);
355
356 maxEvents = confReader->GetInt("::MaxEvents", 0);
357 skipEvents = confReader->GetInt("::SkipEvents", 0);
358
359 if(maxEvents < 0)
360 {
361 throw runtime_error("MaxEvents must be zero or positive");
362 }
363
364 if(skipEvents < 0)
365 {
366 throw runtime_error("SkipEvents must be zero or positive");
367 }
368
369 modularDelphes = new Delphes("Delphes");
370 modularDelphes->SetConfReader(confReader);
371 modularDelphes->SetTreeWriter(treeWriter);
372
373 factory = modularDelphes->GetFactory();
374 allParticleOutputArray = modularDelphes->ExportArray("allParticles");
375 stableParticleOutputArray = modularDelphes->ExportArray("stableParticles");
376 partonOutputArray = modularDelphes->ExportArray("partons");
377
378 modularDelphes->InitTask();
379
380 for(i = 3; i < argc && !interrupted; ++i)
381 {
382 cout << "** Reading " << argv[i] << endl;
383
384 inputFile = TFile::Open(argv[i]);
385
386 if(inputFile == NULL)
387 {
388 message << "can't open " << argv[i] << endl;
389 throw runtime_error(message.str());
390 }
391
392 fwlite::Event event(inputFile);
393
394 numberOfEvents = event.size();
395
396 if(numberOfEvents <= 0) continue;
397
398 // ExRootProgressBar progressBar(numberOfEvents - 1);
399 ExRootProgressBar progressBar(-1);
400
401 // Loop over all objects
402 eventCounter = 0;
403 modularDelphes->Clear();
404 treeWriter->Clear();
405
406 for(event.toBegin(); !event.atEnd() && !interrupted && (maxEvents <= 0 || eventCounter-skipEvents < maxEvents); ++event)
407 {
408 if(eventCounter >= skipEvents){
409 ConvertInput(event, eventCounter, branchEvent, branchWeight, factory,
410 allParticleOutputArray, stableParticleOutputArray, partonOutputArray, firstEvent);
411 modularDelphes->ProcessTask();
412
413 firstEvent = kFALSE;
414
415 treeWriter->Fill();
416
417 modularDelphes->Clear();
418 treeWriter->Clear();
419 }
420
421 progressBar.Update(eventCounter, eventCounter);
422 ++eventCounter;
423 }
424
425 progressBar.Update(eventCounter, eventCounter, kTRUE);
426 progressBar.Finish();
427
428 inputFile->Close();
429 }
430
431 modularDelphes->FinishTask();
432 treeWriter->Write();
433
434 cout << "** Exiting..." << endl;
435
436 delete modularDelphes;
437 delete confReader;
438 delete treeWriter;
439 delete outputFile;
440
441 return 0;
442 }
443 catch(runtime_error &e)
444 {
445 if(treeWriter) delete treeWriter;
446 if(outputFile) delete outputFile;
447 cerr << "** ERROR: " << e.what() << endl;
448 return 1;
449 }
450}
Note: See TracBrowser for help on using the repository browser.