Fork me on GitHub

source: svn/trunk/readers/DelphesProMC.cpp@ 1183

Last change on this file since 1183 was 1182, checked in by Pavel Demin, 11 years ago

add stop watches

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