Fork me on GitHub

source: git/classes/DelphesLHEFReader.cc @ b23e468

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

remove pdgParticle->Stable() from readers

  • Property mode set to 100644
File size: 6.8 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  = 1024;
51
52//---------------------------------------------------------------------------
53
54DelphesLHEFReader::DelphesLHEFReader() :
55  fInputFile(0), fBuffer(0), fPDG(0),
56  fEventReady(kFALSE), fEventCounter(-1), fParticleCounter(-1)
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;
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: " << "invalid event format" << endl;
126      return kFALSE;
127    }
128
129    --fEventCounter;
130  }
131  else if(fParticleCounter > 0)
132  {
133    DelphesStream bufferStream(fBuffer);
134
135    rc = bufferStream.ReadInt(fPID)
136      && bufferStream.ReadInt(fStatus)
137      && bufferStream.ReadInt(fM1)
138      && bufferStream.ReadInt(fM2)
139      && bufferStream.ReadInt(fC1)
140      && bufferStream.ReadInt(fC2)
141      && bufferStream.ReadDbl(fPx)
142      && bufferStream.ReadDbl(fPy)
143      && bufferStream.ReadDbl(fPz)
144      && bufferStream.ReadDbl(fE)
145      && bufferStream.ReadDbl(fMass);
146
147    if(!rc)
148    {
149      cerr << "** ERROR: " << "invalid particle format" << endl;
150      return kFALSE;
151    }
152
153    AnalyzeParticle(factory, allParticleOutputArray,
154      stableParticleOutputArray, partonOutputArray);
155
156    --fParticleCounter;
157  }
158  else if(strstr(fBuffer, "<wgt"))
159  {
160    pch = strpbrk(fBuffer, "\"'");
161    if(!pch)
162    {
163      cerr << "** ERROR: " << "invalid weight format" << endl;
164      return kFALSE;
165    }
166
167    DelphesStream idStream(pch + 1);
168    rc = idStream.ReadInt(id);
169
170    pch = strchr(fBuffer, '>');
171    if(!pch)
172    {
173      cerr << "** ERROR: " << "invalid weight format" << endl;
174      return kFALSE;
175    }
176
177    DelphesStream weightStream(pch + 1);
178    rc = weightStream.ReadDbl(weight);
179
180    if(!rc)
181    {
182      cerr << "** ERROR: " << "invalid weight format" << endl;
183      return kFALSE;
184    }
185
186    fWeightList.push_back(make_pair(id, weight));
187  }
188  else if(strstr(fBuffer, "</event>"))
189  {
190    fEventReady = kTRUE;
191  }
192
193  return kTRUE;
194}
195
196//---------------------------------------------------------------------------
197
198void DelphesLHEFReader::AnalyzeEvent(ExRootTreeBranch *branch, long long eventNumber,
199  TStopwatch *readStopWatch, TStopwatch *procStopWatch)
200{
201  LHEFEvent *element;
202
203  element = static_cast<LHEFEvent *>(branch->NewEntry());
204  element->Number = eventNumber;
205
206  element->ProcessID = fProcessID;
207  element->Weight = fWeight;
208  element->ScalePDF = fScalePDF;
209  element->AlphaQED = fAlphaQED;
210  element->AlphaQCD = fAlphaQCD;
211
212  element->ReadTime = readStopWatch->RealTime();
213  element->ProcTime = procStopWatch->RealTime();
214}
215
216//---------------------------------------------------------------------------
217
218void DelphesLHEFReader::AnalyzeWeight(ExRootTreeBranch *branch)
219{
220  LHEFWeight *element;
221  vector< pair< int, double > >::const_iterator itWeightList;
222
223  for(itWeightList = fWeightList.begin(); itWeightList != fWeightList.end(); ++itWeightList)
224  {
225    element = static_cast<LHEFWeight *>(branch->NewEntry());
226
227    element->ID = itWeightList->first;
228    element->Weight = itWeightList->second;
229  }
230}
231
232//---------------------------------------------------------------------------
233
234void DelphesLHEFReader::AnalyzeParticle(DelphesFactory *factory,
235  TObjArray *allParticleOutputArray,
236  TObjArray *stableParticleOutputArray,
237  TObjArray *partonOutputArray)
238{
239  Candidate *candidate;
240  TParticlePDG *pdgParticle;
241  int pdgCode;
242
243  candidate = factory->NewCandidate();
244
245  candidate->PID = fPID;
246  pdgCode = TMath::Abs(candidate->PID);
247
248  candidate->Status = fStatus;
249
250  pdgParticle = fPDG->GetParticle(fPID);
251  candidate->Charge = pdgParticle ? int(pdgParticle->Charge()/3.0) : -999;
252  candidate->Mass = fMass;
253
254  candidate->Momentum.SetPxPyPzE(fPx, fPy, fPz, fE);
255  candidate->Position.SetXYZT(0.0, 0.0, 0.0, 0.0);
256
257  candidate->M1 = fM1 - 1;
258  candidate->M2 = fM2 - 1;
259
260  candidate->D1 = -1;
261  candidate->D2 = -1;
262
263  allParticleOutputArray->Add(candidate);
264
265  if(!pdgParticle) return;
266
267  if(fStatus == 1)
268  {
269    stableParticleOutputArray->Add(candidate);
270  }
271  else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
272  {
273    partonOutputArray->Add(candidate);
274  }
275}
276
277//---------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.