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