Fork me on GitHub

source: git/classes/DelphesHepMCReader.cc @ cf22deb

ImprovedOutputFileTimingdual_readoutllp
Last change on this file since cf22deb was cf22deb, checked in by pavel <pavel@…>, 7 years ago

read momentum and position units from HepMC file

  • Property mode set to 100644
File size: 10.0 KB
Line 
1
2/** \class DelphesHepMCReader
3 *
4 *  Reads HepMC file
5 *
6 *
7 *  $Date$
8 *  $Revision$
9 *
10 *
11 *  \author P. Demin - UCL, Louvain-la-Neuve
12 *
13 */
14
15#include "classes/DelphesHepMCReader.h"
16
17#include <stdexcept>
18#include <iostream>
19#include <sstream>
20
21#include <map>
22#include <vector>
23
24#include <stdio.h>
25
26#include "TObjArray.h"
27#include "TStopwatch.h"
28#include "TDatabasePDG.h"
29#include "TParticlePDG.h"
30#include "TLorentzVector.h"
31
32#include "classes/DelphesClasses.h"
33#include "classes/DelphesFactory.h"
34#include "classes/DelphesStream.h"
35
36#include "ExRootAnalysis/ExRootTreeBranch.h"
37
38using namespace std;
39
40static const int kBufferSize  = 1024;
41
42//---------------------------------------------------------------------------
43
44DelphesHepMCReader::DelphesHepMCReader() :
45  fInputFile(0), fBuffer(0), fPDG(0),
46  fVertexCounter(-1), fInCounter(-1), fOutCounter(-1),
47  fParticleCounter(0)
48{
49  fBuffer = new char[kBufferSize];
50
51  fPDG = TDatabasePDG::Instance();
52}
53
54//---------------------------------------------------------------------------
55
56DelphesHepMCReader::~DelphesHepMCReader()
57{
58  if(fBuffer) delete[] fBuffer;
59}
60
61//---------------------------------------------------------------------------
62
63void DelphesHepMCReader::SetInputFile(FILE *inputFile)
64{
65  fInputFile = inputFile;
66}
67
68//---------------------------------------------------------------------------
69
70void DelphesHepMCReader::Clear()
71{
72  fStateSize = 0;
73  fState.clear();
74  fWeightSize = 0;
75  fWeight.clear();
76  fMomentumCoefficient = 1.0;
77  fPositionCoefficient = 1.0;
78  fVertexCounter = -1;
79  fInCounter = -1;
80  fOutCounter = -1;
81  fMotherMap.clear();
82  fDaughterMap.clear();
83  fParticleCounter = 0;
84}
85
86//---------------------------------------------------------------------------
87
88bool DelphesHepMCReader::EventReady()
89{
90  return (fVertexCounter == 0) && (fInCounter == 0) && (fOutCounter == 0);
91}
92
93//---------------------------------------------------------------------------
94
95bool DelphesHepMCReader::ReadBlock(DelphesFactory *factory,
96  TObjArray *allParticleOutputArray,
97  TObjArray *stableParticleOutputArray,
98  TObjArray *partonOutputArray)
99{
100  map< int, pair< int, int > >::iterator itMotherMap;
101  map< int, pair< int, int > >::iterator itDaughterMap;
102  char key, momentumUnit[4], positionUnit[3];
103  int i, rc, state;
104  double weight;
105
106  if(!fgets(fBuffer, kBufferSize, fInputFile)) return kFALSE;
107
108  DelphesStream bufferStream(fBuffer + 1);
109
110  key = fBuffer[0];
111
112  if(key == 'E')
113  {
114    Clear();
115
116    rc = bufferStream.ReadInt(fEventNumber)
117      && bufferStream.ReadInt(fMPI)
118      && bufferStream.ReadDbl(fScale)
119      && bufferStream.ReadDbl(fAlphaQCD)
120      && bufferStream.ReadDbl(fAlphaQED)
121      && bufferStream.ReadInt(fProcessID)
122      && bufferStream.ReadInt(fSignalCode)
123      && bufferStream.ReadInt(fVertexCounter)
124      && bufferStream.ReadInt(fBeamCode[0])
125      && bufferStream.ReadInt(fBeamCode[1])
126      && bufferStream.ReadInt(fStateSize);
127
128    if(!rc)
129    {
130      cerr << "** ERROR: " << "invalid event format" << endl;
131      return kFALSE;
132    }
133
134    for(i = 0; i < fStateSize; ++i)
135    {
136      rc = rc && bufferStream.ReadInt(state);
137      fState.push_back(state);
138    }
139
140    rc = rc && bufferStream.ReadInt(fWeightSize);
141
142    if(!rc)
143    {
144      cerr << "** ERROR: " << "invalid event format" << endl;
145      return kFALSE;
146    }
147
148    for(i = 0; i < fWeightSize; ++i)
149    {
150      rc = rc && bufferStream.ReadDbl(weight);
151      fWeight.push_back(weight);
152    }
153
154    if(!rc)
155    {
156      cerr << "** ERROR: " << "invalid event format" << endl;
157      return kFALSE;
158    }
159  }
160  else if(key == 'U')
161  {
162    rc = sscanf(fBuffer + 1, "%3s %2s", momentumUnit, positionUnit);
163
164    if(rc != 2)
165    {
166      cerr << "** ERROR: " << "invalid units format" << endl;
167      return kFALSE;
168    }
169
170    if(strncmp(momentumUnit, "GEV", 3) == 0)
171    {
172      fMomentumCoefficient = 1.0;
173    }
174    else if(strncmp(momentumUnit, "MEV", 3) == 0)
175    {
176      fMomentumCoefficient = 0.001;
177    }
178   
179    if(strncmp(positionUnit, "MM", 3) == 0)
180    {
181      fPositionCoefficient = 1.0;
182    }
183    else if(strncmp(positionUnit, "CM", 3) == 0)
184    {
185      fPositionCoefficient = 10.0;
186    }
187  }
188  else if(key == 'F')
189  {
190    rc = bufferStream.ReadInt(fID1)
191      && bufferStream.ReadInt(fID2)
192      && bufferStream.ReadDbl(fX1)
193      && bufferStream.ReadDbl(fX2)
194      && bufferStream.ReadDbl(fScalePDF)
195      && bufferStream.ReadDbl(fPDF1)
196      && bufferStream.ReadDbl(fPDF2);
197
198    if(!rc)
199    {
200      cerr << "** ERROR: " << "invalid PDF format" << endl;
201      return kFALSE;
202    }
203  }
204  else if(key == 'V' && fVertexCounter > 0)
205  {
206    rc = bufferStream.ReadInt(fOutVertexCode)
207      && bufferStream.ReadInt(fVertexID)
208      && bufferStream.ReadDbl(fX)
209      && bufferStream.ReadDbl(fY)
210      && bufferStream.ReadDbl(fZ)
211      && bufferStream.ReadDbl(fT)
212      && bufferStream.ReadInt(fInCounter)
213      && bufferStream.ReadInt(fOutCounter);
214
215    if(!rc)
216    {
217      cerr << "** ERROR: " << "invalid vertex format" << endl;
218      return kFALSE;
219    }
220    --fVertexCounter;
221  }
222  else if(key == 'P' && fOutCounter > 0)
223  {
224    rc = bufferStream.ReadInt(fParticleCode)
225      && bufferStream.ReadInt(fPID)
226      && bufferStream.ReadDbl(fPx)
227      && bufferStream.ReadDbl(fPy)
228      && bufferStream.ReadDbl(fPz)
229      && bufferStream.ReadDbl(fE)
230      && bufferStream.ReadDbl(fMass)
231      && bufferStream.ReadInt(fStatus)
232      && bufferStream.ReadDbl(fTheta)
233      && bufferStream.ReadDbl(fPhi)
234      && bufferStream.ReadInt(fInVertexCode);
235
236    if(!rc)
237    {
238      cerr << "** ERROR: " << "invalid particle format" << endl;
239      return kFALSE;
240    }
241
242    if(fInVertexCode < 0)
243    {
244      itMotherMap = fMotherMap.find(fInVertexCode);
245      if(itMotherMap == fMotherMap.end())
246      {
247        fMotherMap[fInVertexCode] = make_pair(fParticleCounter, -1);
248      }
249      else
250      {
251        itMotherMap->second.second = fParticleCounter;
252      }
253    }
254
255    if(fInCounter <= 0)
256    {
257      itDaughterMap = fDaughterMap.find(fOutVertexCode);
258      if(itDaughterMap == fDaughterMap.end())
259      {
260        fDaughterMap[fOutVertexCode] = make_pair(fParticleCounter, fParticleCounter);
261      }
262      else
263      {
264        itDaughterMap->second.second = fParticleCounter;
265      }
266    }
267
268    AnalyzeParticle(factory, allParticleOutputArray,
269      stableParticleOutputArray, partonOutputArray);
270
271    if(fInCounter > 0)
272    {
273      --fInCounter;
274    }
275    else
276    {
277      --fOutCounter;
278    }
279
280    ++fParticleCounter;
281  }
282
283  if(EventReady())
284  {
285    FinalizeParticles(allParticleOutputArray);
286  }
287
288  return kTRUE;
289}
290
291//---------------------------------------------------------------------------
292
293void DelphesHepMCReader::AnalyzeEvent(ExRootTreeBranch *branch, long long eventNumber,
294  TStopwatch *readStopWatch, TStopwatch *procStopWatch)
295{
296  HepMCEvent *element;
297
298  element = static_cast<HepMCEvent *>(branch->NewEntry());
299  element->Number = fEventNumber;
300
301  element->ProcessID = fProcessID;
302  element->MPI = fMPI;
303  element->Weight = fWeight.size() > 0 ? fWeight[0] : 1.0;
304  element->Scale = fScale;
305  element->AlphaQED = fAlphaQED;
306  element->AlphaQCD = fAlphaQCD;
307
308  element->ID1 = fID1;
309  element->ID2 = fID2;
310  element->X1 = fX1;
311  element->X2 = fX2;
312  element->ScalePDF = fScalePDF;
313  element->PDF1 = fPDF1;
314  element->PDF2 = fPDF2;
315
316  element->ReadTime = readStopWatch->RealTime();
317  element->ProcTime = procStopWatch->RealTime();
318}
319
320//---------------------------------------------------------------------------
321
322void DelphesHepMCReader::AnalyzeParticle(DelphesFactory *factory,
323  TObjArray *allParticleOutputArray,
324  TObjArray *stableParticleOutputArray,
325  TObjArray *partonOutputArray)
326{
327  Candidate *candidate;
328  TParticlePDG *pdgParticle;
329  int pdgCode;
330
331  candidate = factory->NewCandidate();
332
333  candidate->PID = fPID;
334  pdgCode = TMath::Abs(candidate->PID);
335
336  candidate->Status = fStatus;
337
338  pdgParticle = fPDG->GetParticle(fPID);
339  candidate->Charge = pdgParticle ? int(pdgParticle->Charge()/3.0) : -999;
340  candidate->Mass = pdgParticle ? pdgParticle->Mass() : -999.9;
341
342  candidate->Momentum.SetPxPyPzE(fPx, fPy, fPz, fE);
343  if(fMomentumCoefficient != 1.0)
344  {
345    candidate->Momentum *= fMomentumCoefficient;
346  }
347
348  candidate->M2 = 1;
349  candidate->D2 = 1;
350  if(fInCounter > 0)
351  {
352    candidate->M1 = 1;
353    candidate->Position.SetXYZT(0.0, 0.0, 0.0, 0.0);
354  }
355  else
356  {
357    candidate->M1 = fOutVertexCode;
358    candidate->Position.SetXYZT(fX, fY, fZ, fT);
359    if(fPositionCoefficient != 1.0)
360    {
361      candidate->Position *= fPositionCoefficient;
362    }
363  }
364  if(fInVertexCode < 0)
365  {
366    candidate->D1 = fInVertexCode;
367  }
368  else
369  {
370    candidate->D1 = 1;
371  }
372
373  allParticleOutputArray->Add(candidate);
374
375  if(!pdgParticle) return;
376
377  if(fStatus == 1 && pdgParticle->Stable())
378  {
379    stableParticleOutputArray->Add(candidate);
380  }
381  else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
382  {
383    partonOutputArray->Add(candidate);
384  }
385}
386
387//---------------------------------------------------------------------------
388
389void DelphesHepMCReader::FinalizeParticles(TObjArray *allParticleOutputArray)
390{
391  Candidate *candidate;
392  map< int, pair< int, int > >::iterator itMotherMap;
393  map< int, pair< int, int > >::iterator itDaughterMap;
394  int i;
395
396  for(i = 0; i < allParticleOutputArray->GetEntriesFast(); ++i)
397  {
398    candidate = static_cast<Candidate *>(allParticleOutputArray->At(i));
399
400    if(candidate->M1 > 0)
401    {
402      candidate->M1 = -1;
403      candidate->M2 = -1;
404    }
405    else
406    {
407      itMotherMap = fMotherMap.find(candidate->M1);
408      if(itMotherMap == fMotherMap.end())
409      {
410        candidate->M1 = -1;
411        candidate->M2 = -1;
412      }
413      else
414      {
415        candidate->M1 = itMotherMap->second.first;
416        candidate->M2 = itMotherMap->second.second;
417      }
418    }
419    if(candidate->D1 > 0)
420    {
421      candidate->D1 = -1;
422      candidate->D2 = -1;
423    }
424    else
425    {
426      itDaughterMap = fDaughterMap.find(candidate->D1);
427      if(itDaughterMap == fDaughterMap.end())
428      {
429        candidate->D1 = -1;
430        candidate->D2 = -1;
431      }
432      else
433      {
434        candidate->D1 = itDaughterMap->second.first;
435        candidate->D2 = itDaughterMap->second.second;
436      }
437    }
438  }
439}
440
441//---------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.