Fork me on GitHub

source: git/classes/DelphesHepMCReader.cc @ cab38f6

ImprovedOutputFileTimingdual_readoutllp
Last change on this file since cab38f6 was cab38f6, checked in by Pavel Demin <pavel.demin@…>, 6 years ago

remove svn tags and fix formatting

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