Fork me on GitHub

source: git/classes/DelphesHepMCReader.cc @ 1fa50c2

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

fix GPLv3 header

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