Fork me on GitHub

source: svn/trunk/readers/DelphesLHEF.cpp

Last change on this file was 1265, checked in by Pavel Demin, 11 years ago

read reweighting information from LHEF

File size: 5.6 KB
Line 
1#include <stdexcept>
2#include <iostream>
3#include <sstream>
4
5#include <signal.h>
6
7#include "TROOT.h"
8#include "TApplication.h"
9
10#include "TFile.h"
11#include "TObjArray.h"
12#include "TStopwatch.h"
13#include "TDatabasePDG.h"
14#include "TParticlePDG.h"
15#include "TLorentzVector.h"
16
17#include "modules/Delphes.h"
18#include "classes/DelphesClasses.h"
19#include "classes/DelphesFactory.h"
20#include "classes/DelphesLHEFReader.h"
21
22#include "ExRootAnalysis/ExRootTreeWriter.h"
23#include "ExRootAnalysis/ExRootTreeBranch.h"
24#include "ExRootAnalysis/ExRootProgressBar.h"
25
26using namespace std;
27
28//---------------------------------------------------------------------------
29
30static bool interrupted = false;
31
32void SignalHandler(int sig)
33{
34 interrupted = true;
35}
36
37//---------------------------------------------------------------------------
38
39int main(int argc, char *argv[])
40{
41 char appName[] = "DelphesLHEF";
42 stringstream message;
43 FILE *inputFile = 0;
44 TFile *outputFile = 0;
45 TStopwatch readStopWatch, procStopWatch;
46 ExRootTreeWriter *treeWriter = 0;
47 ExRootTreeBranch *branchEvent = 0, *branchRwgt = 0;
48 ExRootConfReader *confReader = 0;
49 Delphes *modularDelphes = 0;
50 DelphesFactory *factory = 0;
51 TObjArray *stableParticleOutputArray = 0, *allParticleOutputArray = 0, *partonOutputArray = 0;
52 DelphesLHEFReader *reader = 0;
53 Int_t i, maxEvents, skipEvents;
54 Long64_t length, eventCounter;
55
56 if(argc < 3)
57 {
58 cout << " Usage: " << appName << " config_file" << " output_file" << " [input_file(s)]" << endl;
59 cout << " config_file - configuration file in Tcl format," << endl;
60 cout << " output_file - output file in ROOT format," << endl;
61 cout << " input_file(s) - input file(s) in LHEF format," << endl;
62 cout << " with no input_file, or when input_file is -, read standard input." << endl;
63 return 1;
64 }
65
66 signal(SIGINT, SignalHandler);
67
68 gROOT->SetBatch();
69
70 int appargc = 1;
71 char *appargv[] = {appName};
72 TApplication app(appName, &appargc, appargv);
73
74 try
75 {
76 outputFile = TFile::Open(argv[2], "CREATE");
77
78 if(outputFile == NULL)
79 {
80 message << "can't create output file " << argv[2];
81 throw runtime_error(message.str());
82 }
83
84 treeWriter = new ExRootTreeWriter(outputFile, "Delphes");
85
86 branchEvent = treeWriter->NewBranch("Event", LHEFEvent::Class());
87 branchRwgt = treeWriter->NewBranch("Rwgt", Weight::Class());
88
89 confReader = new ExRootConfReader;
90 confReader->ReadFile(argv[1]);
91
92 maxEvents = confReader->GetInt("::MaxEvents", 0);
93 skipEvents = confReader->GetInt("::SkipEvents", 0);
94
95 if(maxEvents < 0)
96 {
97 throw runtime_error("MaxEvents must be zero or positive");
98 }
99
100 if(skipEvents < 0)
101 {
102 throw runtime_error("SkipEvents must be zero or positive");
103 }
104
105 modularDelphes = new Delphes("Delphes");
106 modularDelphes->SetConfReader(confReader);
107 modularDelphes->SetTreeWriter(treeWriter);
108
109 factory = modularDelphes->GetFactory();
110 allParticleOutputArray = modularDelphes->ExportArray("allParticles");
111 stableParticleOutputArray = modularDelphes->ExportArray("stableParticles");
112 partonOutputArray = modularDelphes->ExportArray("partons");
113
114 reader = new DelphesLHEFReader;
115
116 modularDelphes->InitTask();
117
118 i = 3;
119 do
120 {
121 if(interrupted) break;
122
123 if(i == argc || strncmp(argv[i], "-", 2) == 0)
124 {
125 cout << "** Reading standard input" << endl;
126 inputFile = stdin;
127 length = -1;
128 }
129 else
130 {
131 cout << "** Reading " << argv[i] << endl;
132 inputFile = fopen(argv[i], "r");
133
134 if(inputFile == NULL)
135 {
136 message << "can't open " << argv[i];
137 throw runtime_error(message.str());
138 }
139
140 fseek(inputFile, 0L, SEEK_END);
141 length = ftello(inputFile);
142 fseek(inputFile, 0L, SEEK_SET);
143
144 if(length <= 0)
145 {
146 fclose(inputFile);
147 ++i;
148 continue;
149 }
150 }
151
152 reader->SetInputFile(inputFile);
153
154 ExRootProgressBar progressBar(length);
155
156 // Loop over all objects
157 eventCounter = 0;
158 treeWriter->Clear();
159 modularDelphes->Clear();
160 reader->Clear();
161 readStopWatch.Start();
162 while((maxEvents <= 0 || eventCounter - skipEvents < maxEvents) &&
163 reader->ReadBlock(factory, allParticleOutputArray,
164 stableParticleOutputArray, partonOutputArray) && !interrupted)
165 {
166 if(reader->EventReady())
167 {
168 ++eventCounter;
169
170 readStopWatch.Stop();
171
172 if(eventCounter > skipEvents)
173 {
174 readStopWatch.Stop();
175 procStopWatch.Start();
176 modularDelphes->ProcessTask();
177 procStopWatch.Stop();
178
179 reader->AnalyzeEvent(branchEvent, eventCounter, &readStopWatch, &procStopWatch);
180 reader->AnalyzeRwgt(branchRwgt);
181
182 treeWriter->Fill();
183
184 treeWriter->Clear();
185 }
186
187 modularDelphes->Clear();
188 reader->Clear();
189
190 readStopWatch.Start();
191 }
192 progressBar.Update(ftello(inputFile), eventCounter);
193 }
194
195 fseek(inputFile, 0L, SEEK_END);
196 progressBar.Update(ftello(inputFile), eventCounter, kTRUE);
197 progressBar.Finish();
198
199 if(inputFile != stdin) fclose(inputFile);
200
201 ++i;
202 }
203 while(i < argc);
204
205 modularDelphes->FinishTask();
206 treeWriter->Write();
207
208 cout << "** Exiting..." << endl;
209
210 delete reader;
211 delete modularDelphes;
212 delete confReader;
213 delete treeWriter;
214 delete outputFile;
215
216 return 0;
217 }
218 catch(runtime_error &e)
219 {
220 if(treeWriter) delete treeWriter;
221 if(outputFile) delete outputFile;
222 cerr << "** ERROR: " << e.what() << endl;
223 return 1;
224 }
225}
Note: See TracBrowser for help on using the repository browser.