Fork me on GitHub

source: git/classes/DelphesLHEFReader.cc @ b443089

ImprovedOutputFileTimingdual_readoutllp
Last change on this file since b443089 was b443089, checked in by Pavel Demin <pavel.demin@…>, 6 years ago

fix EOL characters in GPLv3 header

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