Fork me on GitHub

source: svn/trunk/readers/DelphesCMSFWLite.cpp

Last change on this file was 1395, checked in by Pavel Demin, 10 years ago

synchronize with GitHub

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