Fork me on GitHub

source: git/classes/DelphesLHEFReader.cc @ be42bf4

ImprovedOutputFileTimingdual_readoutllp
Last change on this file since be42bf4 was be42bf4, checked in by pavel <pavel@…>, 7 years ago

replace <event with <event>

  • Property mode set to 100644
File size: 4.8 KB
Line 
1
2/** \class DelphesLHEFReader
3 *
4 *  Reads LHEF file
5 *
6 *
7 *  $Date$
8 *  $Revision$
9 *
10 *
11 *  \author P. Demin - UCL, Louvain-la-Neuve
12 *
13 */
14
15#include "classes/DelphesLHEFReader.h"
16
17#include <stdexcept>
18#include <iostream>
19#include <sstream>
20
21#include <stdio.h>
22
23#include "TObjArray.h"
24#include "TStopwatch.h"
25#include "TDatabasePDG.h"
26#include "TParticlePDG.h"
27#include "TLorentzVector.h"
28
29#include "classes/DelphesClasses.h"
30#include "classes/DelphesFactory.h"
31#include "classes/DelphesStream.h"
32
33#include "ExRootAnalysis/ExRootTreeBranch.h"
34
35using namespace std;
36
37static const int kBufferSize  = 1024;
38
39//---------------------------------------------------------------------------
40
41DelphesLHEFReader::DelphesLHEFReader() :
42  fInputFile(0), fBuffer(0), fPDG(0),
43  fEventCounter(-1), fParticleCounter(-1)
44{
45  fBuffer = new char[kBufferSize];
46
47  fPDG = TDatabasePDG::Instance();
48}
49
50//---------------------------------------------------------------------------
51
52DelphesLHEFReader::~DelphesLHEFReader()
53{
54  if(fBuffer) delete[] fBuffer;
55}
56
57//---------------------------------------------------------------------------
58
59void DelphesLHEFReader::SetInputFile(FILE *inputFile)
60{
61  fInputFile = inputFile;
62}
63
64//---------------------------------------------------------------------------
65
66void DelphesLHEFReader::Clear()
67{
68  fEventCounter = -1;
69  fParticleCounter = -1;
70}
71
72//---------------------------------------------------------------------------
73
74bool DelphesLHEFReader::EventReady()
75{
76  return (fEventCounter == 0) && (fParticleCounter == 0);
77}
78
79//---------------------------------------------------------------------------
80
81bool DelphesLHEFReader::ReadBlock(DelphesFactory *factory,
82  TObjArray *allParticleOutputArray,
83  TObjArray *stableParticleOutputArray,
84  TObjArray *partonOutputArray)
85{
86  int rc;
87
88  if(!fgets(fBuffer, kBufferSize, fInputFile)) return kFALSE;
89
90  DelphesStream bufferStream(fBuffer);
91
92  if(strstr(fBuffer, "<event>"))
93  {
94    Clear();
95    fEventCounter = 1;
96  }
97  else if(fEventCounter > 0)
98  {
99    rc = bufferStream.ReadInt(fParticleCounter)
100      && bufferStream.ReadInt(fProcessID)
101      && bufferStream.ReadDbl(fWeight)
102      && bufferStream.ReadDbl(fScalePDF)
103      && bufferStream.ReadDbl(fAlphaQED)
104      && bufferStream.ReadDbl(fAlphaQCD);
105
106    if(!rc)
107    {
108      cerr << "** ERROR: " << "invalid event format" << endl;
109      return kFALSE;
110    }
111
112    --fEventCounter;
113  }
114  else if(fParticleCounter > 0)
115  {
116    rc = bufferStream.ReadInt(fPID)
117      && bufferStream.ReadInt(fStatus)
118      && bufferStream.ReadInt(fM1)
119      && bufferStream.ReadInt(fM2)
120      && bufferStream.ReadInt(fC1)
121      && bufferStream.ReadInt(fC2)
122      && bufferStream.ReadDbl(fPx)
123      && bufferStream.ReadDbl(fPy)
124      && bufferStream.ReadDbl(fPz)
125      && bufferStream.ReadDbl(fE)
126      && bufferStream.ReadDbl(fMass);
127
128    if(!rc)
129    {
130      cerr << "** ERROR: " << "invalid particle format" << endl;
131      return kFALSE;
132    }
133
134    AnalyzeParticle(factory, allParticleOutputArray,
135      stableParticleOutputArray, partonOutputArray);
136
137    --fParticleCounter;
138  }
139
140  return kTRUE;
141}
142
143//---------------------------------------------------------------------------
144
145void DelphesLHEFReader::AnalyzeEvent(ExRootTreeBranch *branch, long long eventNumber,
146  TStopwatch *readStopWatch, TStopwatch *procStopWatch)
147{
148  LHEFEvent *element;
149
150  element = static_cast<LHEFEvent *>(branch->NewEntry());
151  element->Number = eventNumber;
152
153  element->ProcessID = fProcessID;
154  element->Weight = fWeight;
155  element->ScalePDF = fScalePDF;
156  element->AlphaQED = fAlphaQED;
157  element->AlphaQCD = fAlphaQCD;
158
159  element->ReadTime = readStopWatch->RealTime();
160  element->ProcTime = procStopWatch->RealTime();
161}
162
163//---------------------------------------------------------------------------
164
165void DelphesLHEFReader::AnalyzeParticle(DelphesFactory *factory,
166  TObjArray *allParticleOutputArray,
167  TObjArray *stableParticleOutputArray,
168  TObjArray *partonOutputArray)
169{
170  Candidate *candidate;
171  TParticlePDG *pdgParticle;
172  int pdgCode;
173
174  candidate = factory->NewCandidate();
175
176  candidate->PID = fPID;
177  pdgCode = TMath::Abs(candidate->PID);
178
179  candidate->Status = fStatus;
180
181  pdgParticle = fPDG->GetParticle(fPID);
182  candidate->Charge = pdgParticle ? int(pdgParticle->Charge()/3.0) : -999;
183  candidate->Mass = pdgParticle ? pdgParticle->Mass() : -999.9;
184
185  candidate->Momentum.SetPxPyPzE(fPx, fPy, fPz, fE);
186  candidate->Position.SetXYZT(0.0, 0.0, 0.0, 0.0);
187
188  candidate->M1 = fM1 - 1;
189  candidate->M2 = fM2 - 1;
190
191  candidate->D1 = -1;
192  candidate->D2 = -1;
193
194  allParticleOutputArray->Add(candidate);
195
196  if(!pdgParticle) return;
197
198  if(fStatus == 1 && pdgParticle->Stable())
199  {
200    stableParticleOutputArray->Add(candidate);
201  }
202  else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
203  {
204    partonOutputArray->Add(candidate);
205  }
206}
207
208//---------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.