Fork me on GitHub

source: git/classes/DelphesHepMCReader.cc@ 7993cad

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 7993cad was b23e468, checked in by Pavel Demin <pavel.demin@…>, 9 years ago

remove pdgParticle->Stable() from readers
  • Property mode set to 100644
File size: 10.7 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 DelphesHepMCReader
21 *
22 * Reads HepMC file
23 *
24 * \author P. Demin - UCL, Louvain-la-Neuve
25 *
26 */
27
28#include "classes/DelphesHepMCReader.h"
29
30#include <stdexcept>
31#include <iostream>
32#include <sstream>
33
34#include <map>
35#include <vector>
36
37#include <stdio.h>
38
39#include "TObjArray.h"
40#include "TStopwatch.h"
41#include "TDatabasePDG.h"
42#include "TParticlePDG.h"
43#include "TLorentzVector.h"
44
45#include "classes/DelphesClasses.h"
46#include "classes/DelphesFactory.h"
47#include "classes/DelphesStream.h"
48
49#include "ExRootAnalysis/ExRootTreeBranch.h"
50
51using namespace std;
52
53static const int kBufferSize = 1024;
54
55//---------------------------------------------------------------------------
56
57DelphesHepMCReader::DelphesHepMCReader() :
58 fInputFile(0), fBuffer(0), fPDG(0),
59 fVertexCounter(-1), fInCounter(-1), fOutCounter(-1),
60 fParticleCounter(0)
61{
62 fBuffer = new char[kBufferSize];
63
64 fPDG = TDatabasePDG::Instance();
65}
66
67//---------------------------------------------------------------------------
68
69DelphesHepMCReader::~DelphesHepMCReader()
70{
71 if(fBuffer) delete[] fBuffer;
72}
73
74//---------------------------------------------------------------------------
75
76void DelphesHepMCReader::SetInputFile(FILE *inputFile)
77{
78 fInputFile = inputFile;
79}
80
81//---------------------------------------------------------------------------
82
83void DelphesHepMCReader::Clear()
84{
85 fStateSize = 0;
86 fState.clear();
87 fWeightSize = 0;
88 fWeight.clear();
89 fMomentumCoefficient = 1.0;
90 fPositionCoefficient = 1.0;
91 fVertexCounter = -1;
92 fInCounter = -1;
93 fOutCounter = -1;
94 fMotherMap.clear();
95 fDaughterMap.clear();
96 fParticleCounter = 0;
97}
98
99//---------------------------------------------------------------------------
100
101bool DelphesHepMCReader::EventReady()
102{
103 return (fVertexCounter == 0) && (fInCounter == 0) && (fOutCounter == 0);
104}
105
106//---------------------------------------------------------------------------
107
108bool DelphesHepMCReader::ReadBlock(DelphesFactory *factory,
109 TObjArray *allParticleOutputArray,
110 TObjArray *stableParticleOutputArray,
111 TObjArray *partonOutputArray)
112{
113 map< int, pair< int, int > >::iterator itMotherMap;
114 map< int, pair< int, int > >::iterator itDaughterMap;
115 char key, momentumUnit[4], positionUnit[3];
116 int i, rc, state;
117 double weight;
118
119 if(!fgets(fBuffer, kBufferSize, fInputFile)) return kFALSE;
120
121 DelphesStream bufferStream(fBuffer + 1);
122
123 key = fBuffer[0];
124
125 if(key == 'E')
126 {
127 Clear();
128
129 rc = bufferStream.ReadInt(fEventNumber)
130 && bufferStream.ReadInt(fMPI)
131 && bufferStream.ReadDbl(fScale)
132 && bufferStream.ReadDbl(fAlphaQCD)
133 && bufferStream.ReadDbl(fAlphaQED)
134 && bufferStream.ReadInt(fProcessID)
135 && bufferStream.ReadInt(fSignalCode)
136 && bufferStream.ReadInt(fVertexCounter)
137 && bufferStream.ReadInt(fBeamCode[0])
138 && bufferStream.ReadInt(fBeamCode[1])
139 && bufferStream.ReadInt(fStateSize);
140
141 if(!rc)
142 {
143 cerr << "** ERROR: " << "invalid event format" << endl;
144 return kFALSE;
145 }
146
147 for(i = 0; i < fStateSize; ++i)
148 {
149 rc = rc && bufferStream.ReadInt(state);
150 fState.push_back(state);
151 }
152
153 rc = rc && bufferStream.ReadInt(fWeightSize);
154
155 if(!rc)
156 {
157 cerr << "** ERROR: " << "invalid event format" << endl;
158 return kFALSE;
159 }
160
161 for(i = 0; i < fWeightSize; ++i)
162 {
163 rc = rc && bufferStream.ReadDbl(weight);
164 fWeight.push_back(weight);
165 }
166
167 if(!rc)
168 {
169 cerr << "** ERROR: " << "invalid event format" << endl;
170 return kFALSE;
171 }
172 }
173 else if(key == 'U')
174 {
175 rc = sscanf(fBuffer + 1, "%3s %2s", momentumUnit, positionUnit);
176
177 if(rc != 2)
178 {
179 cerr << "** ERROR: " << "invalid units format" << endl;
180 return kFALSE;
181 }
182
183 if(strncmp(momentumUnit, "GEV", 3) == 0)
184 {
185 fMomentumCoefficient = 1.0;
186 }
187 else if(strncmp(momentumUnit, "MEV", 3) == 0)
188 {
189 fMomentumCoefficient = 0.001;
190 }
191
192 if(strncmp(positionUnit, "MM", 3) == 0)
193 {
194 fPositionCoefficient = 1.0;
195 }
196 else if(strncmp(positionUnit, "CM", 3) == 0)
197 {
198 fPositionCoefficient = 10.0;
199 }
200 }
201 else if(key == 'F')
202 {
203 rc = bufferStream.ReadInt(fID1)
204 && bufferStream.ReadInt(fID2)
205 && bufferStream.ReadDbl(fX1)
206 && bufferStream.ReadDbl(fX2)
207 && bufferStream.ReadDbl(fScalePDF)
208 && bufferStream.ReadDbl(fPDF1)
209 && bufferStream.ReadDbl(fPDF2);
210
211 if(!rc)
212 {
213 cerr << "** ERROR: " << "invalid PDF format" << endl;
214 return kFALSE;
215 }
216 }
217 else if(key == 'V' && fVertexCounter > 0)
218 {
219 rc = bufferStream.ReadInt(fOutVertexCode)
220 && bufferStream.ReadInt(fVertexID)
221 && bufferStream.ReadDbl(fX)
222 && bufferStream.ReadDbl(fY)
223 && bufferStream.ReadDbl(fZ)
224 && bufferStream.ReadDbl(fT)
225 && bufferStream.ReadInt(fInCounter)
226 && bufferStream.ReadInt(fOutCounter);
227
228 if(!rc)
229 {
230 cerr << "** ERROR: " << "invalid vertex format" << endl;
231 return kFALSE;
232 }
233 --fVertexCounter;
234 }
235 else if(key == 'P' && fOutCounter > 0)
236 {
237 rc = bufferStream.ReadInt(fParticleCode)
238 && bufferStream.ReadInt(fPID)
239 && bufferStream.ReadDbl(fPx)
240 && bufferStream.ReadDbl(fPy)
241 && bufferStream.ReadDbl(fPz)
242 && bufferStream.ReadDbl(fE)
243 && bufferStream.ReadDbl(fMass)
244 && bufferStream.ReadInt(fStatus)
245 && bufferStream.ReadDbl(fTheta)
246 && bufferStream.ReadDbl(fPhi)
247 && bufferStream.ReadInt(fInVertexCode);
248
249 if(!rc)
250 {
251 cerr << "** ERROR: " << "invalid particle format" << endl;
252 return kFALSE;
253 }
254
255 if(fInVertexCode < 0)
256 {
257 itMotherMap = fMotherMap.find(fInVertexCode);
258 if(itMotherMap == fMotherMap.end())
259 {
260 fMotherMap[fInVertexCode] = make_pair(fParticleCounter, -1);
261 }
262 else
263 {
264 itMotherMap->second.second = fParticleCounter;
265 }
266 }
267
268 if(fInCounter <= 0)
269 {
270 itDaughterMap = fDaughterMap.find(fOutVertexCode);
271 if(itDaughterMap == fDaughterMap.end())
272 {
273 fDaughterMap[fOutVertexCode] = make_pair(fParticleCounter, fParticleCounter);
274 }
275 else
276 {
277 itDaughterMap->second.second = fParticleCounter;
278 }
279 }
280
281 AnalyzeParticle(factory, allParticleOutputArray,
282 stableParticleOutputArray, partonOutputArray);
283
284 if(fInCounter > 0)
285 {
286 --fInCounter;
287 }
288 else
289 {
290 --fOutCounter;
291 }
292
293 ++fParticleCounter;
294 }
295
296 if(EventReady())
297 {
298 FinalizeParticles(allParticleOutputArray);
299 }
300
301 return kTRUE;
302}
303
304//---------------------------------------------------------------------------
305
306void DelphesHepMCReader::AnalyzeEvent(ExRootTreeBranch *branch, long long eventNumber,
307 TStopwatch *readStopWatch, TStopwatch *procStopWatch)
308{
309 HepMCEvent *element;
310
311 element = static_cast<HepMCEvent *>(branch->NewEntry());
312 element->Number = fEventNumber;
313
314 element->ProcessID = fProcessID;
315 element->MPI = fMPI;
316 element->Weight = fWeight.size() > 0 ? fWeight[0] : 1.0;
317 element->Scale = fScale;
318 element->AlphaQED = fAlphaQED;
319 element->AlphaQCD = fAlphaQCD;
320
321 element->ID1 = fID1;
322 element->ID2 = fID2;
323 element->X1 = fX1;
324 element->X2 = fX2;
325 element->ScalePDF = fScalePDF;
326 element->PDF1 = fPDF1;
327 element->PDF2 = fPDF2;
328
329 element->ReadTime = readStopWatch->RealTime();
330 element->ProcTime = procStopWatch->RealTime();
331}
332
333//---------------------------------------------------------------------------
334
335void DelphesHepMCReader::AnalyzeParticle(DelphesFactory *factory,
336 TObjArray *allParticleOutputArray,
337 TObjArray *stableParticleOutputArray,
338 TObjArray *partonOutputArray)
339{
340 Candidate *candidate;
341 TParticlePDG *pdgParticle;
342 int pdgCode;
343
344 candidate = factory->NewCandidate();
345
346 candidate->PID = fPID;
347 pdgCode = TMath::Abs(candidate->PID);
348
349 candidate->Status = fStatus;
350
351 pdgParticle = fPDG->GetParticle(fPID);
352 candidate->Charge = pdgParticle ? int(pdgParticle->Charge()/3.0) : -999;
353 candidate->Mass = fMass;
354
355 candidate->Momentum.SetPxPyPzE(fPx, fPy, fPz, fE);
356 if(fMomentumCoefficient != 1.0)
357 {
358 candidate->Momentum *= fMomentumCoefficient;
359 }
360
361 candidate->M2 = 1;
362 candidate->D2 = 1;
363 if(fInCounter > 0)
364 {
365 candidate->M1 = 1;
366 candidate->Position.SetXYZT(0.0, 0.0, 0.0, 0.0);
367 }
368 else
369 {
370 candidate->M1 = fOutVertexCode;
371 candidate->Position.SetXYZT(fX, fY, fZ, fT);
372 if(fPositionCoefficient != 1.0)
373 {
374 candidate->Position *= fPositionCoefficient;
375 }
376 }
377 if(fInVertexCode < 0)
378 {
379 candidate->D1 = fInVertexCode;
380 }
381 else
382 {
383 candidate->D1 = 1;
384 }
385
386 allParticleOutputArray->Add(candidate);
387
388 if(!pdgParticle) return;
389
390 if(fStatus == 1)
391 {
392 stableParticleOutputArray->Add(candidate);
393 }
394 else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
395 {
396 partonOutputArray->Add(candidate);
397 }
398}
399
400//---------------------------------------------------------------------------
401
402void DelphesHepMCReader::FinalizeParticles(TObjArray *allParticleOutputArray)
403{
404 Candidate *candidate;
405 map< int, pair< int, int > >::iterator itMotherMap;
406 map< int, pair< int, int > >::iterator itDaughterMap;
407 int i;
408
409 for(i = 0; i < allParticleOutputArray->GetEntriesFast(); ++i)
410 {
411 candidate = static_cast<Candidate *>(allParticleOutputArray->At(i));
412
413 if(candidate->M1 > 0)
414 {
415 candidate->M1 = -1;
416 candidate->M2 = -1;
417 }
418 else
419 {
420 itMotherMap = fMotherMap.find(candidate->M1);
421 if(itMotherMap == fMotherMap.end())
422 {
423 candidate->M1 = -1;
424 candidate->M2 = -1;
425 }
426 else
427 {
428 candidate->M1 = itMotherMap->second.first;
429 candidate->M2 = itMotherMap->second.second;
430 }
431 }
432 if(candidate->D1 > 0)
433 {
434 candidate->D1 = -1;
435 candidate->D2 = -1;
436 }
437 else
438 {
439 itDaughterMap = fDaughterMap.find(candidate->D1);
440 if(itDaughterMap == fDaughterMap.end())
441 {
442 candidate->D1 = -1;
443 candidate->D2 = -1;
444 }
445 else
446 {
447 candidate->D1 = itDaughterMap->second.first;
448 candidate->D2 = itDaughterMap->second.second;
449 }
450 }
451 }
452}
453
454//---------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.