Fork me on GitHub

source: git/classes/DelphesHepMCReader.cc@ b8b29f1

ImprovedOutputFile Timing dual_readout llp
Last change on this file since b8b29f1 was d203a37, checked in by Pavel Demin <pavel-demin@…>, 7 years ago

add branch Weight to DelphesHepMC.cpp
  • Property mode set to 100644
File size: 11.1 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 = 16384;
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::AnalyzeWeight(ExRootTreeBranch *branch)
336{
337 Weight *element;
338 vector< double >::const_iterator itWeight;
339
340 for(itWeight = fWeight.begin(); itWeight != fWeight.end(); ++itWeight)
341 {
342 element = static_cast<Weight *>(branch->NewEntry());
343
344 element->Weight = *itWeight;
345 }
346}
347
348//---------------------------------------------------------------------------
349
350void DelphesHepMCReader::AnalyzeParticle(DelphesFactory *factory,
351 TObjArray *allParticleOutputArray,
352 TObjArray *stableParticleOutputArray,
353 TObjArray *partonOutputArray)
354{
355 Candidate *candidate;
356 TParticlePDG *pdgParticle;
357 int pdgCode;
358
359 candidate = factory->NewCandidate();
360
361 candidate->PID = fPID;
362 pdgCode = TMath::Abs(candidate->PID);
363
364 candidate->Status = fStatus;
365
366 pdgParticle = fPDG->GetParticle(fPID);
367 candidate->Charge = pdgParticle ? int(pdgParticle->Charge()/3.0) : -999;
368 candidate->Mass = fMass;
369
370 candidate->Momentum.SetPxPyPzE(fPx, fPy, fPz, fE);
371 if(fMomentumCoefficient != 1.0)
372 {
373 candidate->Momentum *= fMomentumCoefficient;
374 }
375
376 candidate->M2 = 1;
377 candidate->D2 = 1;
378 if(fInCounter > 0)
379 {
380 candidate->M1 = 1;
381 candidate->Position.SetXYZT(0.0, 0.0, 0.0, 0.0);
382 }
383 else
384 {
385 candidate->M1 = fOutVertexCode;
386 candidate->Position.SetXYZT(fX, fY, fZ, fT);
387 if(fPositionCoefficient != 1.0)
388 {
389 candidate->Position *= fPositionCoefficient;
390 }
391 }
392 if(fInVertexCode < 0)
393 {
394 candidate->D1 = fInVertexCode;
395 }
396 else
397 {
398 candidate->D1 = 1;
399 }
400
401 allParticleOutputArray->Add(candidate);
402
403 if(!pdgParticle) return;
404
405 if(fStatus == 1)
406 {
407 stableParticleOutputArray->Add(candidate);
408 }
409 else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
410 {
411 partonOutputArray->Add(candidate);
412 }
413}
414
415//---------------------------------------------------------------------------
416
417void DelphesHepMCReader::FinalizeParticles(TObjArray *allParticleOutputArray)
418{
419 Candidate *candidate;
420 map< int, pair< int, int > >::iterator itMotherMap;
421 map< int, pair< int, int > >::iterator itDaughterMap;
422 int i;
423
424 for(i = 0; i < allParticleOutputArray->GetEntriesFast(); ++i)
425 {
426 candidate = static_cast<Candidate *>(allParticleOutputArray->At(i));
427
428 if(candidate->M1 > 0)
429 {
430 candidate->M1 = -1;
431 candidate->M2 = -1;
432 }
433 else
434 {
435 itMotherMap = fMotherMap.find(candidate->M1);
436 if(itMotherMap == fMotherMap.end())
437 {
438 candidate->M1 = -1;
439 candidate->M2 = -1;
440 }
441 else
442 {
443 candidate->M1 = itMotherMap->second.first;
444 candidate->M2 = itMotherMap->second.second;
445 }
446 }
447 if(candidate->D1 > 0)
448 {
449 candidate->D1 = -1;
450 candidate->D2 = -1;
451 }
452 else
453 {
454 itDaughterMap = fDaughterMap.find(candidate->D1);
455 if(itDaughterMap == fDaughterMap.end())
456 {
457 candidate->D1 = -1;
458 candidate->D2 = -1;
459 }
460 else
461 {
462 candidate->D1 = itDaughterMap->second.first;
463 candidate->D2 = itDaughterMap->second.second;
464 }
465 }
466 }
467}
468
469//---------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.