Fork me on GitHub

source: git/converters/root2lhco.cpp@ f3fcb6e

ImprovedOutputFile Timing dual_readout llp
Last change on this file since f3fcb6e was 7c0fcd5, checked in by Pavel Demin <demin@…>, 10 years ago

delete duplicate license file and prepend GPLv3 header to all source code files

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