source: HiSusy/trunk/Delphes-3.0.0/converters/root2lhco.cpp @ 1

Last change on this file since 1 was 1, checked in by zerwas, 11 years ago

first import of structure, PYTHIA8 and DELPHES

File size: 10.4 KB
Line 
1#include <stdexcept>
2#include <iostream>
3#include <fstream>
4#include <sstream>
5#include <string>
6
7#include <stdlib.h>
8#include <signal.h>
9#include <stdio.h>
10
11#include "TROOT.h"
12#include "TApplication.h"
13
14#include "TFile.h"
15#include "TClonesArray.h"
16
17#include "classes/DelphesClasses.h"
18
19#include "ExRootAnalysis/ExRootTreeReader.h"
20#include "ExRootAnalysis/ExRootProgressBar.h"
21
22using namespace std;
23
24/*
25LHC Olympics format discription from http://www.jthaler.net/olympicswiki/doku.php?id=lhc_olympics:data_file_format
26
27    * The first column of each row is just a counter that labels the object.
28    * 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).
29    * 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].
30    * The next three columns give the pseudorapidity, the azimuthal angle, and the transverse momentum of the object.
31    * The sixth column gives the invariant mass of the object.
32    * 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.
33    * 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.
34    * 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.
35*/
36
37class LHCOWriter
38{
39public:
40  LHCOWriter(ExRootTreeReader *treeReader, FILE *outputFile);
41  ~LHCOWriter();
42
43  void ProcessEvent();
44
45private:
46
47  void Reset();
48  void Write();
49
50  void AnalyseEvent();
51
52  void AnalysePhotons();
53  void AnalyseElectrons();
54  void AnalyseMuons();
55  void AnalyseTauJets();
56  void AnalyseJets();
57
58  void AnalyseMissingET();
59
60  enum {kIntParamSize = 2, kDblParamSize = 9};
61  Int_t fIntParam[kIntParamSize];
62  Double_t fDblParam[kDblParamSize];
63 
64  Long64_t fTriggerWord, fEventNumber;
65
66  ExRootTreeReader *fTreeReader;
67  FILE *fOutputFile;
68
69  TClonesArray *fBranchEvent;
70  TClonesArray *fBranchPhoton;
71  TClonesArray *fBranchElectron;
72  TClonesArray *fBranchMuon;
73  TClonesArray *fBranchTauJet;
74  TClonesArray *fBranchJet;
75  TClonesArray *fBranchMissingET;
76
77  TIterator *fItPhoton;
78  TIterator *fItElectron;
79  TIterator *fItMuon;
80  TIterator *fItTauJet;
81  TIterator *fItJet;
82
83};
84
85//------------------------------------------------------------------------------
86
87LHCOWriter::LHCOWriter(ExRootTreeReader *treeReader, FILE *outputFile) :
88  fTriggerWord(0), fEventNumber(1), fTreeReader(0), fOutputFile(0)
89{
90  fTreeReader = treeReader;
91  fOutputFile = outputFile;
92
93  // information about reconstructed event
94  fBranchEvent = fTreeReader->UseBranch("Event");
95  // reconstructed photons
96  fBranchPhoton = fTreeReader->UseBranch("Photon");
97  fItPhoton = fBranchPhoton->MakeIterator();
98  // reconstructed electrons
99  fBranchElectron = fTreeReader->UseBranch("Electron");
100  fItElectron = fBranchElectron->MakeIterator();
101  // reconstructed muons
102  fBranchMuon = fTreeReader->UseBranch("Muon");
103  fItMuon = fBranchMuon->MakeIterator();
104  // reconstructed hadronically-decaying tau leptons
105  fBranchTauJet = fTreeReader->UseBranch("TauJet");
106  fItTauJet = fBranchTauJet->MakeIterator();
107  // reconstructed jets
108  fBranchJet = fTreeReader->UseBranch("Jet");
109  fItJet = fBranchJet->MakeIterator();
110  // missing transverse energy
111  fBranchMissingET = fTreeReader->UseBranch("MissingET");
112}
113
114//------------------------------------------------------------------------------
115
116LHCOWriter::~LHCOWriter()
117{
118}
119
120//---------------------------------------------------------------------------
121
122void LHCOWriter::ProcessEvent()
123{
124  fIntParam[0] = 0;
125
126  AnalyseEvent();
127
128  AnalysePhotons();
129  AnalyseElectrons();
130  AnalyseMuons();
131  AnalyseTauJets();
132  AnalyseJets();
133
134  AnalyseMissingET();
135}
136
137//---------------------------------------------------------------------------
138
139void LHCOWriter::Reset()
140{
141  int i;
142  for(i = 1; i < kIntParamSize; ++i)
143  {
144    fIntParam[i] = 0;
145  }
146
147  for(i = 0; i < kDblParamSize; ++i)
148  {
149    fDblParam[i] = 0.0;
150  }
151}
152
153//---------------------------------------------------------------------------
154
155void LHCOWriter::Write()
156{
157  fprintf(fOutputFile, "%d\t%d\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\n",
158    fIntParam[0], fIntParam[1], fDblParam[0], fDblParam[1], fDblParam[2],
159    fDblParam[3], fDblParam[4], fDblParam[5], fDblParam[6], fDblParam[7], fDblParam[8]);
160
161  ++fIntParam[0];
162}
163
164//---------------------------------------------------------------------------
165
166void LHCOWriter::AnalyseEvent()
167{
168  Event *element;
169
170  element = static_cast<Event*>(fBranchEvent->At(0));
171
172  fprintf(fOutputFile, "0\t%lld\t0\n", element->Number);
173
174  ++fIntParam[0];
175}
176
177//---------------------------------------------------------------------------
178
179void LHCOWriter::AnalysePhotons()
180{
181  Photon *element;
182
183  fItPhoton->Reset();
184  while((element = static_cast<Photon*>(fItPhoton->Next())))
185  {
186    Reset();
187
188    fIntParam[1] = 0;
189
190    fDblParam[0] = element->Eta;
191    fDblParam[1] = element->Phi;
192    fDblParam[2] = element->PT;
193
194    fDblParam[6] = 0.0; // element->EhadOverEem;
195
196    Write();
197  }
198}
199
200//---------------------------------------------------------------------------
201
202void LHCOWriter::AnalyseElectrons()
203{
204  Electron *element;
205
206  fItElectron->Reset();
207  while((element = static_cast<Electron*>(fItElectron->Next())))
208  {
209    Reset();
210
211    fIntParam[1] = 1;
212
213    fDblParam[0] = element->Eta;
214    fDblParam[1] = element->Phi;
215    fDblParam[2] = element->PT;
216
217    fDblParam[4] = 0.0; // element->Ntrk * element->Charge;
218
219    fDblParam[6] = 0.0; // element->EhadOverEem;
220
221    Write();
222  }
223}
224
225//---------------------------------------------------------------------------
226
227void LHCOWriter::AnalyseMuons()
228{
229  Muon *element;
230
231  fItMuon->Reset();
232  while((element = static_cast<Muon*>(fItMuon->Next())))
233  {
234    Reset();
235
236    fIntParam[1] = 2;
237
238    fDblParam[0] = element->Eta;
239    fDblParam[1] = element->Phi;
240    fDblParam[2] = element->PT;
241   
242    fDblParam[4] = 0.0; // element->Ntrk * element->Charge;
243
244    fDblParam[5] = 0.0; // element->JetIndex;
245
246    fDblParam[6] = 0.0; // element->PTiso + element->ETiso;
247
248    Write();
249  }
250}
251
252//---------------------------------------------------------------------------
253
254void LHCOWriter::AnalyseTauJets()
255{
256  TauJet *element;
257
258  fItTauJet->Reset();
259  while((element = static_cast<TauJet*>(fItTauJet->Next())))
260  {
261    Reset();
262
263    fIntParam[1] = 3;
264
265    fDblParam[0] = element->Eta;
266    fDblParam[1] = element->Phi;
267    fDblParam[2] = element->PT;
268   
269    fDblParam[4] = 0.0; // element->Ntrk * element->Charge;
270
271    fDblParam[6] = 0.0; // element->EhadOverEem;
272   
273    Write();
274  }
275}
276
277//---------------------------------------------------------------------------
278
279void LHCOWriter::AnalyseJets()
280{
281  Jet *element;
282
283  fItJet->Reset();
284  while((element = static_cast<Jet*>(fItJet->Next())))
285  {
286    Reset();
287
288    fIntParam[1] = 4;
289
290    fDblParam[0] = element->Eta;
291    fDblParam[1] = element->Phi;
292    fDblParam[2] = element->PT;
293   
294    fDblParam[3] = element->Mass;
295
296    fDblParam[4] = element->Ntrk;
297   
298    fDblParam[5] = element->BTag;
299
300    fDblParam[6] = 0.0; // element->EhadOverEem;
301   
302    Write();
303  }
304}
305
306//---------------------------------------------------------------------------
307
308void LHCOWriter::AnalyseMissingET()
309{
310  MissingET *element;
311
312  element = static_cast<MissingET*>(fBranchMissingET->At(0));
313
314  Reset();
315
316  fIntParam[1] = 6;
317
318  fDblParam[1] = element->Phi;
319  fDblParam[2] = element->MET;
320
321  Write();
322}
323
324//---------------------------------------------------------------------------
325
326static bool interrupted = false;
327
328void SignalHandler(int sig)
329{
330        interrupted = true;
331}
332
333//---------------------------------------------------------------------------
334
335int main(int argc, char *argv[])
336{
337  char appName[] = "root2lhco";
338  stringstream message;
339  FILE *outputFile = 0;
340  TChain *inputChain = 0;
341  LHCOWriter *writer = 0;
342  ExRootTreeReader *treeReader = 0;
343  Long64_t entry, allEntries;
344
345  if(argc < 2 || argc > 3)
346  {
347    cout << " Usage: " << appName << " input_file" << " [output_file]" << endl;
348    cout << " input_file - input file in ROOT format," << endl;
349    cout << " output_file - output file in LHCO format," << endl;
350    cout << " with no output_file, or when output_file is -, write to standard output." << endl;
351    return 1;
352  }
353
354  signal(SIGINT, SignalHandler);
355 
356  gROOT->SetBatch();
357 
358  int appargc = 1;
359  char *appargv[] = {appName};
360  TApplication app(appName, &appargc, appargv);
361
362  try
363  {
364    cout << "** Reading " << argv[1] << endl;
365    inputChain = new TChain("Delphes");
366    inputChain->Add(argv[1]);
367
368    ExRootTreeReader *treeReader = new ExRootTreeReader(inputChain);
369
370    if(argc == 2 || strcmp(argv[2], "-") == 0)
371    {
372      outputFile = stdout;
373    }
374    else
375    {
376      outputFile = fopen(argv[2], "w");
377
378      if(outputFile == NULL)
379      {
380        message << "can't open " << argv[2];
381        throw runtime_error(message.str());
382      }
383    }
384
385    allEntries = treeReader->GetEntries();
386    cerr << "** Input file contains " << allEntries << " events" << endl;
387
388    if(allEntries > 0)
389    {
390      // Create LHC Olympics converter:
391      writer = new LHCOWriter(treeReader, outputFile);
392
393      ExRootProgressBar progressBar(allEntries - 1);
394      // Loop over all events
395      for(entry = 0; entry < allEntries && !interrupted; ++entry)
396      {
397        if(!treeReader->ReadEntry(entry))
398        {
399          cout << "** ERROR: cannot read event " << entry << endl;
400          break;
401        }
402   
403        writer->ProcessEvent();
404
405        progressBar.Update(entry);
406      }
407      progressBar.Finish();
408     
409      delete writer;
410    }
411
412    cout << "** Exiting..." << endl;
413
414    if(outputFile != stdout) fclose(outputFile);
415    delete treeReader;
416    delete inputChain;
417    return 0;
418  }
419  catch(runtime_error &e)
420  {
421    if(writer) delete writer;
422    if(treeReader) delete treeReader;
423    if(inputChain) delete inputChain;
424    cerr << "** ERROR: " << e.what() << endl;
425    return 1;
426  }
427}
428
429
Note: See TracBrowser for help on using the repository browser.