Fork me on GitHub

source: git/converters/root2lhco.cpp@ ca9e119

ImprovedOutputFile Timing dual_readout llp
Last change on this file since ca9e119 was c62695e, checked in by Michele Selvaggi <michele.selvaggi@…>, 9 years ago

force ntrk = 1 for tau's in root2lhco converter

  • 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 = 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
438int main(int argc, char *argv[])
439{
440 char appName[] = "root2lhco";
441 stringstream message;
442 FILE *outputFile = 0;
443 TChain *inputChain = 0;
444 LHCOWriter *writer = 0;
445 ExRootTreeReader *treeReader = 0;
446 Long64_t entry, allEntries;
447
448 if(argc < 2 || argc > 3)
449 {
450 cerr << " Usage: " << appName << " input_file" << " [output_file]" << endl;
451 cerr << " input_file - input file in ROOT format," << endl;
452 cerr << " output_file - output file in LHCO format," << endl;
453 cerr << " with no output_file, or when output_file is -, write to standard output." << endl;
454 return 1;
455 }
456
457 signal(SIGINT, SignalHandler);
458
459 gROOT->SetBatch();
460
461 int appargc = 1;
462 char *appargv[] = {appName};
463 TApplication app(appName, &appargc, appargv);
464
465 try
466 {
467 cerr << "** Reading " << argv[1] << endl;
468 inputChain = new TChain("Delphes");
469 inputChain->Add(argv[1]);
470
471 ExRootTreeReader *treeReader = new ExRootTreeReader(inputChain);
472
473 if(argc == 2 || strcmp(argv[2], "-") == 0)
474 {
475 outputFile = stdout;
476 }
477 else
478 {
479 outputFile = fopen(argv[2], "w");
480
481 if(outputFile == NULL)
482 {
483 message << "can't open " << argv[2];
484 throw runtime_error(message.str());
485 }
486 }
487
488 fprintf(outputFile, " # typ eta phi pt jmas ntrk btag had/em dum1 dum2\n");
489
490 allEntries = treeReader->GetEntries();
491 cerr << "** Input file contains " << allEntries << " events" << endl;
492
493 if(allEntries > 0)
494 {
495 // Create LHC Olympics converter:
496 writer = new LHCOWriter(treeReader, outputFile);
497
498 ExRootProgressBar progressBar(allEntries - 1);
499 // Loop over all events
500 for(entry = 0; entry < allEntries && !interrupted; ++entry)
501 {
502 if(!treeReader->ReadEntry(entry))
503 {
504 cerr << "** ERROR: cannot read event " << entry << endl;
505 break;
506 }
507
508 writer->ProcessEvent();
509
510 progressBar.Update(entry);
511 }
512 progressBar.Finish();
513
514 delete writer;
515 }
516
517 cerr << "** Exiting..." << endl;
518
519 if(outputFile != stdout) fclose(outputFile);
520 delete treeReader;
521 delete inputChain;
522 return 0;
523 }
524 catch(runtime_error &e)
525 {
526 if(writer) delete writer;
527 if(treeReader) delete treeReader;
528 if(inputChain) delete inputChain;
529 cerr << "** ERROR: " << e.what() << endl;
530 return 1;
531 }
532}
533
534
Note: See TracBrowser for help on using the repository browser.