Fork me on GitHub

source: git/classes/DelphesLHEFReader.cc@ 0518688

ImprovedOutputFile Timing llp
Last change on this file since 0518688 was 1d9c62a, checked in by Pavel Demin <pavel-demin@…>, 6 years ago

fix DelphesLHEFReader.cc compilation under Ubuntu 14.04

  • Property mode set to 100644
File size: 7.4 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
20/** \class DelphesLHEFReader
21 *
22 * Reads LHEF file
23 *
24 * \author P. Demin - UCL, Louvain-la-Neuve
25 *
26 */
27
28#include "classes/DelphesLHEFReader.h"
29
30#include <stdexcept>
31#include <iostream>
32#include <sstream>
33
34#include <stdio.h>
35
36#include "TObjArray.h"
37#include "TStopwatch.h"
38#include "TDatabasePDG.h"
39#include "TParticlePDG.h"
40#include "TLorentzVector.h"
41
42#include "classes/DelphesClasses.h"
43#include "classes/DelphesFactory.h"
44#include "classes/DelphesStream.h"
45
46#include "ExRootAnalysis/ExRootTreeBranch.h"
47
48using namespace std;
49
50static const int kBufferSize = 16384;
51
52//---------------------------------------------------------------------------
53
54DelphesLHEFReader::DelphesLHEFReader() :
55 fInputFile(0), fBuffer(0), fPDG(0),
56 fEventReady(kFALSE), fEventCounter(-1), fParticleCounter(-1), fCrossSection(1)
57
58{
59 fBuffer = new char[kBufferSize];
60
61 fPDG = TDatabasePDG::Instance();
62}
63
64//---------------------------------------------------------------------------
65
66DelphesLHEFReader::~DelphesLHEFReader()
67{
68 if(fBuffer) delete[] fBuffer;
69}
70
71//---------------------------------------------------------------------------
72
73void DelphesLHEFReader::SetInputFile(FILE *inputFile)
74{
75 fInputFile = inputFile;
76}
77
78//---------------------------------------------------------------------------
79
80void DelphesLHEFReader::Clear()
81{
82 fEventReady = kFALSE;
83 fEventCounter = -1;
84 fParticleCounter = -1;
85 fWeightList.clear();
86}
87
88//---------------------------------------------------------------------------
89
90bool DelphesLHEFReader::EventReady()
91{
92 return fEventReady;
93}
94
95//---------------------------------------------------------------------------
96
97bool DelphesLHEFReader::ReadBlock(DelphesFactory *factory,
98 TObjArray *allParticleOutputArray,
99 TObjArray *stableParticleOutputArray,
100 TObjArray *partonOutputArray)
101{
102 int rc, id;
103 char *pch;
104 double weight, xsec;
105
106 if(!fgets(fBuffer, kBufferSize, fInputFile)) return kFALSE;
107
108 if(strstr(fBuffer, "<event>"))
109 {
110 Clear();
111 fEventCounter = 1;
112 }
113 else if(fEventCounter > 0)
114 {
115 DelphesStream bufferStream(fBuffer);
116
117 rc = bufferStream.ReadInt(fParticleCounter)
118 && bufferStream.ReadInt(fProcessID)
119 && bufferStream.ReadDbl(fWeight)
120 && bufferStream.ReadDbl(fScalePDF)
121 && bufferStream.ReadDbl(fAlphaQED)
122 && bufferStream.ReadDbl(fAlphaQCD);
123
124 if(!rc)
125 {
126 cerr << "** ERROR: " << "invalid event format" << endl;
127 return kFALSE;
128 }
129
130 --fEventCounter;
131 }
132 else if(fParticleCounter > 0)
133 {
134 DelphesStream bufferStream(fBuffer);
135
136 rc = bufferStream.ReadInt(fPID)
137 && bufferStream.ReadInt(fStatus)
138 && bufferStream.ReadInt(fM1)
139 && bufferStream.ReadInt(fM2)
140 && bufferStream.ReadInt(fC1)
141 && bufferStream.ReadInt(fC2)
142 && bufferStream.ReadDbl(fPx)
143 && bufferStream.ReadDbl(fPy)
144 && bufferStream.ReadDbl(fPz)
145 && bufferStream.ReadDbl(fE)
146 && bufferStream.ReadDbl(fMass);
147
148 if(!rc)
149 {
150 cerr << "** ERROR: " << "invalid particle format" << endl;
151 return kFALSE;
152 }
153
154 AnalyzeParticle(factory, allParticleOutputArray,
155 stableParticleOutputArray, partonOutputArray);
156
157 --fParticleCounter;
158 }
159 else if(strstr(fBuffer, "<wgt"))
160 {
161 pch = strpbrk(fBuffer, "\"'");
162 if(!pch)
163 {
164 cerr << "** ERROR: " << "invalid weight format" << endl;
165 return kFALSE;
166 }
167
168 DelphesStream idStream(pch + 1);
169 rc = idStream.ReadInt(id);
170
171 pch = strchr(fBuffer, '>');
172 if(!pch)
173 {
174 cerr << "** ERROR: " << "invalid weight format" << endl;
175 return kFALSE;
176 }
177
178 DelphesStream weightStream(pch + 1);
179 rc = weightStream.ReadDbl(weight);
180
181 if(!rc)
182 {
183 cerr << "** ERROR: " << "invalid weight format" << endl;
184 return kFALSE;
185 }
186
187 fWeightList.push_back(make_pair(id, weight));
188 }
189 else if(strstr(fBuffer, "<xsecinfo"))
190 {
191 pch = strstr(fBuffer, "totxsec");
192 if(!pch)
193 {
194 cerr << "** ERROR: " << "invalid cross section format" << endl;
195 return kFALSE;
196 }
197
198 pch = strpbrk(pch + 1, "\"'");
199 if(!pch)
200 {
201 cerr << "** ERROR: " << "invalid cross section format" << endl;
202 return kFALSE;
203 }
204
205 DelphesStream xsecStream(pch + 1);
206 rc = xsecStream.ReadDbl(xsec);
207
208 if(!rc)
209 {
210 cerr << "** ERROR: " << "invalid cross section format" << endl;
211 return kFALSE;
212 }
213
214 fCrossSection = xsec;
215 }
216 else if(strstr(fBuffer, "</event>"))
217 {
218 fEventReady = kTRUE;
219 }
220
221 return kTRUE;
222}
223
224//---------------------------------------------------------------------------
225
226void DelphesLHEFReader::AnalyzeEvent(ExRootTreeBranch *branch, long long eventNumber,
227 TStopwatch *readStopWatch, TStopwatch *procStopWatch)
228{
229 LHEFEvent *element;
230
231 element = static_cast<LHEFEvent *>(branch->NewEntry());
232 element->Number = eventNumber;
233
234 element->ProcessID = fProcessID;
235 element->Weight = fWeight;
236 element->CrossSection = fCrossSection;
237
238 element->ScalePDF = fScalePDF;
239 element->AlphaQED = fAlphaQED;
240 element->AlphaQCD = fAlphaQCD;
241
242 element->ReadTime = readStopWatch->RealTime();
243 element->ProcTime = procStopWatch->RealTime();
244}
245
246//---------------------------------------------------------------------------
247
248void DelphesLHEFReader::AnalyzeWeight(ExRootTreeBranch *branch)
249{
250 LHEFWeight *element;
251 vector< pair< int, double > >::const_iterator itWeightList;
252
253 for(itWeightList = fWeightList.begin(); itWeightList != fWeightList.end(); ++itWeightList)
254 {
255 element = static_cast<LHEFWeight *>(branch->NewEntry());
256
257 element->ID = itWeightList->first;
258 element->Weight = itWeightList->second;
259 }
260}
261
262//---------------------------------------------------------------------------
263
264void DelphesLHEFReader::AnalyzeParticle(DelphesFactory *factory,
265 TObjArray *allParticleOutputArray,
266 TObjArray *stableParticleOutputArray,
267 TObjArray *partonOutputArray)
268{
269 Candidate *candidate;
270 TParticlePDG *pdgParticle;
271 int pdgCode;
272
273 candidate = factory->NewCandidate();
274
275 candidate->PID = fPID;
276 pdgCode = TMath::Abs(candidate->PID);
277
278 candidate->Status = fStatus;
279
280 pdgParticle = fPDG->GetParticle(fPID);
281 candidate->Charge = pdgParticle ? int(pdgParticle->Charge()/3.0) : -999;
282 candidate->Mass = fMass;
283
284 candidate->Momentum.SetPxPyPzE(fPx, fPy, fPz, fE);
285 candidate->Position.SetXYZT(0.0, 0.0, 0.0, 0.0);
286
287 candidate->M1 = fM1 - 1;
288 candidate->M2 = fM2 - 1;
289
290 candidate->D1 = -1;
291 candidate->D2 = -1;
292
293 allParticleOutputArray->Add(candidate);
294
295 if(!pdgParticle) return;
296
297 if(fStatus == 1)
298 {
299 stableParticleOutputArray->Add(candidate);
300 }
301 else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
302 {
303 partonOutputArray->Add(candidate);
304 }
305}
306
307//---------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.