Fork me on GitHub

source: git/classes/DelphesLHEFReader.cc @ e79c954

Last change on this file since e79c954 was 77e9ae1, checked in by Pavel Demin <pavel-demin@…>, 19 months 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.