Fork me on GitHub

source: git/classes/DelphesHepMCReader.cc@ f50b2aa

ImprovedOutputFile Timing dual_readout llp
Last change on this file since f50b2aa was edeb0f0, checked in by Michele Selvaggi <michele.selvaggi@…>, 6 years ago

added cross section and error in hepmc event class and reader

  • Property mode set to 100644
File size: 11.3 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
202 else if(key == 'C')
203 {
204 rc = bufferStream.ReadDbl(fCrossSection)
205 && bufferStream.ReadDbl(fCrossSectionError);
206 }
207
208 else if(key == 'F')
209 {
210 rc = bufferStream.ReadInt(fID1)
211 && bufferStream.ReadInt(fID2)
212 && bufferStream.ReadDbl(fX1)
213 && bufferStream.ReadDbl(fX2)
214 && bufferStream.ReadDbl(fScalePDF)
215 && bufferStream.ReadDbl(fPDF1)
216 && bufferStream.ReadDbl(fPDF2);
217
218 if(!rc)
219 {
220 cerr << "** ERROR: " << "invalid PDF format" << endl;
221 return kFALSE;
222 }
223 }
224 else if(key == 'V' && fVertexCounter > 0)
225 {
226 rc = bufferStream.ReadInt(fOutVertexCode)
227 && bufferStream.ReadInt(fVertexID)
228 && bufferStream.ReadDbl(fX)
229 && bufferStream.ReadDbl(fY)
230 && bufferStream.ReadDbl(fZ)
231 && bufferStream.ReadDbl(fT)
232 && bufferStream.ReadInt(fInCounter)
233 && bufferStream.ReadInt(fOutCounter);
234
235 if(!rc)
236 {
237 cerr << "** ERROR: " << "invalid vertex format" << endl;
238 return kFALSE;
239 }
240 --fVertexCounter;
241 }
242 else if(key == 'P' && fOutCounter > 0)
243 {
244 rc = bufferStream.ReadInt(fParticleCode)
245 && bufferStream.ReadInt(fPID)
246 && bufferStream.ReadDbl(fPx)
247 && bufferStream.ReadDbl(fPy)
248 && bufferStream.ReadDbl(fPz)
249 && bufferStream.ReadDbl(fE)
250 && bufferStream.ReadDbl(fMass)
251 && bufferStream.ReadInt(fStatus)
252 && bufferStream.ReadDbl(fTheta)
253 && bufferStream.ReadDbl(fPhi)
254 && bufferStream.ReadInt(fInVertexCode);
255
256 if(!rc)
257 {
258 cerr << "** ERROR: " << "invalid particle format" << endl;
259 return kFALSE;
260 }
261
262 if(fInVertexCode < 0)
263 {
264 itMotherMap = fMotherMap.find(fInVertexCode);
265 if(itMotherMap == fMotherMap.end())
266 {
267 fMotherMap[fInVertexCode] = make_pair(fParticleCounter, -1);
268 }
269 else
270 {
271 itMotherMap->second.second = fParticleCounter;
272 }
273 }
274
275 if(fInCounter <= 0)
276 {
277 itDaughterMap = fDaughterMap.find(fOutVertexCode);
278 if(itDaughterMap == fDaughterMap.end())
279 {
280 fDaughterMap[fOutVertexCode] = make_pair(fParticleCounter, fParticleCounter);
281 }
282 else
283 {
284 itDaughterMap->second.second = fParticleCounter;
285 }
286 }
287
288 AnalyzeParticle(factory, allParticleOutputArray,
289 stableParticleOutputArray, partonOutputArray);
290
291 if(fInCounter > 0)
292 {
293 --fInCounter;
294 }
295 else
296 {
297 --fOutCounter;
298 }
299
300 ++fParticleCounter;
301 }
302
303 if(EventReady())
304 {
305 FinalizeParticles(allParticleOutputArray);
306 }
307
308 return kTRUE;
309}
310
311//---------------------------------------------------------------------------
312
313void DelphesHepMCReader::AnalyzeEvent(ExRootTreeBranch *branch, long long eventNumber,
314 TStopwatch *readStopWatch, TStopwatch *procStopWatch)
315{
316 HepMCEvent *element;
317
318 element = static_cast<HepMCEvent *>(branch->NewEntry());
319 element->Number = fEventNumber;
320
321 element->ProcessID = fProcessID;
322 element->MPI = fMPI;
323 element->Weight = fWeight.size() > 0 ? fWeight[0] : 1.0;
324 element->CrossSection = fCrossSection;
325 element->CrossSectionError = fCrossSectionError;
326 element->Scale = fScale;
327 element->AlphaQED = fAlphaQED;
328 element->AlphaQCD = fAlphaQCD;
329
330 element->ID1 = fID1;
331 element->ID2 = fID2;
332 element->X1 = fX1;
333 element->X2 = fX2;
334 element->ScalePDF = fScalePDF;
335 element->PDF1 = fPDF1;
336 element->PDF2 = fPDF2;
337
338 element->ReadTime = readStopWatch->RealTime();
339 element->ProcTime = procStopWatch->RealTime();
340}
341
342//---------------------------------------------------------------------------
343
344void DelphesHepMCReader::AnalyzeWeight(ExRootTreeBranch *branch)
345{
346 Weight *element;
347 vector< double >::const_iterator itWeight;
348
349 for(itWeight = fWeight.begin(); itWeight != fWeight.end(); ++itWeight)
350 {
351 element = static_cast<Weight *>(branch->NewEntry());
352
353 element->Weight = *itWeight;
354 }
355}
356
357//---------------------------------------------------------------------------
358
359void DelphesHepMCReader::AnalyzeParticle(DelphesFactory *factory,
360 TObjArray *allParticleOutputArray,
361 TObjArray *stableParticleOutputArray,
362 TObjArray *partonOutputArray)
363{
364 Candidate *candidate;
365 TParticlePDG *pdgParticle;
366 int pdgCode;
367
368 candidate = factory->NewCandidate();
369
370 candidate->PID = fPID;
371 pdgCode = TMath::Abs(candidate->PID);
372
373 candidate->Status = fStatus;
374
375 pdgParticle = fPDG->GetParticle(fPID);
376 candidate->Charge = pdgParticle ? int(pdgParticle->Charge()/3.0) : -999;
377 candidate->Mass = fMass;
378
379 candidate->Momentum.SetPxPyPzE(fPx, fPy, fPz, fE);
380 if(fMomentumCoefficient != 1.0)
381 {
382 candidate->Momentum *= fMomentumCoefficient;
383 }
384
385 candidate->M2 = 1;
386 candidate->D2 = 1;
387 if(fInCounter > 0)
388 {
389 candidate->M1 = 1;
390 candidate->Position.SetXYZT(0.0, 0.0, 0.0, 0.0);
391 }
392 else
393 {
394 candidate->M1 = fOutVertexCode;
395 candidate->Position.SetXYZT(fX, fY, fZ, fT);
396 if(fPositionCoefficient != 1.0)
397 {
398 candidate->Position *= fPositionCoefficient;
399 }
400 }
401 if(fInVertexCode < 0)
402 {
403 candidate->D1 = fInVertexCode;
404 }
405 else
406 {
407 candidate->D1 = 1;
408 }
409
410 allParticleOutputArray->Add(candidate);
411
412 if(!pdgParticle) return;
413
414 if(fStatus == 1)
415 {
416 stableParticleOutputArray->Add(candidate);
417 }
418 else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
419 {
420 partonOutputArray->Add(candidate);
421 }
422}
423
424//---------------------------------------------------------------------------
425
426void DelphesHepMCReader::FinalizeParticles(TObjArray *allParticleOutputArray)
427{
428 Candidate *candidate;
429 map< int, pair< int, int > >::iterator itMotherMap;
430 map< int, pair< int, int > >::iterator itDaughterMap;
431 int i;
432
433 for(i = 0; i < allParticleOutputArray->GetEntriesFast(); ++i)
434 {
435 candidate = static_cast<Candidate *>(allParticleOutputArray->At(i));
436
437 if(candidate->M1 > 0)
438 {
439 candidate->M1 = -1;
440 candidate->M2 = -1;
441 }
442 else
443 {
444 itMotherMap = fMotherMap.find(candidate->M1);
445 if(itMotherMap == fMotherMap.end())
446 {
447 candidate->M1 = -1;
448 candidate->M2 = -1;
449 }
450 else
451 {
452 candidate->M1 = itMotherMap->second.first;
453 candidate->M2 = itMotherMap->second.second;
454 }
455 }
456 if(candidate->D1 > 0)
457 {
458 candidate->D1 = -1;
459 candidate->D2 = -1;
460 }
461 else
462 {
463 itDaughterMap = fDaughterMap.find(candidate->D1);
464 if(itDaughterMap == fDaughterMap.end())
465 {
466 candidate->D1 = -1;
467 candidate->D2 = -1;
468 }
469 else
470 {
471 candidate->D1 = itDaughterMap->second.first;
472 candidate->D2 = itDaughterMap->second.second;
473 }
474 }
475 }
476}
477
478//---------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.