Fork me on GitHub

source: git/classes/DelphesHepMCReader.cc @ 59abd43

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

add weight for HepMCEvent

  • Property mode set to 100644
File size: 9.2 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  fVertexCounter = -1;
77  fInCounter = -1;
78  fOutCounter = -1;
79  fMotherMap.clear();
80  fDaughterMap.clear();
81  fParticleCounter = 0;
82}
83
84//---------------------------------------------------------------------------
85
86bool DelphesHepMCReader::EventReady()
87{
88  return (fVertexCounter == 0) && (fInCounter == 0) && (fOutCounter == 0);
89}
90
91//---------------------------------------------------------------------------
92
93bool DelphesHepMCReader::ReadBlock(DelphesFactory *factory,
94  TObjArray *allParticleOutputArray,
95  TObjArray *stableParticleOutputArray,
96  TObjArray *partonOutputArray)
97{
98  map< int, pair< int, int > >::iterator itMotherMap;
99  map< int, pair< int, int > >::iterator itDaughterMap;
100  char key;
101  int i, rc, state;
102  double weight;
103
104  if(!fgets(fBuffer, kBufferSize, fInputFile)) return kFALSE;
105
106  DelphesStream bufferStream(fBuffer + 1);
107
108  key = fBuffer[0];
109
110  if(key == 'E')
111  {
112    Clear();
113
114    rc = bufferStream.ReadInt(fEventNumber)
115      && bufferStream.ReadInt(fMPI)
116      && bufferStream.ReadDbl(fScale)
117      && bufferStream.ReadDbl(fAlphaQCD)
118      && bufferStream.ReadDbl(fAlphaQED)
119      && bufferStream.ReadInt(fProcessID)
120      && bufferStream.ReadInt(fSignalCode)
121      && bufferStream.ReadInt(fVertexCounter)
122      && bufferStream.ReadInt(fBeamCode[0])
123      && bufferStream.ReadInt(fBeamCode[1])
124      && bufferStream.ReadInt(fStateSize);
125
126    if(!rc)
127    {
128      cerr << "** ERROR: " << "invalid event format" << endl;
129      return kFALSE;
130    }
131
132    for(i = 0; i < fStateSize; ++i)
133    {
134      rc = rc && bufferStream.ReadInt(state);
135      fState.push_back(state);
136    }
137
138    rc = rc && bufferStream.ReadInt(fWeightSize);
139
140    if(!rc)
141    {
142      cerr << "** ERROR: " << "invalid event format" << endl;
143      return kFALSE;
144    }
145
146    for(i = 0; i < fWeightSize; ++i)
147    {
148      rc = rc && bufferStream.ReadDbl(weight);
149      fWeight.push_back(weight);
150    }
151
152    if(!rc)
153    {
154      cerr << "** ERROR: " << "invalid event format" << endl;
155      return kFALSE;
156    }
157  }
158  else if(key == 'F')
159  {
160    rc = bufferStream.ReadInt(fID1)
161      && bufferStream.ReadInt(fID2)
162      && bufferStream.ReadDbl(fX1)
163      && bufferStream.ReadDbl(fX2)
164      && bufferStream.ReadDbl(fScalePDF)
165      && bufferStream.ReadDbl(fPDF1)
166      && bufferStream.ReadDbl(fPDF2);
167
168    if(!rc)
169    {
170      cerr << "** ERROR: " << "invalid PDF format" << endl;
171      return kFALSE;
172    }
173  }
174  else if(key == 'V' && fVertexCounter > 0)
175  {
176    rc = bufferStream.ReadInt(fOutVertexCode)
177      && bufferStream.ReadInt(fVertexID)
178      && bufferStream.ReadDbl(fX)
179      && bufferStream.ReadDbl(fY)
180      && bufferStream.ReadDbl(fZ)
181      && bufferStream.ReadDbl(fT)
182      && bufferStream.ReadInt(fInCounter)
183      && bufferStream.ReadInt(fOutCounter);
184
185    if(!rc)
186    {
187      cerr << "** ERROR: " << "invalid vertex format" << endl;
188      return kFALSE;
189    }
190    --fVertexCounter;
191  }
192  else if(key == 'P' && fOutCounter > 0)
193  {
194    rc = bufferStream.ReadInt(fParticleCode)
195      && bufferStream.ReadInt(fPID)
196      && bufferStream.ReadDbl(fPx)
197      && bufferStream.ReadDbl(fPy)
198      && bufferStream.ReadDbl(fPz)
199      && bufferStream.ReadDbl(fE)
200      && bufferStream.ReadDbl(fMass)
201      && bufferStream.ReadInt(fStatus)
202      && bufferStream.ReadDbl(fTheta)
203      && bufferStream.ReadDbl(fPhi)
204      && bufferStream.ReadInt(fInVertexCode);
205
206    if(!rc)
207    {
208      cerr << "** ERROR: " << "invalid particle format" << endl;
209      return kFALSE;
210    }
211
212    if(fInVertexCode < 0)
213    {
214      itMotherMap = fMotherMap.find(fInVertexCode);
215      if(itMotherMap == fMotherMap.end())
216      {
217        fMotherMap[fInVertexCode] = make_pair(fParticleCounter, -1);
218      }
219      else
220      {
221        itMotherMap->second.second = fParticleCounter;
222      }
223    }
224
225    if(fInCounter <= 0)
226    {
227      itDaughterMap = fDaughterMap.find(fOutVertexCode);
228      if(itDaughterMap == fDaughterMap.end())
229      {
230        fDaughterMap[fOutVertexCode] = make_pair(fParticleCounter, fParticleCounter);
231      }
232      else
233      {
234        itDaughterMap->second.second = fParticleCounter;
235      }
236    }
237
238    AnalyzeParticle(factory, allParticleOutputArray,
239      stableParticleOutputArray, partonOutputArray);
240
241    if(fInCounter > 0)
242    {
243      --fInCounter;
244    }
245    else
246    {
247      --fOutCounter;
248    }
249
250    ++fParticleCounter;
251  }
252
253  if(EventReady())
254  {
255    FinalizeParticles(allParticleOutputArray);
256  }
257
258  return kTRUE;
259}
260
261//---------------------------------------------------------------------------
262
263void DelphesHepMCReader::AnalyzeEvent(ExRootTreeBranch *branch, long long eventNumber,
264  TStopwatch *readStopWatch, TStopwatch *procStopWatch)
265{
266  HepMCEvent *element;
267
268  element = static_cast<HepMCEvent *>(branch->NewEntry());
269  element->Number = fEventNumber;
270
271  element->ProcessID = fProcessID;
272  element->MPI = fMPI;
273  element->Weight = fWeight.size() > 0 ? fWeight[0] : 1.0;
274  element->Scale = fScale;
275  element->AlphaQED = fAlphaQED;
276  element->AlphaQCD = fAlphaQCD;
277
278  element->ID1 = fID1;
279  element->ID2 = fID2;
280  element->X1 = fX1;
281  element->X2 = fX2;
282  element->ScalePDF = fScalePDF;
283  element->PDF1 = fPDF1;
284  element->PDF2 = fPDF2;
285
286  element->ReadTime = readStopWatch->RealTime();
287  element->ProcTime = procStopWatch->RealTime();
288}
289
290//---------------------------------------------------------------------------
291
292void DelphesHepMCReader::AnalyzeParticle(DelphesFactory *factory,
293  TObjArray *allParticleOutputArray,
294  TObjArray *stableParticleOutputArray,
295  TObjArray *partonOutputArray)
296{
297  Candidate *candidate;
298  TParticlePDG *pdgParticle;
299  int pdgCode;
300
301  candidate = factory->NewCandidate();
302
303  candidate->PID = fPID;
304  pdgCode = TMath::Abs(candidate->PID);
305
306  candidate->Status = fStatus;
307
308  pdgParticle = fPDG->GetParticle(fPID);
309  candidate->Charge = pdgParticle ? int(pdgParticle->Charge()/3.0) : -999;
310  candidate->Mass = pdgParticle ? pdgParticle->Mass() : -999.9;
311
312  candidate->Momentum.SetPxPyPzE(fPx, fPy, fPz, fE);
313
314  candidate->M2 = 1;
315  candidate->D2 = 1;
316  if(fInCounter > 0)
317  {
318    candidate->M1 = 1;
319    candidate->Position.SetXYZT(0.0, 0.0, 0.0, 0.0);
320  }
321  else
322  {
323    candidate->M1 = fOutVertexCode;
324    candidate->Position.SetXYZT(fX, fY, fZ, fT);
325  }
326  if(fInVertexCode < 0)
327  {
328    candidate->D1 = fInVertexCode;
329  }
330  else
331  {
332    candidate->D1 = 1;
333  }
334
335  allParticleOutputArray->Add(candidate);
336
337  if(!pdgParticle) return;
338
339  if(fStatus == 1 && pdgParticle->Stable())
340  {
341    stableParticleOutputArray->Add(candidate);
342  }
343  else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
344  {
345    partonOutputArray->Add(candidate);
346  }
347}
348
349//---------------------------------------------------------------------------
350
351void DelphesHepMCReader::FinalizeParticles(TObjArray *allParticleOutputArray)
352{
353  Candidate *candidate;
354  map< int, pair< int, int > >::iterator itMotherMap;
355  map< int, pair< int, int > >::iterator itDaughterMap;
356  int i;
357
358  for(i = 0; i < allParticleOutputArray->GetEntriesFast(); ++i)
359  {
360    candidate = static_cast<Candidate *>(allParticleOutputArray->At(i));
361
362    if(candidate->M1 > 0)
363    {
364      candidate->M1 = -1;
365      candidate->M2 = -1;
366    }
367    else
368    {
369      itMotherMap = fMotherMap.find(candidate->M1);
370      if(itMotherMap == fMotherMap.end())
371      {
372        candidate->M1 = -1;
373        candidate->M2 = -1;
374      }
375      else
376      {
377        candidate->M1 = itMotherMap->second.first;
378        candidate->M2 = itMotherMap->second.second;
379      }
380    }
381    if(candidate->D1 > 0)
382    {
383      candidate->D1 = -1;
384      candidate->D2 = -1;
385    }
386    else
387    {
388      itDaughterMap = fDaughterMap.find(candidate->D1);
389      if(itDaughterMap == fDaughterMap.end())
390      {
391        candidate->D1 = -1;
392        candidate->D2 = -1;
393      }
394      else
395      {
396        candidate->D1 = itDaughterMap->second.first;
397        candidate->D2 = itDaughterMap->second.second;
398      }
399    }
400  }
401}
402
403//---------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.