Fork me on GitHub

source: git/classes/DelphesHepMCReader.cc @ d7d2da3

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

move branches/ModularDelphes to trunk

  • Property mode set to 100644
File size: 8.4 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
23#include <stdio.h>
24
25#include "TObjArray.h"
26#include "TStopwatch.h"
27#include "TDatabasePDG.h"
28#include "TParticlePDG.h"
29#include "TLorentzVector.h"
30
31#include "classes/DelphesClasses.h"
32#include "classes/DelphesFactory.h"
33#include "classes/DelphesStream.h"
34
35#include "ExRootAnalysis/ExRootTreeBranch.h"
36
37using namespace std;
38
39static const int kBufferSize  = 1024;
40
41//---------------------------------------------------------------------------
42
43DelphesHepMCReader::DelphesHepMCReader() :
44  fInputFile(0), fBuffer(0), fPDG(0),
45  fVertexCounter(-1), fInCounter(-1), fOutCounter(-1),
46  fParticleCounter(0)
47{
48  fBuffer = new char[kBufferSize];
49
50  fPDG = TDatabasePDG::Instance();
51}
52
53//---------------------------------------------------------------------------
54
55DelphesHepMCReader::~DelphesHepMCReader()
56{
57  if(fBuffer) delete[] fBuffer;
58}
59
60//---------------------------------------------------------------------------
61
62void DelphesHepMCReader::SetInputFile(FILE *inputFile)
63{
64  fInputFile = inputFile;
65}
66
67//---------------------------------------------------------------------------
68
69void DelphesHepMCReader::Clear()
70{
71  fVertexCounter = -1;
72  fInCounter = -1;
73  fOutCounter = -1;
74  fMotherMap.clear();
75  fDaughterMap.clear();
76  fParticleCounter = 0;
77}
78
79//---------------------------------------------------------------------------
80
81bool DelphesHepMCReader::EventReady()
82{
83  return (fVertexCounter == 0) && (fInCounter == 0) && (fOutCounter == 0);
84}
85
86//---------------------------------------------------------------------------
87
88bool DelphesHepMCReader::ReadBlock(DelphesFactory *factory,
89  TObjArray *allParticleOutputArray,
90  TObjArray *stableParticleOutputArray,
91  TObjArray *partonOutputArray)
92{
93  map< int, pair< int, int > >::iterator itMotherMap;
94  map< int, pair< int, int > >::iterator itDaughterMap;
95  char key;
96  int rc;
97
98  if(!fgets(fBuffer, kBufferSize, fInputFile)) return kFALSE;
99
100  DelphesStream bufferStream(fBuffer + 1);
101
102  key = fBuffer[0];
103
104  if(key == 'E')
105  {
106    Clear();
107
108    rc = bufferStream.ReadInt(fEventNumber)
109      && bufferStream.ReadInt(fMPI)
110      && bufferStream.ReadDbl(fScale)
111      && bufferStream.ReadDbl(fAlphaQCD)
112      && bufferStream.ReadDbl(fAlphaQED)
113      && bufferStream.ReadInt(fProcessID)
114      && bufferStream.ReadInt(fSignalCode)
115      && bufferStream.ReadInt(fVertexCounter);
116
117    if(!rc)
118    {
119      cerr << "** ERROR: " << "invalid event format" << endl;
120      return kFALSE;
121    }
122  }
123  else if(key == 'F')
124  {
125    rc = bufferStream.ReadInt(fID1)
126      && bufferStream.ReadInt(fID2)
127      && bufferStream.ReadDbl(fX1)
128      && bufferStream.ReadDbl(fX2)
129      && bufferStream.ReadDbl(fScalePDF)
130      && bufferStream.ReadDbl(fPDF1)
131      && bufferStream.ReadDbl(fPDF2);
132
133    if(!rc)
134    {
135      cerr << "** ERROR: " << "invalid PDF format" << endl;
136      return kFALSE;
137    }
138  }
139  else if(key == 'V' && fVertexCounter > 0)
140  {
141    rc = bufferStream.ReadInt(fOutVertexCode)
142      && bufferStream.ReadInt(fVertexID)
143      && bufferStream.ReadDbl(fX)
144      && bufferStream.ReadDbl(fY)
145      && bufferStream.ReadDbl(fZ)
146      && bufferStream.ReadDbl(fT)
147      && bufferStream.ReadInt(fInCounter)
148      && bufferStream.ReadInt(fOutCounter);
149
150    if(!rc)
151    {
152      cerr << "** ERROR: " << "invalid vertex format" << endl;
153      return kFALSE;
154    }
155    --fVertexCounter;
156  }
157  else if(key == 'P' && fOutCounter > 0)
158  {
159    rc = bufferStream.ReadInt(fParticleCode)
160      && bufferStream.ReadInt(fPID)
161      && bufferStream.ReadDbl(fPx)
162      && bufferStream.ReadDbl(fPy)
163      && bufferStream.ReadDbl(fPz)
164      && bufferStream.ReadDbl(fE)
165      && bufferStream.ReadDbl(fMass)
166      && bufferStream.ReadInt(fStatus)
167      && bufferStream.ReadDbl(fTheta)
168      && bufferStream.ReadDbl(fPhi)
169      && bufferStream.ReadInt(fInVertexCode);
170
171    if(!rc)
172    {
173      cerr << "** ERROR: " << "invalid particle format" << endl;
174      return kFALSE;
175    }
176
177    if(fInVertexCode < 0)
178    {
179      itMotherMap = fMotherMap.find(fInVertexCode);
180      if(itMotherMap == fMotherMap.end())
181      {
182        fMotherMap[fInVertexCode] = make_pair(fParticleCounter, -1);
183      }
184      else
185      {
186        itMotherMap->second.second = fParticleCounter;
187      }
188    }
189
190    if(fInCounter <= 0)
191    {
192      itDaughterMap = fDaughterMap.find(fOutVertexCode);
193      if(itDaughterMap == fDaughterMap.end())
194      {
195        fDaughterMap[fOutVertexCode] = make_pair(fParticleCounter, fParticleCounter);
196      }
197      else
198      {
199        itDaughterMap->second.second = fParticleCounter;
200      }
201    }
202
203    AnalyzeParticle(factory, allParticleOutputArray,
204      stableParticleOutputArray, partonOutputArray);
205
206    if(fInCounter > 0)
207    {
208      --fInCounter;
209    }
210    else
211    {
212      --fOutCounter;
213    }
214
215    ++fParticleCounter;
216  }
217
218  if(EventReady())
219  {
220    FinalizeParticles(allParticleOutputArray);
221  }
222
223  return kTRUE;
224}
225
226//---------------------------------------------------------------------------
227
228void DelphesHepMCReader::AnalyzeEvent(ExRootTreeBranch *branch, long long eventNumber,
229  TStopwatch *readStopWatch, TStopwatch *procStopWatch)
230{
231  HepMCEvent *element;
232
233  element = static_cast<HepMCEvent *>(branch->NewEntry());
234  element->Number = fEventNumber;
235
236  element->ProcessID = fProcessID;
237  element->MPI = fMPI;
238  element->Scale = fScale;
239  element->AlphaQED = fAlphaQED;
240  element->AlphaQCD = fAlphaQCD;
241
242  element->ID1 = fID1;
243  element->ID2 = fID2;
244  element->X1 = fX1;
245  element->X2 = fX2;
246  element->ScalePDF = fScalePDF;
247  element->PDF1 = fPDF1;
248  element->PDF2 = fPDF2;
249
250  element->ReadTime = readStopWatch->RealTime();
251  element->ProcTime = procStopWatch->RealTime();
252}
253
254//---------------------------------------------------------------------------
255
256void DelphesHepMCReader::AnalyzeParticle(DelphesFactory *factory,
257  TObjArray *allParticleOutputArray,
258  TObjArray *stableParticleOutputArray,
259  TObjArray *partonOutputArray)
260{
261  Candidate *candidate;
262  TParticlePDG *pdgParticle;
263  int pdgCode;
264
265  candidate = factory->NewCandidate();
266
267  candidate->PID = fPID;
268  pdgCode = TMath::Abs(candidate->PID);
269
270  candidate->Status = fStatus;
271
272  pdgParticle = fPDG->GetParticle(fPID);
273  candidate->Charge = pdgParticle ? int(pdgParticle->Charge()/3.0) : -999;
274  candidate->Mass = pdgParticle ? pdgParticle->Mass() : -999.9;
275
276  candidate->Momentum.SetPxPyPzE(fPx, fPy, fPz, fE);
277
278  candidate->M2 = 1;
279  candidate->D2 = 1;
280  if(fInCounter > 0)
281  {
282    candidate->M1 = 1;
283    candidate->Position.SetXYZT(0.0, 0.0, 0.0, 0.0);
284  }
285  else
286  {
287    candidate->M1 = fOutVertexCode;
288    candidate->Position.SetXYZT(fX, fY, fZ, fT);
289  }
290  if(fInVertexCode < 0)
291  {
292    candidate->D1 = fInVertexCode;
293  }
294  else
295  {
296    candidate->D1 = 1;
297  }
298
299  allParticleOutputArray->Add(candidate);
300
301  if(!pdgParticle) return;
302
303  if(fStatus == 1 && pdgParticle->Stable())
304  {
305    stableParticleOutputArray->Add(candidate);
306  }
307  else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
308  {
309    partonOutputArray->Add(candidate);
310  }
311}
312
313//---------------------------------------------------------------------------
314
315void DelphesHepMCReader::FinalizeParticles(TObjArray *allParticleOutputArray)
316{
317  Candidate *candidate;
318  map< int, pair< int, int > >::iterator itMotherMap;
319  map< int, pair< int, int > >::iterator itDaughterMap;
320  int i;
321
322  for(i = 0; i < allParticleOutputArray->GetEntriesFast(); ++i)
323  {
324    candidate = static_cast<Candidate *>(allParticleOutputArray->At(i));
325
326    if(candidate->M1 > 0)
327    {
328      candidate->M1 = -1;
329      candidate->M2 = -1;
330    }
331    else
332    {
333      itMotherMap = fMotherMap.find(candidate->M1);
334      if(itMotherMap == fMotherMap.end())
335      {
336        candidate->M1 = -1;
337        candidate->M2 = -1;
338      }
339      else
340      {
341        candidate->M1 = itMotherMap->second.first;
342        candidate->M2 = itMotherMap->second.second;
343      }
344    }
345    if(candidate->D1 > 0)
346    {
347      candidate->D1 = -1;
348      candidate->D2 = -1;
349    }
350    else
351    {
352      itDaughterMap = fDaughterMap.find(candidate->D1);
353      if(itDaughterMap == fDaughterMap.end())
354      {
355        candidate->D1 = -1;
356        candidate->D2 = -1;
357      }
358      else
359      {
360        candidate->D1 = itDaughterMap->second.first;
361        candidate->D2 = itDaughterMap->second.second;
362      }
363    }
364  }
365}
366
367//---------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.