Fork me on GitHub

source: git/readers/DelphesROOT.cpp@ 3d19203

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

remove broken MaxEvent param in DelphesROOT

  • Property mode set to 100644
File size: 7.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 "TClonesArray.h"
37#include "TObjArray.h"
38#include "TStopwatch.h"
39#include "TDatabasePDG.h"
40#include "TParticlePDG.h"
41#include "TLorentzVector.h"
42
43#include "modules/Delphes.h"
44#include "classes/DelphesStream.h"
45#include "classes/DelphesClasses.h"
46#include "classes/DelphesFactory.h"
47
48#include "ExRootAnalysis/ExRootTreeWriter.h"
49#include "ExRootAnalysis/ExRootTreeReader.h"
50#include "ExRootAnalysis/ExRootTreeBranch.h"
51#include "ExRootAnalysis/ExRootProgressBar.h"
52
53
54
55using namespace std;
56
57//---------------------------------------------------------------------------
58
59
60//---------------------------------------------------------------------------
61
62static bool interrupted = false;
63
64void SignalHandler(int sig)
65{
66 interrupted = true;
67}
68
69//---------------------------------------------------------------------------
70
71int main(int argc, char *argv[])
72{
73 char appName[] = "DelphesROOT";
74 stringstream message;
75 TFile *inputFile = 0;
76 TFile *outputFile = 0;
77 TStopwatch eventStopWatch;
78 ExRootTreeWriter *treeWriter = 0;
79 ExRootTreeBranch *branchEvent = 0;
80 ExRootConfReader *confReader = 0;
81 Delphes *modularDelphes = 0;
82 DelphesFactory *factory = 0;
83 GenParticle *gen;
84 HepMCEvent *element, *eve;
85 Candidate *candidate;
86 Int_t pdgCode;
87 Int_t maxEvents;
88
89 const Double_t c_light = 2.99792458E8;
90
91 TObjArray *allParticleOutputArray = 0, *stableParticleOutputArray = 0, *partonOutputArray = 0;
92 Int_t i;
93 Long64_t eventCounter, numberOfEvents;
94
95 if(argc < 4)
96 {
97 cout << " Usage: " << appName << " config_file" << " output_file" << " input_file(s)" << endl;
98 cout << " config_file - configuration file in Tcl format," << endl;
99 cout << " output_file - output file in ROOT format," << endl;
100 cout << " input_file(s) - input file(s) in ROOT format." << endl;
101 return 1;
102 }
103
104 signal(SIGINT, SignalHandler);
105
106 gROOT->SetBatch();
107
108 int appargc = 1;
109 char *appargv[] = {appName};
110 TApplication app(appName, &appargc, appargv);
111
112 try
113 {
114 outputFile = TFile::Open(argv[2], "CREATE");
115
116 if(outputFile == NULL)
117 {
118 message << "can't open " << argv[2] << endl;
119 throw runtime_error(message.str());
120 }
121
122 treeWriter = new ExRootTreeWriter(outputFile, "Delphes");
123
124 branchEvent = treeWriter->NewBranch("Event", HepMCEvent::Class());
125
126 confReader = new ExRootConfReader;
127 confReader->ReadFile(argv[1]);
128
129 if(maxEvents < 0)
130 {
131 throw runtime_error("MaxEvents must be zero or positive");
132 }
133
134 modularDelphes = new Delphes("Delphes");
135 modularDelphes->SetConfReader(confReader);
136 modularDelphes->SetTreeWriter(treeWriter);
137
138 TChain *chain = new TChain("Delphes");
139
140 factory = modularDelphes->GetFactory();
141 allParticleOutputArray = modularDelphes->ExportArray("allParticles");
142 stableParticleOutputArray = modularDelphes->ExportArray("stableParticles");
143 partonOutputArray = modularDelphes->ExportArray("partons");
144
145 modularDelphes->InitTask();
146
147 for(i = 3; i < argc && !interrupted; ++i)
148 {
149 cout << "** Reading " << argv[i] << endl;
150
151 chain->Add(argv[i]);
152 ExRootTreeReader *treeReader = new ExRootTreeReader(chain);
153
154 inputFile = TFile::Open(argv[i]);
155
156 if(inputFile == NULL)
157 {
158 message << "can't open " << argv[i] << endl;
159 throw runtime_error(message.str());
160 }
161
162 numberOfEvents = treeReader->GetEntries();
163 TClonesArray *branchParticle = treeReader->UseBranch("Particle");
164 TClonesArray *branchHepMCEvent = treeReader->UseBranch("Event");
165
166 if(numberOfEvents <= 0) continue;
167
168 // ExRootProgressBar progressBar(numberOfEvents - 1);
169 ExRootProgressBar progressBar(-1);
170
171 // Loop over all objects
172 eventCounter = 0;
173 modularDelphes->Clear();
174 treeWriter->Clear();
175 for(Int_t entry = 0; entry < numberOfEvents && !interrupted; ++entry)
176 {
177
178 treeReader->ReadEntry(entry);
179
180 // -- TBC need also to include event weights --
181
182 eve = (HepMCEvent*) branchHepMCEvent->At(0);
183 element = static_cast<HepMCEvent *>(branchEvent->NewEntry());
184
185 element->Number = eventCounter;
186
187 element->ProcessID = eve->ProcessID;
188 element->MPI = eve->MPI;
189 element->Weight = eve->Weight;
190 element->Scale = eve->Scale;
191 element->AlphaQED = eve->AlphaQED;
192 element->AlphaQCD = eve->AlphaQCD;
193
194 element->ID1 = eve->ID1;
195 element->ID2 = eve->ID2;
196 element->X1 = eve->X1;
197 element->X2 = eve->X2;
198 element->ScalePDF = eve->ScalePDF;
199 element->PDF1 = eve->PDF1;
200 element->PDF2 = eve->PDF2;
201
202 element->ReadTime = eve->ReadTime;
203 element->ProcTime = eve->ProcTime;
204
205 for(Int_t j=0; j < branchParticle->GetEntriesFast(); j++)
206 {
207
208 gen = (GenParticle*) branchParticle->At(j);
209 candidate = factory->NewCandidate();
210
211 candidate->Momentum = gen->P4();
212 candidate->Position.SetXYZT(gen->X, gen->Y, gen->Z, gen->T*1.0E3*c_light);
213
214 candidate->PID = gen->PID;
215 candidate->Status = gen->Status;
216
217 candidate->M1 = gen->M1;
218 candidate->M2 = gen->M2;
219
220 candidate->D1 = gen->D1;
221 candidate->D2 = gen->D2;
222
223 candidate->Charge = gen->Charge;
224 candidate->Mass = gen->Mass;
225
226 allParticleOutputArray->Add(candidate);
227
228 pdgCode = TMath::Abs(gen->PID);
229
230 if(gen->Status == 1)
231 {
232 stableParticleOutputArray->Add(candidate);
233 }
234 else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
235 {
236 partonOutputArray->Add(candidate);
237 }
238 }
239
240 modularDelphes->ProcessTask();
241
242 treeWriter->Fill();
243
244 modularDelphes->Clear();
245 treeWriter->Clear();
246
247 progressBar.Update(eventCounter, eventCounter);
248 ++eventCounter;
249 }
250
251 progressBar.Update(eventCounter, eventCounter, kTRUE);
252 progressBar.Finish();
253
254 inputFile->Close();
255
256 delete treeReader;
257
258 }
259
260 modularDelphes->FinishTask();
261 treeWriter->Write();
262
263 cout << "** Exiting..." << endl;
264
265 delete modularDelphes;
266 delete confReader;
267 delete treeWriter;
268 delete outputFile;
269 delete chain;
270
271 return 0;
272 }
273 catch(runtime_error &e)
274 {
275 if(treeWriter) delete treeWriter;
276 if(outputFile) delete outputFile;
277 cerr << "** ERROR: " << e.what() << endl;
278 return 1;
279 }
280}
Note: See TracBrowser for help on using the repository browser.