Fork me on GitHub

source: git/classes/DelphesHepMC2Reader.cc

Last change on this file was cc8716b, checked in by GitHub <noreply@…>, 2 years ago

Update to handle CMS endcap muon detector showers for long-lived particles (#103)

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