Fork me on GitHub

source: svn/trunk/classes/DelphesLHEFReader.cc

Last change on this file was 1265, checked in by Pavel Demin, 11 years ago

read reweighting information from LHEF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id Revision Date
File size: 5.7 KB
Line 
1
2/** \class DelphesLHEFReader
3 *
4 * Reads LHEF file
5 *
6 *
7 * $Date: 2013-08-16 14:29:25 +0000 (Fri, 16 Aug 2013) $
8 * $Revision: 1265 $
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 fEventReady(kFALSE), 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 fEventReady = kFALSE;
69 fEventCounter = -1;
70 fParticleCounter = -1;
71 fRwgtList.clear();
72}
73
74//---------------------------------------------------------------------------
75
76bool DelphesLHEFReader::EventReady()
77{
78 return fEventReady;
79}
80
81//---------------------------------------------------------------------------
82
83bool DelphesLHEFReader::ReadBlock(DelphesFactory *factory,
84 TObjArray *allParticleOutputArray,
85 TObjArray *stableParticleOutputArray,
86 TObjArray *partonOutputArray)
87{
88 int rc;
89 char *pch;
90 double weight;
91
92 if(!fgets(fBuffer, kBufferSize, fInputFile)) return kFALSE;
93
94 if(strstr(fBuffer, "<event>"))
95 {
96 Clear();
97 fEventCounter = 1;
98 }
99 else if(fEventCounter > 0)
100 {
101 DelphesStream bufferStream(fBuffer);
102
103 rc = bufferStream.ReadInt(fParticleCounter)
104 && bufferStream.ReadInt(fProcessID)
105 && bufferStream.ReadDbl(fWeight)
106 && bufferStream.ReadDbl(fScalePDF)
107 && bufferStream.ReadDbl(fAlphaQED)
108 && bufferStream.ReadDbl(fAlphaQCD);
109
110 if(!rc)
111 {
112 cerr << "** ERROR: " << "invalid event format" << endl;
113 return kFALSE;
114 }
115
116 --fEventCounter;
117 }
118 else if(fParticleCounter > 0)
119 {
120 DelphesStream bufferStream(fBuffer);
121
122 rc = bufferStream.ReadInt(fPID)
123 && bufferStream.ReadInt(fStatus)
124 && bufferStream.ReadInt(fM1)
125 && bufferStream.ReadInt(fM2)
126 && bufferStream.ReadInt(fC1)
127 && bufferStream.ReadInt(fC2)
128 && bufferStream.ReadDbl(fPx)
129 && bufferStream.ReadDbl(fPy)
130 && bufferStream.ReadDbl(fPz)
131 && bufferStream.ReadDbl(fE)
132 && bufferStream.ReadDbl(fMass);
133
134 if(!rc)
135 {
136 cerr << "** ERROR: " << "invalid particle format" << endl;
137 return kFALSE;
138 }
139
140 AnalyzeParticle(factory, allParticleOutputArray,
141 stableParticleOutputArray, partonOutputArray);
142
143 --fParticleCounter;
144 }
145 else if(strstr(fBuffer, "<wgt"))
146 {
147 pch = strstr(fBuffer, ">");
148 if(!pch)
149 {
150 cerr << "** ERROR: " << "invalid weight format" << endl;
151 return kFALSE;
152 }
153
154 DelphesStream bufferStream(pch + 1);
155 rc = bufferStream.ReadDbl(weight);
156
157 if(!rc)
158 {
159 cerr << "** ERROR: " << "invalid weight format" << endl;
160 return kFALSE;
161 }
162
163 fRwgtList.push_back(weight);
164 }
165 else if(strstr(fBuffer, "</event>"))
166 {
167 fEventReady = kTRUE;
168 }
169
170 return kTRUE;
171}
172
173//---------------------------------------------------------------------------
174
175void DelphesLHEFReader::AnalyzeEvent(ExRootTreeBranch *branch, long long eventNumber,
176 TStopwatch *readStopWatch, TStopwatch *procStopWatch)
177{
178 LHEFEvent *element;
179
180 element = static_cast<LHEFEvent *>(branch->NewEntry());
181 element->Number = eventNumber;
182
183 element->ProcessID = fProcessID;
184 element->Weight = fWeight;
185 element->ScalePDF = fScalePDF;
186 element->AlphaQED = fAlphaQED;
187 element->AlphaQCD = fAlphaQCD;
188
189 element->ReadTime = readStopWatch->RealTime();
190 element->ProcTime = procStopWatch->RealTime();
191}
192
193//---------------------------------------------------------------------------
194
195void DelphesLHEFReader::AnalyzeRwgt(ExRootTreeBranch *branch)
196{
197 Weight *element;
198 vector<double>::const_iterator itRwgtList;
199
200 for(itRwgtList = fRwgtList.begin(); itRwgtList != fRwgtList.end(); ++itRwgtList)
201 {
202 element = static_cast<Weight *>(branch->NewEntry());
203
204 element->Weight = *itRwgtList;
205 }
206}
207
208//---------------------------------------------------------------------------
209
210void DelphesLHEFReader::AnalyzeParticle(DelphesFactory *factory,
211 TObjArray *allParticleOutputArray,
212 TObjArray *stableParticleOutputArray,
213 TObjArray *partonOutputArray)
214{
215 Candidate *candidate;
216 TParticlePDG *pdgParticle;
217 int pdgCode;
218
219 candidate = factory->NewCandidate();
220
221 candidate->PID = fPID;
222 pdgCode = TMath::Abs(candidate->PID);
223
224 candidate->Status = fStatus;
225
226 pdgParticle = fPDG->GetParticle(fPID);
227 candidate->Charge = pdgParticle ? int(pdgParticle->Charge()/3.0) : -999;
228 candidate->Mass = fMass;
229
230 candidate->Momentum.SetPxPyPzE(fPx, fPy, fPz, fE);
231 candidate->Position.SetXYZT(0.0, 0.0, 0.0, 0.0);
232
233 candidate->M1 = fM1 - 1;
234 candidate->M2 = fM2 - 1;
235
236 candidate->D1 = -1;
237 candidate->D2 = -1;
238
239 allParticleOutputArray->Add(candidate);
240
241 if(!pdgParticle) return;
242
243 if(fStatus == 1 && pdgParticle->Stable())
244 {
245 stableParticleOutputArray->Add(candidate);
246 }
247 else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
248 {
249 partonOutputArray->Add(candidate);
250 }
251}
252
253//---------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.