Fork me on GitHub

source: git/converters/root2lhco.cpp@ ff37d75

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

check if ROOT file contains all required branches (close #340)

  • Property mode set to 100644
File size: 13.3 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 <fstream>
22#include <sstream>
23#include <string>
24
25#include <stdlib.h>
26#include <signal.h>
27#include <stdio.h>
28
29#include "TROOT.h"
30#include "TApplication.h"
31
32#include "TFile.h"
33#include "TClonesArray.h"
34
35#include "classes/DelphesClasses.h"
36
37#include "ExRootAnalysis/ExRootTreeReader.h"
38#include "ExRootAnalysis/ExRootProgressBar.h"
39
40using namespace std;
41
42/*
43LHC Olympics format discription from http://www.jthaler.net/olympicswiki/doku.php?id=lhc_olympics:data_file_format
44
45 * The first column of each row is just a counter that labels the object.
46 * 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).
47 * 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].
48 * The next three columns give the pseudorapidity, the azimuthal angle, and the transverse momentum of the object.
49 * The sixth column gives the invariant mass of the object.
50 * 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.
51 * 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.
52 * 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.
53*/
54
55class LHCOWriter
56{
57public:
58 LHCOWriter(ExRootTreeReader *treeReader, FILE *outputFile);
59 ~LHCOWriter();
60
61 void ProcessEvent();
62
63private:
64
65 void Reset();
66 void Write();
67
68 void AnalyseEvent();
69
70 void AnalysePhotons();
71 void AnalyseElectrons();
72 void AnalyseMuons();
73 void AnalyseTauJets();
74 void AnalyseJets();
75
76 void AnalyseMissingET();
77
78 enum {kIntParamSize = 2, kDblParamSize = 9};
79 Int_t fIntParam[kIntParamSize];
80 Double_t fDblParam[kDblParamSize];
81
82 Long64_t fTriggerWord, fEventNumber;
83
84 ExRootTreeReader *fTreeReader;
85 FILE *fOutputFile;
86
87 TClonesArray *fBranchEvent;
88
89 TClonesArray *fBranchTrack;
90 TClonesArray *fBranchTower;
91
92 TClonesArray *fBranchPhoton;
93 TClonesArray *fBranchElectron;
94 TClonesArray *fBranchMuon;
95 TClonesArray *fBranchJet;
96 TClonesArray *fBranchMissingET;
97
98 TIterator *fItTrack;
99 TIterator *fItTower;
100
101 TIterator *fItPhoton;
102 TIterator *fItElectron;
103 TIterator *fItMuon;
104 TIterator *fItJet;
105};
106
107//------------------------------------------------------------------------------
108
109LHCOWriter::LHCOWriter(ExRootTreeReader *treeReader, FILE *outputFile) :
110 fTriggerWord(0), fEventNumber(1), fTreeReader(0), fOutputFile(0),
111 fBranchEvent(0), fBranchTrack(0), fBranchTower(0), fBranchPhoton(0),
112 fBranchElectron(0), fBranchMuon(0), fBranchJet(0), fBranchMissingET(0)
113{
114 fTreeReader = treeReader;
115 fOutputFile = outputFile;
116
117 // information about reconstructed event
118 fBranchEvent = fTreeReader->UseBranch("Event");
119 // reconstructed tracks
120 fBranchTrack = fTreeReader->UseBranch("Track");
121 // calorimeter towers
122 fBranchTower = fTreeReader->UseBranch("Tower");
123 // reconstructed photons
124 fBranchPhoton = fTreeReader->UseBranch("Photon");
125 // reconstructed electrons
126 fBranchElectron = fTreeReader->UseBranch("Electron");
127 // reconstructed muons
128 fBranchMuon = fTreeReader->UseBranch("Muon");
129 // reconstructed jets
130 fBranchJet = fTreeReader->UseBranch("Jet");
131 // missing transverse energy
132 fBranchMissingET = fTreeReader->UseBranch("MissingET");
133
134 if(!fBranchEvent || !fBranchTrack || !fBranchTower || !fBranchPhoton ||
135 !fBranchElectron || !fBranchMuon || !fBranchJet || !fBranchMissingET)
136 {
137 throw runtime_error("ROOT file doesn't contain all required branches");
138 }
139
140 fItTrack = fBranchTrack->MakeIterator();
141 fItTower = fBranchTower->MakeIterator();
142 fItPhoton = fBranchPhoton->MakeIterator();
143 fItElectron = fBranchElectron->MakeIterator();
144 fItMuon = fBranchMuon->MakeIterator();
145 fItJet = fBranchJet->MakeIterator();
146}
147
148//------------------------------------------------------------------------------
149
150LHCOWriter::~LHCOWriter()
151{
152}
153
154//---------------------------------------------------------------------------
155
156void LHCOWriter::ProcessEvent()
157{
158 fIntParam[0] = 0;
159
160 AnalyseEvent();
161
162 AnalysePhotons();
163 AnalyseElectrons();
164 AnalyseMuons();
165 AnalyseTauJets();
166 AnalyseJets();
167
168 AnalyseMissingET();
169}
170
171//---------------------------------------------------------------------------
172
173void LHCOWriter::Reset()
174{
175 int i;
176 for(i = 1; i < kIntParamSize; ++i)
177 {
178 fIntParam[i] = 0;
179 }
180
181 for(i = 0; i < kDblParamSize; ++i)
182 {
183 fDblParam[i] = 0.0;
184 }
185}
186
187//---------------------------------------------------------------------------
188
189void LHCOWriter::Write()
190{
191 fprintf(fOutputFile, "%4d %4d %8.3f %8.3f %7.2f %7.2f %6.1f %6.1f %7.2f %6.1f %6.1f\n",
192 fIntParam[0], fIntParam[1], fDblParam[0], fDblParam[1], fDblParam[2],
193 fDblParam[3], fDblParam[4], fDblParam[5], fDblParam[6], fDblParam[7], fDblParam[8]);
194
195 ++fIntParam[0];
196}
197
198//---------------------------------------------------------------------------
199
200void LHCOWriter::AnalyseEvent()
201{
202 Event *element;
203
204 element = static_cast<Event*>(fBranchEvent->At(0));
205
206 fprintf(fOutputFile, "%4d %13lld %8d\n", 0, element->Number, 0);
207
208 ++fIntParam[0];
209}
210
211//---------------------------------------------------------------------------
212
213void LHCOWriter::AnalysePhotons()
214{
215 Photon *element;
216
217 fItPhoton->Reset();
218 while((element = static_cast<Photon*>(fItPhoton->Next())))
219 {
220 Reset();
221
222 fIntParam[1] = 0;
223
224 fDblParam[0] = element->Eta;
225 fDblParam[1] = element->Phi;
226 fDblParam[2] = element->PT;
227
228 fDblParam[6] = element->EhadOverEem;
229
230 Write();
231 }
232}
233
234//---------------------------------------------------------------------------
235
236void LHCOWriter::AnalyseElectrons()
237{
238 Electron *element;
239
240 fItElectron->Reset();
241 while((element = static_cast<Electron*>(fItElectron->Next())))
242 {
243 Reset();
244
245 fIntParam[1] = 1;
246
247 fDblParam[0] = element->Eta;
248 fDblParam[1] = element->Phi;
249 fDblParam[2] = element->PT;
250
251 fDblParam[4] = element->Charge;
252
253 fDblParam[6] = element->EhadOverEem;
254
255 Write();
256 }
257}
258
259//---------------------------------------------------------------------------
260
261void LHCOWriter::AnalyseMuons()
262{
263 Muon *element;
264 Track *track;
265 Tower *tower;
266 Jet *jet;
267 Int_t muonCounter, tauCounter, jetCounter, minIndex;
268 Float_t sumPT, sumET, ratET, jetDR, minDR;
269
270 muonCounter = 0;
271 fItMuon->Reset();
272 while((element = static_cast<Muon*>(fItMuon->Next())))
273 {
274 Reset();
275
276 sumPT = 0.0;
277 fItTrack->Reset();
278 while((track = static_cast<Track*>(fItTrack->Next())))
279 {
280 if(element->P4().DeltaR(track->P4()) < 0.5) sumPT += track->PT;
281 }
282
283 sumET = 0.0;
284 fItTower->Reset();
285 while((tower = static_cast<Tower*>(fItTower->Next())))
286 {
287 if(element->P4().DeltaR(tower->P4()) < 0.5) sumET += tower->ET;
288 }
289
290 tauCounter = 0;
291 jetCounter = 0;
292 minIndex = -1;
293 minDR = 1.0E9;
294 fItJet->Reset();
295 while((jet = static_cast<Jet*>(fItJet->Next())))
296 {
297 if(jet->TauTag != 0)
298 {
299 ++tauCounter;
300 continue;
301 }
302
303 jetDR = element->P4().DeltaR(jet->P4());
304 if(jetDR < minDR)
305 {
306 minIndex = jetCounter;
307 minDR = jetDR;
308 }
309 ++jetCounter;
310 }
311
312 fIntParam[1] = 2;
313
314 fDblParam[0] = element->Eta;
315 fDblParam[1] = element->Phi;
316 fDblParam[2] = element->PT;
317
318 fDblParam[3] = 0.11;
319
320 fDblParam[4] = element->Charge;
321
322 if(minIndex >= 0)
323 {
324 fDblParam[5] = fIntParam[0] + fBranchMuon->GetEntriesFast() - muonCounter + tauCounter + minIndex;
325 }
326
327 ratET = sumET/element->PT;
328 fDblParam[6] = Float_t(TMath::Nint(sumPT)) + (ratET < 1.0 ? ratET : 0.99);
329
330 Write();
331 ++muonCounter;
332 }
333}
334
335//---------------------------------------------------------------------------
336
337void LHCOWriter::AnalyseTauJets()
338{
339 Jet *element;
340 Track *track;
341 Int_t counter;
342
343 fItJet->Reset();
344 while((element = static_cast<Jet*>(fItJet->Next())))
345 {
346 if(element->TauTag == 0) continue;
347
348 Reset();
349
350 counter = 0;
351 fItTrack->Reset();
352 while((track = static_cast<Track*>(fItTrack->Next())))
353 {
354 if(element->P4().DeltaR(track->P4()) < 0.5) ++counter;
355 }
356
357 fIntParam[1] = 3;
358
359 fDblParam[0] = element->Eta;
360 fDblParam[1] = element->Phi;
361 fDblParam[2] = element->PT;
362 fDblParam[3] = element->Mass;
363 fDblParam[4] = counter * element->Charge;
364
365 fDblParam[6] = element->EhadOverEem;
366
367 Write();
368 }
369}
370
371//---------------------------------------------------------------------------
372
373void LHCOWriter::AnalyseJets()
374{
375 Jet *element;
376 Track *track;
377 Int_t counter;
378
379 fItJet->Reset();
380 while((element = static_cast<Jet*>(fItJet->Next())))
381 {
382 if(element->TauTag != 0) continue;
383
384 Reset();
385
386 counter = 0;
387 fItTrack->Reset();
388 while((track = static_cast<Track*>(fItTrack->Next())))
389 {
390 if(element->P4().DeltaR(track->P4()) < 0.5) ++counter;
391 }
392
393 fIntParam[1] = 4;
394
395 fDblParam[0] = element->Eta;
396 fDblParam[1] = element->Phi;
397 fDblParam[2] = element->PT;
398 fDblParam[3] = element->Mass;
399 fDblParam[4] = counter;
400 fDblParam[5] = element->BTag;
401 fDblParam[6] = element->EhadOverEem;
402
403 Write();
404 }
405}
406
407//---------------------------------------------------------------------------
408
409void LHCOWriter::AnalyseMissingET()
410{
411 MissingET *element;
412
413 element = static_cast<MissingET*>(fBranchMissingET->At(0));
414
415 Reset();
416
417 fIntParam[1] = 6;
418
419 fDblParam[1] = element->Phi;
420 fDblParam[2] = element->MET;
421
422 Write();
423}
424
425//---------------------------------------------------------------------------
426
427static bool interrupted = false;
428
429void SignalHandler(int sig)
430{
431 interrupted = true;
432}
433
434//---------------------------------------------------------------------------
435
436int main(int argc, char *argv[])
437{
438 char appName[] = "root2lhco";
439 stringstream message;
440 FILE *outputFile = 0;
441 TChain *inputChain = 0;
442 LHCOWriter *writer = 0;
443 ExRootTreeReader *treeReader = 0;
444 Long64_t entry, allEntries;
445
446 if(argc < 2 || argc > 3)
447 {
448 cerr << " Usage: " << appName << " input_file" << " [output_file]" << endl;
449 cerr << " input_file - input file in ROOT format," << endl;
450 cerr << " output_file - output file in LHCO format," << endl;
451 cerr << " with no output_file, or when output_file is -, write to standard output." << endl;
452 return 1;
453 }
454
455 signal(SIGINT, SignalHandler);
456
457 gROOT->SetBatch();
458
459 int appargc = 1;
460 char *appargv[] = {appName};
461 TApplication app(appName, &appargc, appargv);
462
463 try
464 {
465 cerr << "** Reading " << argv[1] << endl;
466 inputChain = new TChain("Delphes");
467 inputChain->Add(argv[1]);
468
469 ExRootTreeReader *treeReader = new ExRootTreeReader(inputChain);
470
471 if(argc == 2 || strcmp(argv[2], "-") == 0)
472 {
473 outputFile = stdout;
474 }
475 else
476 {
477 outputFile = fopen(argv[2], "w");
478
479 if(outputFile == NULL)
480 {
481 message << "can't open " << argv[2];
482 throw runtime_error(message.str());
483 }
484 }
485
486 fprintf(outputFile, " # typ eta phi pt jmas ntrk btag had/em dum1 dum2\n");
487
488 allEntries = treeReader->GetEntries();
489 cerr << "** Input file contains " << allEntries << " events" << endl;
490
491 if(allEntries > 0)
492 {
493 // Create LHC Olympics converter:
494 writer = new LHCOWriter(treeReader, outputFile);
495
496 ExRootProgressBar progressBar(allEntries - 1);
497 // Loop over all events
498 for(entry = 0; entry < allEntries && !interrupted; ++entry)
499 {
500 if(!treeReader->ReadEntry(entry))
501 {
502 cerr << "** ERROR: cannot read event " << entry << endl;
503 break;
504 }
505
506 writer->ProcessEvent();
507
508 progressBar.Update(entry);
509 }
510 progressBar.Finish();
511
512 delete writer;
513 }
514
515 cerr << "** Exiting..." << endl;
516
517 if(outputFile != stdout) fclose(outputFile);
518 delete treeReader;
519 delete inputChain;
520 return 0;
521 }
522 catch(runtime_error &e)
523 {
524 if(writer) delete writer;
525 if(treeReader) delete treeReader;
526 if(inputChain) delete inputChain;
527 cerr << "** ERROR: " << e.what() << endl;
528 return 1;
529 }
530}
531
532
Note: See TracBrowser for help on using the repository browser.