Fork me on GitHub

source: git/readers/DelphesCMSFWLite.cpp@ ebe36ed

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

correct LHEEventProduct

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