Fork me on GitHub

source: git/classes/DelphesLHEFReader.cc @ 1d9c62a

ImprovedOutputFileTimingdual_readoutllp
Last change on this file since 1d9c62a was 1d9c62a, checked in by Pavel Demin <pavel-demin@…>, 2 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.