Fork me on GitHub

source: git/classes/DelphesHepMCReader.cc @ edeb0f0

ImprovedOutputFileTimingdual_readoutllp
Last change on this file since edeb0f0 was edeb0f0, checked in by Michele Selvaggi <michele.selvaggi@…>, 2 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.