Fork me on GitHub

source: git/converters/root2lhco.cpp@ 0882413

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 0882413 was c3547ed, checked in by Pavel Demin <pavel-demin@…>, 6 years ago

improve argument parsing in root2lhco.cpp

  • Property mode set to 100644
File size: 14.2 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, string JetBranchName);
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, string jetBranchName) :
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(jetBranchName.c_str());
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 = 1;
351
352 /*
353 fItTrack->Reset();
354 while((track = static_cast<Track*>(fItTrack->Next())))
355 {
356 if(element->P4().DeltaR(track->P4()) < 0.5) ++counter;
357 }
358 */
359 fIntParam[1] = 3;
360
361 fDblParam[0] = element->Eta;
362 fDblParam[1] = element->Phi;
363 fDblParam[2] = element->PT;
364 fDblParam[3] = element->Mass;
365 fDblParam[4] = counter * element->Charge;
366
367 fDblParam[6] = element->EhadOverEem;
368
369 Write();
370 }
371}
372
373//---------------------------------------------------------------------------
374
375void LHCOWriter::AnalyseJets()
376{
377 Jet *element;
378 Track *track;
379 Int_t counter;
380
381 fItJet->Reset();
382 while((element = static_cast<Jet*>(fItJet->Next())))
383 {
384 if(element->TauTag != 0) continue;
385
386 Reset();
387
388 counter = 0;
389 fItTrack->Reset();
390 while((track = static_cast<Track*>(fItTrack->Next())))
391 {
392 if(element->P4().DeltaR(track->P4()) < 0.5) ++counter;
393 }
394
395 fIntParam[1] = 4;
396
397 fDblParam[0] = element->Eta;
398 fDblParam[1] = element->Phi;
399 fDblParam[2] = element->PT;
400 fDblParam[3] = element->Mass;
401 fDblParam[4] = counter;
402 fDblParam[5] = element->BTag;
403 fDblParam[6] = element->EhadOverEem;
404
405 Write();
406 }
407}
408
409//---------------------------------------------------------------------------
410
411void LHCOWriter::AnalyseMissingET()
412{
413 MissingET *element;
414
415 element = static_cast<MissingET*>(fBranchMissingET->At(0));
416
417 Reset();
418
419 fIntParam[1] = 6;
420
421 fDblParam[1] = element->Phi;
422 fDblParam[2] = element->MET;
423
424 Write();
425}
426
427//---------------------------------------------------------------------------
428
429static bool interrupted = false;
430
431void SignalHandler(int sig)
432{
433 interrupted = true;
434}
435
436//---------------------------------------------------------------------------
437
438vector<string> ArgSplitter(char *arg)
439{
440 string s = arg;
441 string delimiter = "=";
442 vector<string> result;
443 size_t first = 0, last = 0;
444
445 while((last = s.find(delimiter, first)) != string::npos)
446 {
447 result.push_back(s.substr(first, last - first));
448 first = last + delimiter.length();
449 }
450
451 result.push_back(s.substr(first, last));
452
453 return result;
454}
455
456//---------------------------------------------------------------------------
457
458int main(int argc, char *argv[])
459{
460 int i, j;
461 char appName[] = "root2lhco";
462 stringstream message;
463 FILE *outputFile = 0;
464 TChain *inputChain = 0;
465 LHCOWriter *writer = 0;
466 ExRootTreeReader *treeReader = 0;
467 Long64_t entry, allEntries;
468 string jetBranchName = "Jet";
469 vector<string> result;
470
471 if(argc < 2 || argc > 4)
472 {
473 cerr << " Usage: " << appName << " input_file" << " [output_file] [--jet-branch=Jet]" << endl;
474 cerr << " input_file - input file in ROOT format," << endl;
475 cerr << " output_file - output file in LHCO format," << endl;
476 cerr << " with no output_file, or when output_file is -, write to standard output." << endl;
477 return 1;
478 }
479
480 i = 2;
481 while(i < argc)
482 {
483 result = ArgSplitter(argv[i]);
484
485 if(result.size() == 2 && result[0] == "--jet-branch")
486 {
487 jetBranchName = result[1];
488 for(j = i + 1; j < argc; ++j) argv[j - 1] = argv[j];
489 --argc;
490 }
491 else
492 {
493 ++i;
494 }
495 }
496 cerr << "** Using the jet branch named " << jetBranchName << endl;
497
498 signal(SIGINT, SignalHandler);
499
500 gROOT->SetBatch();
501
502 int appargc = 1;
503 char *appargv[] = {appName};
504 TApplication app(appName, &appargc, appargv);
505
506 try
507 {
508 cerr << "** Reading " << argv[1] << endl;
509 inputChain = new TChain("Delphes");
510 inputChain->Add(argv[1]);
511
512 ExRootTreeReader *treeReader = new ExRootTreeReader(inputChain);
513
514 if(argc == 2 || strcmp(argv[2], "-") == 0)
515 {
516 outputFile = stdout;
517 }
518 else
519 {
520 outputFile = fopen(argv[2], "w");
521
522 if(outputFile == NULL)
523 {
524 message << "can't open " << argv[2];
525 throw runtime_error(message.str());
526 }
527 }
528
529 fprintf(outputFile, " # typ eta phi pt jmas ntrk btag had/em dum1 dum2\n");
530
531 allEntries = treeReader->GetEntries();
532 cerr << "** Input file contains " << allEntries << " events" << endl;
533
534 if(allEntries > 0)
535 {
536 // Create LHC Olympics converter:
537 writer = new LHCOWriter(treeReader, outputFile, jetBranchName);
538
539 ExRootProgressBar progressBar(allEntries - 1);
540 // Loop over all events
541 for(entry = 0; entry < allEntries && !interrupted; ++entry)
542 {
543 if(!treeReader->ReadEntry(entry))
544 {
545 cerr << "** ERROR: cannot read event " << entry << endl;
546 break;
547 }
548
549 writer->ProcessEvent();
550
551 progressBar.Update(entry);
552 }
553 progressBar.Finish();
554
555 delete writer;
556 }
557
558 cerr << "** Exiting..." << endl;
559
560 if(outputFile != stdout) fclose(outputFile);
561 delete treeReader;
562 delete inputChain;
563 return 0;
564 }
565 catch(runtime_error &e)
566 {
567 if(writer) delete writer;
568 if(treeReader) delete treeReader;
569 if(inputChain) delete inputChain;
570 cerr << "** ERROR: " << e.what() << endl;
571 return 1;
572 }
573}
Note: See TracBrowser for help on using the repository browser.