Fork me on GitHub

source: git/readers/DelphesCMSFWLite.cpp@ 91acc76

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 91acc76 was 91acc76, checked in by Michele Selvaggi <michele.selvaggi@…>, 8 years ago

resolved conflict in ModulesLinkDef

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