Fork me on GitHub

source: git/classes/DelphesLHEFReader.cc@ 3db5282

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 3db5282 was cab38f6, checked in by Pavel Demin <pavel.demin@…>, 10 years ago

remove svn tags and fix formatting

  • 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 * \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 fRwgtList.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;
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 = strstr(fBuffer, ">");
161 if(!pch)
162 {
163 cerr << "** ERROR: " << "invalid weight format" << endl;
164 return kFALSE;
165 }
166
167 DelphesStream bufferStream(pch + 1);
168 rc = bufferStream.ReadDbl(weight);
169
170 if(!rc)
171 {
172 cerr << "** ERROR: " << "invalid weight format" << endl;
173 return kFALSE;
174 }
175
176 fRwgtList.push_back(weight);
177 }
178 else if(strstr(fBuffer, "</event>"))
179 {
180 fEventReady = kTRUE;
181 }
182
183 return kTRUE;
184}
185
186//---------------------------------------------------------------------------
187
188void DelphesLHEFReader::AnalyzeEvent(ExRootTreeBranch *branch, long long eventNumber,
189 TStopwatch *readStopWatch, TStopwatch *procStopWatch)
190{
191 LHEFEvent *element;
192
193 element = static_cast<LHEFEvent *>(branch->NewEntry());
194 element->Number = eventNumber;
195
196 element->ProcessID = fProcessID;
197 element->Weight = fWeight;
198 element->ScalePDF = fScalePDF;
199 element->AlphaQED = fAlphaQED;
200 element->AlphaQCD = fAlphaQCD;
201
202 element->ReadTime = readStopWatch->RealTime();
203 element->ProcTime = procStopWatch->RealTime();
204}
205
206//---------------------------------------------------------------------------
207
208void DelphesLHEFReader::AnalyzeRwgt(ExRootTreeBranch *branch)
209{
210 Weight *element;
211 vector<double>::const_iterator itRwgtList;
212
213 for(itRwgtList = fRwgtList.begin(); itRwgtList != fRwgtList.end(); ++itRwgtList)
214 {
215 element = static_cast<Weight *>(branch->NewEntry());
216
217 element->Weight = *itRwgtList;
218 }
219}
220
221//---------------------------------------------------------------------------
222
223void DelphesLHEFReader::AnalyzeParticle(DelphesFactory *factory,
224 TObjArray *allParticleOutputArray,
225 TObjArray *stableParticleOutputArray,
226 TObjArray *partonOutputArray)
227{
228 Candidate *candidate;
229 TParticlePDG *pdgParticle;
230 int pdgCode;
231
232 candidate = factory->NewCandidate();
233
234 candidate->PID = fPID;
235 pdgCode = TMath::Abs(candidate->PID);
236
237 candidate->Status = fStatus;
238
239 pdgParticle = fPDG->GetParticle(fPID);
240 candidate->Charge = pdgParticle ? int(pdgParticle->Charge()/3.0) : -999;
241 candidate->Mass = fMass;
242
243 candidate->Momentum.SetPxPyPzE(fPx, fPy, fPz, fE);
244 candidate->Position.SetXYZT(0.0, 0.0, 0.0, 0.0);
245
246 candidate->M1 = fM1 - 1;
247 candidate->M2 = fM2 - 1;
248
249 candidate->D1 = -1;
250 candidate->D2 = -1;
251
252 allParticleOutputArray->Add(candidate);
253
254 if(!pdgParticle) return;
255
256 if(fStatus == 1 && pdgParticle->Stable())
257 {
258 stableParticleOutputArray->Add(candidate);
259 }
260 else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
261 {
262 partonOutputArray->Add(candidate);
263 }
264}
265
266//---------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.