Fork me on GitHub

source: git/classes/DelphesHepMCReader.cc

Last change on this file was 77e9ae1, checked in by Pavel Demin <pavel-demin@…>, 9 months ago

set Standard to Cpp03 in .clang-format

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