Fork me on GitHub

source: git/classes/DelphesHepMC3Reader.cc@ 3c59f98

Last change on this file since 3c59f98 was 29b722a, checked in by christinaw97 <christina.wang@…>, 3 years ago

Added DelphesCscClusterFormula
  • Property mode set to 100644
File size: 14.0 KB
Line 
1/*
2 * Delphes: a framework for fast simulation of a generic collider experiment
3 * Copyright (C) 2012-2021 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/** \class DelphesHepMC3Reader
20 *
21 * Reads HepMC file
22 *
23 * \author P. Demin - UCL, Louvain-la-Neuve
24 *
25 */
26
27#include "classes/DelphesHepMC3Reader.h"
28
29#include <iostream>
30#include <sstream>
31#include <stdexcept>
32
33#include <map>
34#include <vector>
35
36#include <stdio.h>
37
38#include "TDatabasePDG.h"
39#include "TLorentzVector.h"
40#include "TObjArray.h"
41#include "TParticlePDG.h"
42#include "TStopwatch.h"
43
44#include "classes/DelphesClasses.h"
45#include "classes/DelphesFactory.h"
46#include "classes/DelphesStream.h"
47
48#include "ExRootAnalysis/ExRootTreeBranch.h"
49
50using namespace std;
51
52static const int kBufferSize = 16384;
53
54//---------------------------------------------------------------------------
55
56DelphesHepMC3Reader::DelphesHepMC3Reader() :
57 fInputFile(0), fBuffer(0), fPDG(0),
58 fVertexCounter(-2), fParticleCounter(-1)
59{
60 fBuffer = new char[kBufferSize];
61
62 fPDG = TDatabasePDG::Instance();
63}
64
65//---------------------------------------------------------------------------
66
67DelphesHepMC3Reader::~DelphesHepMC3Reader()
68{
69 if(fBuffer) delete[] fBuffer;
70}
71
72//---------------------------------------------------------------------------
73
74void DelphesHepMC3Reader::SetInputFile(FILE *inputFile)
75{
76 fInputFile = inputFile;
77}
78
79//---------------------------------------------------------------------------
80
81void DelphesHepMC3Reader::Clear()
82{
83 fWeights.clear();
84 fMomentumCoefficient = 1.0;
85 fPositionCoefficient = 1.0;
86 fVertexCounter = -2;
87 fParticleCounter = -1;
88 fVertices.clear();
89 fParticles.clear();
90 fInVertexMap.clear();
91 fOutVertexMap.clear();
92 fMotherMap.clear();
93 fDaughterMap.clear();
94}
95
96//---------------------------------------------------------------------------
97
98bool DelphesHepMC3Reader::EventReady()
99{
100 return (fVertexCounter == -1) && (fParticleCounter == 0);
101}
102
103//---------------------------------------------------------------------------
104
105bool DelphesHepMC3Reader::ReadBlock(DelphesFactory *factory,
106 TObjArray *allParticleOutputArray,
107 TObjArray *stableParticleOutputArray,
108 TObjArray *partonOutputArray)
109{
110 map<int, pair<int, int> >::iterator itDaughterMap;
111 char key, momentumUnit[4], positionUnit[3];
112 int rc, code;
113 double weight;
114
115 if(!fgets(fBuffer, kBufferSize, fInputFile)) return kFALSE;
116
117 DelphesStream bufferStream(fBuffer + 1);
118
119 key = fBuffer[0];
120
121 if(key == 'E')
122 {
123 Clear();
124
125 rc = bufferStream.ReadInt(fEventNumber)
126 && bufferStream.ReadInt(fVertexCounter)
127 && bufferStream.ReadInt(fParticleCounter);
128
129 if(!rc)
130 {
131 cerr << "** ERROR: "
132 << "invalid event format" << endl;
133 return kFALSE;
134 }
135 }
136 else if(key == 'U')
137 {
138 rc = sscanf(fBuffer + 1, "%3s %2s", momentumUnit, positionUnit);
139
140 if(rc != 2)
141 {
142 cerr << "** ERROR: "
143 << "invalid units format" << endl;
144 return kFALSE;
145 }
146
147 if(strncmp(momentumUnit, "GEV", 3) == 0)
148 {
149 fMomentumCoefficient = 1.0;
150 }
151 else if(strncmp(momentumUnit, "MEV", 3) == 0)
152 {
153 fMomentumCoefficient = 0.001;
154 }
155
156 if(strncmp(positionUnit, "MM", 3) == 0)
157 {
158 fPositionCoefficient = 1.0;
159 }
160 else if(strncmp(positionUnit, "CM", 3) == 0)
161 {
162 fPositionCoefficient = 10.0;
163 }
164 }
165 else if(key == 'W')
166 {
167 while(bufferStream.ReadDbl(weight))
168 {
169 fWeights.push_back(weight);
170 }
171 }
172 else if(key == 'A' && bufferStream.FindStr("mpi"))
173 {
174 rc = bufferStream.ReadInt(fMPI);
175
176 if(!rc)
177 {
178 cerr << "** ERROR: "
179 << "invalid MPI format" << endl;
180 return kFALSE;
181 }
182 }
183 else if(key == 'A' && bufferStream.FindStr("signal_process_id"))
184 {
185 rc = bufferStream.ReadInt(fProcessID);
186
187 if(!rc)
188 {
189 cerr << "** ERROR: "
190 << "invalid process ID format" << endl;
191 return kFALSE;
192 }
193 }
194 else if(key == 'A' && bufferStream.FindStr("event_scale"))
195 {
196 rc = bufferStream.ReadDbl(fScale);
197
198 if(!rc)
199 {
200 cerr << "** ERROR: "
201 << "invalid event scale format" << endl;
202 return kFALSE;
203 }
204 }
205 else if(key == 'A' && bufferStream.FindStr("alphaQCD"))
206 {
207 rc = bufferStream.ReadDbl(fAlphaQCD);
208
209 if(!rc)
210 {
211 cerr << "** ERROR: "
212 << "invalid alphaQCD format" << endl;
213 return kFALSE;
214 }
215 }
216 else if(key == 'A' && bufferStream.FindStr("alphaQED"))
217 {
218 rc = bufferStream.ReadDbl(fAlphaQED);
219
220 if(!rc)
221 {
222 cerr << "** ERROR: "
223 << "invalid alphaQED format" << endl;
224 return kFALSE;
225 }
226 }
227 else if(key == 'A' && bufferStream.FindStr("GenCrossSection"))
228 {
229 rc = bufferStream.ReadDbl(fCrossSection)
230 && bufferStream.ReadDbl(fCrossSectionError);
231
232 if(!rc)
233 {
234 cerr << "** ERROR: "
235 << "invalid cross section format" << endl;
236 return kFALSE;
237 }
238 }
239 else if(key == 'A' && bufferStream.FindStr("GenPdfInfo"))
240 {
241 rc = bufferStream.ReadInt(fID1)
242 && bufferStream.ReadInt(fID2)
243 && bufferStream.ReadDbl(fX1)
244 && bufferStream.ReadDbl(fX2)
245 && bufferStream.ReadDbl(fScalePDF)
246 && bufferStream.ReadDbl(fPDF1)
247 && bufferStream.ReadDbl(fPDF2);
248
249 if(!rc)
250 {
251 cerr << "** ERROR: "
252 << "invalid PDF format" << endl;
253 return kFALSE;
254 }
255 }
256 else if(key == 'V')
257 {
258 fParticles.clear();
259
260 fX = 0.0;
261 fY = 0.0;
262 fZ = 0.0;
263 fT = 0.0;
264
265 rc = bufferStream.ReadInt(fVertexCode)
266 && bufferStream.ReadInt(fVertexStatus);
267
268 if(!rc)
269 {
270 cerr << "** ERROR: "
271 << "invalid vertex format" << endl;
272 return kFALSE;
273 }
274
275 rc = bufferStream.FindChr('[');
276
277 if(!rc)
278 {
279 cerr << "** ERROR: "
280 << "invalid vertex format" << endl;
281 return kFALSE;
282 }
283
284 while(bufferStream.ReadInt(code))
285 {
286 fParticles.push_back(code);
287 bufferStream.FindChr(',');
288 }
289
290 if(bufferStream.FindChr('@'))
291 {
292 rc = bufferStream.ReadDbl(fX)
293 && bufferStream.ReadDbl(fY)
294 && bufferStream.ReadDbl(fZ)
295 && bufferStream.ReadDbl(fT);
296
297 if(!rc)
298 {
299 cerr << "** ERROR: "
300 << "invalid vertex format" << endl;
301 return kFALSE;
302 }
303 }
304
305 AnalyzeVertex(factory, fVertexCode);
306 }
307 else if(key == 'P' && fParticleCounter > 0)
308 {
309 --fParticleCounter;
310
311 rc = bufferStream.ReadInt(fParticleCode)
312 && bufferStream.ReadInt(fOutVertexCode)
313 && bufferStream.ReadInt(fPID)
314 && bufferStream.ReadDbl(fPx)
315 && bufferStream.ReadDbl(fPy)
316 && bufferStream.ReadDbl(fPz)
317 && bufferStream.ReadDbl(fE)
318 && bufferStream.ReadDbl(fMass)
319 && bufferStream.ReadInt(fParticleStatus);
320
321 if(!rc)
322 {
323 cerr << "** ERROR: "
324 << "invalid particle format" << endl;
325 return kFALSE;
326 }
327
328 AnalyzeParticle(factory);
329 }
330
331 if(EventReady())
332 {
333 FinalizeParticles(allParticleOutputArray, stableParticleOutputArray, partonOutputArray);
334 }
335
336 return kTRUE;
337}
338
339//---------------------------------------------------------------------------
340
341void DelphesHepMC3Reader::AnalyzeEvent(ExRootTreeBranch *branch, long long eventNumber,
342 TStopwatch *readStopWatch, TStopwatch *procStopWatch)
343{
344 HepMCEvent *element;
345
346 element = static_cast<HepMCEvent *>(branch->NewEntry());
347 element->Number = fEventNumber;
348
349 element->ProcessID = fProcessID;
350 element->MPI = fMPI;
351 element->Weight = fWeights.size() > 0 ? fWeights[0] : 1.0;
352 element->CrossSection = fCrossSection;
353 element->CrossSectionError = fCrossSectionError;
354 element->Scale = fScale;
355 element->AlphaQED = fAlphaQED;
356 element->AlphaQCD = fAlphaQCD;
357
358 element->ID1 = fID1;
359 element->ID2 = fID2;
360 element->X1 = fX1;
361 element->X2 = fX2;
362 element->ScalePDF = fScalePDF;
363 element->PDF1 = fPDF1;
364 element->PDF2 = fPDF2;
365
366 element->ReadTime = readStopWatch->RealTime();
367 element->ProcTime = procStopWatch->RealTime();
368}
369
370//---------------------------------------------------------------------------
371
372void DelphesHepMC3Reader::AnalyzeWeight(ExRootTreeBranch *branch)
373{
374 Weight *element;
375 vector<double>::const_iterator itWeights;
376
377 for(itWeights = fWeights.begin(); itWeights != fWeights.end(); ++itWeights)
378 {
379 element = static_cast<Weight *>(branch->NewEntry());
380
381 element->Weight = *itWeights;
382 }
383}
384
385//---------------------------------------------------------------------------
386
387void DelphesHepMC3Reader::AnalyzeVertex(DelphesFactory *factory, int code, Candidate *candidate)
388{
389 int index;
390 TLorentzVector *position;
391 TObjArray *array;
392 vector<int>::iterator itParticle;
393 map<int, int>::iterator itVertexMap;
394
395 itVertexMap = fOutVertexMap.find(code);
396 if(itVertexMap == fOutVertexMap.end())
397 {
398 --fVertexCounter;
399
400 index = fVertices.size();
401 fOutVertexMap[code] = index;
402 if(candidate && code > 0) fInVertexMap[code] = index;
403
404 position = factory->New<TLorentzVector>();
405 array = factory->NewArray();
406 position->SetXYZT(0.0, 0.0, 0.0, 0.0);
407 fVertices.push_back(make_pair(position, array));
408 }
409 else
410 {
411 index = itVertexMap->second;
412 position = fVertices[index].first;
413 array = fVertices[index].second;
414 }
415
416 if(candidate)
417 {
418 array->Add(candidate);
419 }
420 else
421 {
422 position->SetXYZT(fX, fY, fZ, fT);
423 for(itParticle = fParticles.begin(); itParticle != fParticles.end(); ++itParticle)
424 {
425 fInVertexMap[*itParticle] = index;
426 }
427 }
428}
429
430//---------------------------------------------------------------------------
431
432void DelphesHepMC3Reader::AnalyzeParticle(DelphesFactory *factory)
433{
434 Candidate *candidate;
435
436 candidate = factory->NewCandidate();
437
438 candidate->PID = fPID;
439
440 candidate->Status = fParticleStatus;
441
442 candidate->Mass = fMass;
443
444 candidate->Momentum.SetPxPyPzE(fPx, fPy, fPz, fE);
445
446 candidate->D1 = fParticleCode;
447
448 AnalyzeVertex(factory, fOutVertexCode, candidate);
449}
450
451//---------------------------------------------------------------------------
452
453void DelphesHepMC3Reader::FinalizeParticles(TObjArray *allParticleOutputArray,
454 TObjArray *stableParticleOutputArray,
455 TObjArray *partonOutputArray)
456{
457 TLorentzVector *position;
458 TObjArray *array;
459 Candidate *candidate;
460 Candidate *candidateDaughter;
461 TParticlePDG *pdgParticle;
462 int pdgCode;
463 map<int, int >::iterator itVertexMap;
464 map<int, pair<int, int> >::iterator itMotherMap;
465 map<int, pair<int, int> >::iterator itDaughterMap;
466 int i, j, code, counter;
467
468 counter = 0;
469 for(i = 0; i < fVertices.size(); ++i)
470 {
471 position = fVertices[i].first;
472 array = fVertices[i].second;
473
474 for(j = 0; j < array->GetEntriesFast(); ++j)
475 {
476 candidate = static_cast<Candidate *>(array->At(j));
477
478 candidate->Position = *position;
479 if(fPositionCoefficient != 1.0)
480 {
481 candidate->Position *= fPositionCoefficient;
482 }
483
484 if(fMomentumCoefficient != 1.0)
485 {
486 candidate->Momentum *= fMomentumCoefficient;
487 }
488
489 candidate->M1 = i;
490
491 itDaughterMap = fDaughterMap.find(i);
492 if(itDaughterMap == fDaughterMap.end())
493 {
494 fDaughterMap[i] = make_pair(counter, counter);
495 }
496 else
497 {
498 itDaughterMap->second.second = counter;
499 }
500
501 code = candidate->D1;
502
503 itVertexMap = fInVertexMap.find(code);
504 if(itVertexMap == fInVertexMap.end())
505 {
506 candidate->D1 = -1;
507 }
508 else
509 {
510 code = itVertexMap->second;
511
512 candidate->D1 = code;
513
514 itMotherMap = fMotherMap.find(code);
515 if(itMotherMap == fMotherMap.end())
516 {
517 fMotherMap[code] = make_pair(counter, -1);
518 }
519 else
520 {
521 itMotherMap->second.second = counter;
522 }
523 }
524
525 allParticleOutputArray->Add(candidate);
526
527 ++counter;
528
529 pdgParticle = fPDG->GetParticle(candidate->PID);
530
531 candidate->Charge = pdgParticle ? int(pdgParticle->Charge() / 3.0) : -999;
532
533 if(!pdgParticle) continue;
534
535 pdgCode = TMath::Abs(candidate->PID);
536
537 if(candidate->Status == 1)
538 {
539 stableParticleOutputArray->Add(candidate);
540 }
541 else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
542 {
543 partonOutputArray->Add(candidate);
544 }
545 }
546 }
547
548 for(i = 0; i < allParticleOutputArray->GetEntriesFast(); ++i)
549 {
550 candidate = static_cast<Candidate *>(allParticleOutputArray->At(i));
551
552 itMotherMap = fMotherMap.find(candidate->M1);
553 if(itMotherMap == fMotherMap.end())
554 {
555 candidate->M1 = -1;
556 candidate->M2 = -1;
557 }
558 else
559 {
560 candidate->M1 = itMotherMap->second.first;
561 candidate->M2 = itMotherMap->second.second;
562 }
563
564 if(candidate->D1 < 0)
565 {
566 candidate->D1 = -1;
567 candidate->D2 = -1;
568 }
569 else
570 {
571 itDaughterMap = fDaughterMap.find(candidate->D1);
572 if(itDaughterMap == fDaughterMap.end())
573 {
574 candidate->D1 = -1;
575 candidate->D2 = -1;
576 const TLorentzVector &decayPosition = candidate->Position;
577 candidate->DecayPosition.SetXYZT(decayPosition.X(), decayPosition.Y(), decayPosition.Z(), decayPosition.T());// decay position
578 }
579 else
580 {
581 candidate->D1 = itDaughterMap->second.first;
582 candidate->D2 = itDaughterMap->second.second;
583 candidateDaughter = static_cast<Candidate *>(allParticleOutputArray->At(candidate->D1));
584 const TLorentzVector &decayPosition = candidateDaughter->Position;
585 candidate->DecayPosition.SetXYZT(decayPosition.X(), decayPosition.Y(), decayPosition.Z(), decayPosition.T());// decay position
586 }
587 }
588 }
589}
590
591//---------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.