Fork me on GitHub

source: git/classes/DelphesHepMCReader.cc@ caba091

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

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