Fork me on GitHub

source: git/converters/root2lhco.cpp@ 127644a

Last change on this file since 127644a was 341014c, checked in by Pavel Demin <pavel-demin@…>, 6 years ago

apply .clang-format to all .h, .cc and .cpp files

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