Fork me on GitHub

source: git/readers/DelphesCMSFWLite.cpp@ 97510e6

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 97510e6 was bba4646, checked in by pavel <pavel@…>, 11 years ago

use mass from input file

  • Property mode set to 100644
File size: 6.7 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
39using namespace std;
40
41//---------------------------------------------------------------------------
42
43void ConvertInput(fwlite::Event &event, DelphesFactory *factory, TObjArray *allParticleOutputArray, TObjArray *stableParticleOutputArray, TObjArray *partonOutputArray)
44{
45 fwlite::Handle< vector< reco::GenParticle > > handleParticle;
46 vector< reco::GenParticle >::const_iterator itParticle;
47
48 vector< const reco::Candidate * > vectorCandidate;
49 vector< const reco::Candidate * >::iterator itCandidate;
50
51 handleParticle.getByLabel(event, "genParticles");
52
53 Candidate *candidate;
54 TDatabasePDG *pdg;
55 TParticlePDG *pdgParticle;
56 Int_t pdgCode;
57
58 Int_t pid, status;
59 Double_t px, py, pz, e, mass;
60 Double_t x, y, z;
61
62 pdg = TDatabasePDG::Instance();
63
64 for(itParticle = handleParticle->begin(); itParticle != handleParticle->end(); ++itParticle)
65 {
66 vectorCandidate.push_back(&*itParticle);
67 }
68
69 for(itParticle = handleParticle->begin(); itParticle != handleParticle->end(); ++itParticle)
70 {
71 const reco::GenParticle &particle = *itParticle;
72
73 pid = particle.pdgId();
74 status = particle.status();
75 px = particle.px(); py = particle.py(); pz = particle.pz(); e = particle.energy(); mass = particle.mass();
76 x = particle.vx(); y = particle.vy(); z = particle.vz();
77
78 candidate = factory->NewCandidate();
79
80 candidate->PID = pid;
81 pdgCode = TMath::Abs(candidate->PID);
82
83 candidate->Status = status;
84
85 itCandidate = find(vectorCandidate.begin(), vectorCandidate.end(), particle.daughter(0));
86 if(itCandidate != vectorCandidate.end()) candidate->D1 = distance(vectorCandidate.begin(), itCandidate);
87
88 itCandidate = find(vectorCandidate.begin(), vectorCandidate.end(), particle.daughter(particle.numberOfDaughters() - 1));
89 if(itCandidate != vectorCandidate.end()) candidate->D2 = distance(vectorCandidate.begin(), itCandidate);
90
91 pdgParticle = pdg->GetParticle(pid);
92 candidate->Charge = pdgParticle ? Int_t(pdgParticle->Charge()/3.0) : -999;
93 candidate->Mass = mass;
94
95 candidate->Momentum.SetPxPyPzE(px, py, pz, e);
96
97 candidate->Position.SetXYZT(x, y, z, 0.0);
98
99 allParticleOutputArray->Add(candidate);
100
101 if(!pdgParticle) continue;
102
103 if(status == 1)
104 {
105 stableParticleOutputArray->Add(candidate);
106 }
107 else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
108 {
109 partonOutputArray->Add(candidate);
110 }
111 }
112}
113
114//---------------------------------------------------------------------------
115
116static bool interrupted = false;
117
118void SignalHandler(int sig)
119{
120 interrupted = true;
121}
122
123//---------------------------------------------------------------------------
124
125int main(int argc, char *argv[])
126{
127 char appName[] = "DelphesCMSFWLite";
128 stringstream message;
129 TFile *inputFile = 0;
130 TFile *outputFile = 0;
131 TStopwatch eventStopWatch;
132 ExRootTreeWriter *treeWriter = 0;
133 ExRootConfReader *confReader = 0;
134 Delphes *modularDelphes = 0;
135 DelphesFactory *factory = 0;
136 TObjArray *allParticleOutputArray = 0, *stableParticleOutputArray = 0, *partonOutputArray = 0;
137 Int_t i;
138 Long64_t eventCounter, numberOfEvents;
139
140 if(argc < 4)
141 {
142 cout << " Usage: " << appName << " config_file" << " output_file" << " input_file(s)" << endl;
143 cout << " config_file - configuration file in Tcl format," << endl;
144 cout << " output_file - output file in ROOT format," << endl;
145 cout << " input_file(s) - input file(s) in ROOT format." << endl;
146 return 1;
147 }
148
149 signal(SIGINT, SignalHandler);
150
151 gROOT->SetBatch();
152
153 int appargc = 1;
154 char *appargv[] = {appName};
155 TApplication app(appName, &appargc, appargv);
156
157 AutoLibraryLoader::enable();
158
159 try
160 {
161 outputFile = TFile::Open(argv[2], "CREATE");
162
163 if(outputFile == NULL)
164 {
165 message << "can't open " << argv[2] << endl;
166 throw runtime_error(message.str());
167 }
168
169 treeWriter = new ExRootTreeWriter(outputFile, "Delphes");
170
171 confReader = new ExRootConfReader;
172 confReader->ReadFile(argv[1]);
173
174 modularDelphes = new Delphes("Delphes");
175 modularDelphes->SetConfReader(confReader);
176 modularDelphes->SetTreeWriter(treeWriter);
177
178 factory = modularDelphes->GetFactory();
179 allParticleOutputArray = modularDelphes->ExportArray("allParticles");
180 stableParticleOutputArray = modularDelphes->ExportArray("stableParticles");
181 partonOutputArray = modularDelphes->ExportArray("partons");
182
183 modularDelphes->InitTask();
184
185 for(i = 3; i < argc && !interrupted; ++i)
186 {
187 cout << "** Reading " << argv[i] << endl;
188
189 inputFile = TFile::Open(argv[i]);
190
191 if(inputFile == NULL)
192 {
193 message << "can't open " << argv[i] << endl;
194 throw runtime_error(message.str());
195 }
196
197 fwlite::Event event(inputFile);
198
199 numberOfEvents = event.size();
200
201 if(numberOfEvents <= 0) continue;
202
203 // ExRootProgressBar progressBar(numberOfEvents - 1);
204 ExRootProgressBar progressBar(-1);
205
206 // Loop over all objects
207 eventCounter = 0;
208 modularDelphes->Clear();
209 treeWriter->Clear();
210 for(event.toBegin(); !event.atEnd() && !interrupted; ++event)
211 {
212 ConvertInput(event, factory, allParticleOutputArray, stableParticleOutputArray, partonOutputArray);
213 modularDelphes->ProcessTask();
214
215 treeWriter->Fill();
216
217 modularDelphes->Clear();
218 treeWriter->Clear();
219
220 progressBar.Update(eventCounter, eventCounter);
221 ++eventCounter;
222 }
223
224 progressBar.Update(eventCounter, eventCounter, kTRUE);
225 progressBar.Finish();
226
227 inputFile->Close();
228 }
229
230 modularDelphes->FinishTask();
231 treeWriter->Write();
232
233 cout << "** Exiting..." << endl;
234
235 delete modularDelphes;
236 delete confReader;
237 delete treeWriter;
238 delete outputFile;
239
240 return 0;
241 }
242 catch(runtime_error &e)
243 {
244 if(treeWriter) delete treeWriter;
245 if(outputFile) delete outputFile;
246 cerr << "** ERROR: " << e.what() << endl;
247 return 1;
248 }
249}
Note: See TracBrowser for help on using the repository browser.