Fork me on GitHub

source: git/classes/DelphesLHEFReader.cc@ 224c929

Last change on this file since 224c929 was 77e9ae1, checked in by Pavel Demin <pavel-demin@…>, 6 years ago

set Standard to Cpp03 in .clang-format

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