Fork me on GitHub

source: git/converters/lhco2root.cpp@ 622e3e8

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 622e3e8 was 1fa50c2, checked in by Pavel Demin <pavel.demin@…>, 10 years ago

fix GPLv3 header

  • Property mode set to 100644
File size: 12.7 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 <stdexcept>
20#include <iostream>
21#include <sstream>
22
23#include <stdlib.h>
24#include <signal.h>
25#include <stdio.h>
26
27#include "TROOT.h"
28#include "TApplication.h"
29
30#include "TFile.h"
31#include "TObjArray.h"
32#include "TStopwatch.h"
33#include "TDatabasePDG.h"
34#include "TParticlePDG.h"
35#include "TLorentzVector.h"
36
37#include "modules/Delphes.h"
38#include "classes/DelphesStream.h"
39#include "classes/DelphesClasses.h"
40#include "classes/DelphesFactory.h"
41
42#include "ExRootAnalysis/ExRootTreeWriter.h"
43#include "ExRootAnalysis/ExRootTreeBranch.h"
44#include "ExRootAnalysis/ExRootProgressBar.h"
45
46using namespace std;
47
48static const int kBufferSize = 1024;
49
50/*
51LHC Olympics format discription from http://www.jthaler.net/olympicswiki/doku.php?id=lhc_olympics:data_file_format
52
53 * The first column of each row is just a counter that labels the object.
54 * The event begins with a row labelled "0"; this row contains the event number and the triggering information. The last row of the event is always the missing transverse momentum (MET).
55 * The second column of each row gives the type of object being listed [0, 1, 2, 3, 4, 6 = photon, electron, muon, hadronically-decaying tau, jet, missing transverse energy].
56 * The next three columns give the pseudorapidity, the azimuthal angle, and the transverse momentum of the object.
57 * The sixth column gives the invariant mass of the object.
58 * The seventh column gives the number of tracks associated with the object; in the case of a lepton, this number is multiplied by the charge of the lepton.
59 * The eighth column is 1 or 2 for a jet that has been "tagged" as containing a b-quark (actually a heavy flavor tag that sometimes indicates c-quarks), otherwise it is 0. For muons, the integer part of this number is the identity of the jet (see column 1) that is closest ot this muon in Delta R.
60 * The ninth column is the ratio of the hadronic versus electromagnetic energy deposited in the calorimeter cells associated with the object. For muons to the left of the decimal point is the summed pT in a R=0.4 cone (excluding the muon). To the right of the decimal point is etrat, which is a percentage between .00 and .99. It is the ratio of the transverse energy in a 3x3 grid surrounding the muon to the pT of the muon.
61*/
62
63//------------------------------------------------------------------------------
64
65class LHCOConverter
66{
67public:
68 LHCOConverter(TFile *outputFile);
69 ~LHCOConverter();
70
71 void Write();
72
73 Bool_t ReadLine(FILE *inputFile);
74
75private:
76
77 void AddMissingEvents();
78
79 void AnalyseEvent(ExRootTreeBranch *branch);
80
81 void AnalysePhoton(ExRootTreeBranch *branch);
82 void AnalyseElectron(ExRootTreeBranch *branch);
83 void AnalyseMuon(ExRootTreeBranch *branch);
84 void AnalyseTau(ExRootTreeBranch *branch);
85 void AnalyseJet(ExRootTreeBranch *branch);
86 void AnalyseMissingET(ExRootTreeBranch *branch);
87
88 enum {kIntParamSize = 2, kDblParamSize = 7};
89 Int_t fIntParam[kIntParamSize];
90 Double_t fDblParam[kDblParamSize];
91
92 Bool_t fIsReadyToFill;
93
94 Int_t fTriggerWord, fEventNumber;
95
96 char *fBuffer;
97
98 ExRootTreeWriter *fTreeWriter;
99
100 ExRootTreeBranch *fBranchEvent;
101 ExRootTreeBranch *fBranchTrack;
102 ExRootTreeBranch *fBranchTower;
103 ExRootTreeBranch *fBranchPhoton;
104 ExRootTreeBranch *fBranchElectron;
105 ExRootTreeBranch *fBranchMuon;
106 ExRootTreeBranch *fBranchJet;
107 ExRootTreeBranch *fBranchMissingET;
108
109};
110
111//------------------------------------------------------------------------------
112
113LHCOConverter::LHCOConverter(TFile *outputFile) :
114 fIsReadyToFill(kFALSE),
115 fTriggerWord(0), fEventNumber(1),
116 fBuffer(0), fTreeWriter(0)
117{
118 fBuffer = new char[kBufferSize];
119 fTreeWriter = new ExRootTreeWriter(outputFile, "Delphes");
120
121 // information about reconstructed event
122 fBranchEvent = fTreeWriter->NewBranch("Event", LHCOEvent::Class());
123 // reconstructed tracks
124 fBranchTrack = fTreeWriter->NewBranch("Track", Track::Class());
125 /// calorimeter towers
126 fBranchTower = fTreeWriter->NewBranch("Tower", Tower::Class());
127 // reconstructed photons
128 fBranchPhoton = fTreeWriter->NewBranch("Photon", Photon::Class());
129 // reconstructed electrons
130 fBranchElectron = fTreeWriter->NewBranch("Electron", Electron::Class());
131 // reconstructed muons
132 fBranchMuon = fTreeWriter->NewBranch("Muon", Muon::Class());
133 // reconstructed jets
134 fBranchJet = fTreeWriter->NewBranch("Jet", Jet::Class());
135 // missing transverse energy
136 fBranchMissingET = fTreeWriter->NewBranch("MissingET", MissingET::Class());
137}
138
139//------------------------------------------------------------------------------
140
141LHCOConverter::~LHCOConverter()
142{
143 if(fTreeWriter) delete fTreeWriter;
144 if(fBuffer) delete[] fBuffer;
145}
146
147//------------------------------------------------------------------------------
148
149Bool_t LHCOConverter::ReadLine(FILE *inputFile)
150{
151 int rc;
152
153 if(!fgets(fBuffer, kBufferSize, inputFile)) return kFALSE;
154
155 DelphesStream bufferStream(fBuffer);
156
157 rc = bufferStream.ReadInt(fIntParam[0]);
158
159 if(!rc)
160 {
161 return kTRUE;
162 }
163
164 if(fIntParam[0] == 0)
165 {
166 rc = bufferStream.ReadInt(fEventNumber)
167 && bufferStream.ReadInt(fTriggerWord);
168
169 if(!rc)
170 {
171 cerr << "** ERROR: " << "invalid event format" << endl;
172 return kFALSE;
173 }
174
175 if(fIsReadyToFill && fTreeWriter)
176 {
177 fTreeWriter->Fill();
178 fTreeWriter->Clear();
179 }
180
181 AnalyseEvent(fBranchEvent);
182 fIsReadyToFill = kTRUE;
183 }
184 else
185 {
186 rc = bufferStream.ReadInt(fIntParam[1])
187 && bufferStream.ReadDbl(fDblParam[0])
188 && bufferStream.ReadDbl(fDblParam[1])
189 && bufferStream.ReadDbl(fDblParam[2])
190 && bufferStream.ReadDbl(fDblParam[3])
191 && bufferStream.ReadDbl(fDblParam[4])
192 && bufferStream.ReadDbl(fDblParam[5])
193 && bufferStream.ReadDbl(fDblParam[6]);
194
195 if(!rc)
196 {
197 cerr << "** ERROR: " << "invalid object format" << endl;
198 return kFALSE;
199 }
200
201 switch(fIntParam[1])
202 {
203 case 0: AnalysePhoton(fBranchPhoton); break;
204 case 1: AnalyseElectron(fBranchElectron); break;
205 case 2: AnalyseMuon(fBranchMuon); break;
206 case 3: AnalyseTau(fBranchJet); break;
207 case 4: AnalyseJet(fBranchJet); break;
208 case 6: AnalyseMissingET(fBranchMissingET); break;
209 }
210 }
211
212 return kTRUE;
213}
214
215//---------------------------------------------------------------------------
216
217void LHCOConverter::Write()
218{
219 if(fIsReadyToFill && fTreeWriter) fTreeWriter->Fill();
220 if(fTreeWriter) fTreeWriter->Write();
221 fIsReadyToFill = kFALSE;
222}
223
224//---------------------------------------------------------------------------
225
226void LHCOConverter::AnalyseEvent(ExRootTreeBranch *branch)
227{
228 LHCOEvent *element;
229
230 element = static_cast<LHCOEvent*>(branch->NewEntry());
231
232 element->Number = fEventNumber;
233 element->Trigger = fTriggerWord;
234}
235
236//---------------------------------------------------------------------------
237
238void LHCOConverter::AnalysePhoton(ExRootTreeBranch *branch)
239{
240 Photon *element;
241
242 element = static_cast<Photon*>(branch->NewEntry());
243
244 element->Eta = fDblParam[0];
245 element->Phi = fDblParam[1];
246 element->PT = fDblParam[2];
247 element->EhadOverEem = fDblParam[6];
248}
249
250//---------------------------------------------------------------------------
251
252void LHCOConverter::AnalyseElectron(ExRootTreeBranch *branch)
253{
254 Electron *element;
255
256 element = static_cast<Electron*>(branch->NewEntry());
257
258 element->Eta = fDblParam[0];
259 element->Phi = fDblParam[1];
260 element->PT = fDblParam[2];
261
262 element->Charge = fDblParam[4] < 0.0 ? -1 : 1;
263/*
264 element->Ntrk = TMath::Abs(fDblParam[4]);
265*/
266 element->EhadOverEem = fDblParam[6];
267}
268
269//---------------------------------------------------------------------------
270
271void LHCOConverter::AnalyseMuon(ExRootTreeBranch *branch)
272{
273 Muon *element;
274
275 element = static_cast<Muon*>(branch->NewEntry());
276
277 element->Eta = fDblParam[0];
278 element->Phi = fDblParam[1];
279 element->PT = fDblParam[2];
280
281 element->Charge = fDblParam[4] < 0.0 ? -1 : 1;
282/*
283 element->Ntrk = TMath::Abs(fDblParam[4]);
284
285 element->JetIndex = Int_t(fDblParam[5]);
286
287 element->PTiso = Int_t(fDblParam[6]);
288 element->ETiso = fDblParam[6] - element->PTiso;
289*/
290}
291
292//---------------------------------------------------------------------------
293
294void LHCOConverter::AnalyseTau(ExRootTreeBranch *branch)
295{
296 Jet *element;
297
298 element = static_cast<Jet*>(branch->NewEntry());
299
300 element->Eta = fDblParam[0];
301 element->Phi = fDblParam[1];
302 element->PT = fDblParam[2];
303
304 element->Mass = fDblParam[3];
305
306 element->BTag = 0;
307 element->TauTag = 1;
308
309 element->Charge = fDblParam[4] < 0 ? -1 : 1;
310/*
311 element->Ntrk = TMath::Abs(fDblParam[4]);
312*/
313 element->EhadOverEem = fDblParam[6];
314}
315
316//---------------------------------------------------------------------------
317
318void LHCOConverter::AnalyseJet(ExRootTreeBranch *branch)
319{
320 Jet *element;
321
322 element = static_cast<Jet*>(branch->NewEntry());
323
324 element->Eta = fDblParam[0];
325 element->Phi = fDblParam[1];
326 element->PT = fDblParam[2];
327
328 element->Mass = fDblParam[3];
329/*
330 element->Ntrk = TMath::Abs(Int_t(fDblParam[4]));
331*/
332 element->BTag = Int_t(fDblParam[5]);
333 element->TauTag = 0;
334
335 element->Charge = 0;
336
337 element->EhadOverEem = fDblParam[6];
338/*
339 element->Index = fIntParam[0];
340*/
341}
342
343//---------------------------------------------------------------------------
344
345void LHCOConverter::AnalyseMissingET(ExRootTreeBranch *branch)
346{
347 MissingET *element;
348
349 element = static_cast<MissingET*>(branch->NewEntry());
350
351 element->Phi = fDblParam[1];
352 element->MET = fDblParam[2];
353}
354
355//---------------------------------------------------------------------------
356
357static bool interrupted = false;
358
359void SignalHandler(int sig)
360{
361 interrupted = true;
362}
363
364//---------------------------------------------------------------------------
365
366int main(int argc, char *argv[])
367{
368 char appName[] = "lhco2root";
369 stringstream message;
370 FILE *inputFile = 0;
371 TFile *outputFile = 0;
372 LHCOConverter *converter = 0;
373 Int_t i;
374 Long64_t length, eventCounter;
375
376 if(argc < 2)
377 {
378 cout << " Usage: " << appName << " output_file" << " [input_file(s)]" << endl;
379 cout << " output_file - output file in ROOT format," << endl;
380 cout << " input_file(s) - input file(s) in LHCO format," << endl;
381 cout << " with no input_file, or when input_file is -, read standard input." << endl;
382 return 1;
383 }
384
385 signal(SIGINT, SignalHandler);
386
387 gROOT->SetBatch();
388
389 int appargc = 1;
390 char *appargv[] = {appName};
391 TApplication app(appName, &appargc, appargv);
392
393 try
394 {
395 outputFile = TFile::Open(argv[1], "CREATE");
396
397 if(outputFile == NULL)
398 {
399 message << "can't open " << argv[1];
400 throw runtime_error(message.str());
401 }
402
403 converter = new LHCOConverter(outputFile);
404
405 i = 2;
406 do
407 {
408 if(interrupted) break;
409
410 if(i == argc || strcmp(argv[i], "-") == 0)
411 {
412 cout << "** Reading standard input" << endl;
413 inputFile = stdin;
414 length = -1;
415 }
416 else
417 {
418 cout << "** Reading " << argv[i] << endl;
419 inputFile = fopen(argv[i], "r");
420
421 if(inputFile == NULL)
422 {
423 message << "can't open " << argv[i];
424 throw runtime_error(message.str());
425 }
426
427 fseek(inputFile, 0L, SEEK_END);
428 length = ftello(inputFile);
429 fseek(inputFile, 0L, SEEK_SET);
430
431 if(length <= 0)
432 {
433 fclose(inputFile);
434 ++i;
435 continue;
436 }
437 }
438
439 eventCounter = 0;
440 ExRootProgressBar progressBar(length);
441
442 // Loop over all objects
443 while(converter->ReadLine(inputFile) && !interrupted)
444 {
445 ++eventCounter;
446
447 progressBar.Update(ftello(inputFile), eventCounter);
448 }
449 converter->Write();
450
451 fseek(inputFile, 0L, SEEK_END);
452 progressBar.Update(ftello(inputFile), eventCounter, kTRUE);
453 progressBar.Finish();
454
455 if(inputFile != stdin) fclose(inputFile);
456
457 ++i;
458 }
459 while(i < argc);
460
461 cout << "** Exiting..." << endl;
462
463 delete converter;
464 delete outputFile;
465
466 return 0;
467 }
468 catch(runtime_error &e)
469 {
470 if(converter) delete converter;
471 if(outputFile) delete outputFile;
472 cerr << "** ERROR: " << e.what() << endl;
473 return 1;
474 }
475}
476
477
Note: See TracBrowser for help on using the repository browser.