Fork me on GitHub

source: git/classes/DelphesSTDHEPReader.cc

Last change on this file was 341014c, checked in by Pavel Demin <pavel-demin@…>, 18 months ago

apply .clang-format to all .h, .cc and .cpp files

  • Property mode set to 100644
File size: 11.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/** \class DelphesSTDHEPReader
20 *
21 *  Reads STDHEP file
22 *
23 *  \author P. Demin - UCL, Louvain-la-Neuve
24 *
25 */
26
27#include "classes/DelphesSTDHEPReader.h"
28
29#include <iostream>
30#include <sstream>
31#include <stdexcept>
32
33#include <errno.h>
34#include <stdint.h>
35#include <stdio.h>
36#include <string.h>
37
38#include "TDatabasePDG.h"
39#include "TLorentzVector.h"
40#include "TObjArray.h"
41#include "TParticlePDG.h"
42#include "TStopwatch.h"
43
44#include "classes/DelphesClasses.h"
45#include "classes/DelphesFactory.h"
46#include "classes/DelphesXDRReader.h"
47
48#include "ExRootAnalysis/ExRootTreeBranch.h"
49
50using namespace std;
51
52static const int kBufferSize = 1000000;
53
54//---------------------------------------------------------------------------
55
56DelphesSTDHEPReader::DelphesSTDHEPReader() :
57  fInputFile(0), fBuffer(0), fPDG(0), fBlockType(-1)
58{
59  fBuffer = new uint8_t[kBufferSize * 96 + 24];
60
61  fPDG = TDatabasePDG::Instance();
62}
63
64//---------------------------------------------------------------------------
65
66DelphesSTDHEPReader::~DelphesSTDHEPReader()
67{
68  if(fBuffer) delete fBuffer;
69}
70
71//---------------------------------------------------------------------------
72
73void DelphesSTDHEPReader::SetInputFile(FILE *inputFile)
74{
75  fInputFile = inputFile;
76  fReader[0].SetFile(inputFile);
77}
78
79//---------------------------------------------------------------------------
80
81void DelphesSTDHEPReader::Clear()
82{
83  fBlockType = -1;
84}
85
86//---------------------------------------------------------------------------
87
88bool DelphesSTDHEPReader::EventReady()
89{
90  return (fBlockType == MCFIO_STDHEP) || (fBlockType == MCFIO_STDHEP4);
91}
92
93//---------------------------------------------------------------------------
94
95bool DelphesSTDHEPReader::ReadBlock(DelphesFactory *factory,
96  TObjArray *allParticleOutputArray,
97  TObjArray *stableParticleOutputArray,
98  TObjArray *partonOutputArray)
99{
100  fReader[0].ReadValue(&fBlockType, 4);
101
102  if(feof(fInputFile)) return kFALSE;
103
104  SkipBytes(4);
105
106  if(fBlockType == FILEHEADER)
107  {
108    ReadFileHeader();
109  }
110  else if(fBlockType == EVENTTABLE)
111  {
112    ReadEventTable();
113  }
114  else if(fBlockType == EVENTHEADER)
115  {
116    ReadEventHeader();
117  }
118  else if(fBlockType == MCFIO_STDHEPBEG || fBlockType == MCFIO_STDHEPEND)
119  {
120    ReadSTDCM1();
121  }
122  else if(fBlockType == MCFIO_STDHEP)
123  {
124    ReadSTDHEP();
125    AnalyzeParticles(factory, allParticleOutputArray,
126      stableParticleOutputArray, partonOutputArray);
127  }
128  else if(fBlockType == MCFIO_STDHEP4)
129  {
130    ReadSTDHEP();
131    AnalyzeParticles(factory, allParticleOutputArray,
132      stableParticleOutputArray, partonOutputArray);
133    ReadSTDHEP4();
134  }
135  else
136  {
137    throw runtime_error("Unsupported block type.");
138  }
139
140  return kTRUE;
141}
142
143//---------------------------------------------------------------------------
144
145void DelphesSTDHEPReader::SkipBytes(int size)
146{
147  int rc;
148  int rndup;
149
150  rndup = size % 4;
151  if(rndup > 0)
152  {
153    rndup = 4 - rndup;
154  }
155
156  rc = fseek(fInputFile, size + rndup, SEEK_CUR);
157
158  if(rc != 0 && errno == ESPIPE)
159  {
160    fReader[0].ReadRaw(fBuffer, size);
161  }
162}
163
164//---------------------------------------------------------------------------
165
166void DelphesSTDHEPReader::SkipArray(int elsize)
167{
168  uint32_t size;
169  fReader[0].ReadValue(&size, 4);
170  SkipBytes(size * elsize);
171}
172
173//---------------------------------------------------------------------------
174
175void DelphesSTDHEPReader::ReadFileHeader()
176{
177  uint32_t i;
178  enum STDHEPVersion
179  {
180    UNKNOWN,
181    V1,
182    V2,
183    V21
184  } version;
185
186  // version
187  fReader[0].ReadString(fBuffer, 100);
188  if(fBuffer[0] == '\0' || fBuffer[1] == '\0')
189    version = UNKNOWN;
190  else if(fBuffer[0] == '1')
191    version = V1;
192  else if(strncmp((char *)fBuffer, "2.01", 4) == 0)
193    version = V21;
194  else if(fBuffer[0] == '2')
195    version = V2;
196  else
197    version = UNKNOWN;
198
199  if(version == UNKNOWN)
200  {
201    throw runtime_error("Unknown file format version.");
202  }
203
204  SkipArray(1);
205  SkipArray(1);
206  SkipArray(1);
207
208  if(version == V21)
209  {
210    SkipArray(1);
211  }
212
213  // Expected number of events
214  SkipBytes(4);
215
216  // Number of events
217  fReader[0].ReadValue(&fEntries, 4);
218
219  SkipBytes(8);
220
221  // Number of blocks
222  uint32_t nBlocks = 0;
223  fReader[0].ReadValue(&nBlocks, 4);
224
225  // Number of NTuples
226  uint32_t nNTuples = 0;
227  if(version != V1)
228  {
229    fReader[0].ReadValue(&nNTuples, 4);
230  }
231
232  if(nNTuples != 0)
233  {
234    throw runtime_error("Files containing n-tuples are not supported.");
235  }
236
237  // Processing blocks extraction
238  if(nBlocks != 0)
239  {
240    SkipArray(4);
241
242    for(i = 0; i < nBlocks; i++)
243    {
244      SkipArray(1);
245    }
246  }
247}
248
249//---------------------------------------------------------------------------
250
251void DelphesSTDHEPReader::ReadEventTable()
252{
253  // version
254  fReader[0].ReadString(fBuffer, 100);
255  if(strncmp((char *)fBuffer, "1.00", 4) == 0)
256  {
257    SkipBytes(8);
258
259    SkipArray(4);
260    SkipArray(4);
261    SkipArray(4);
262    SkipArray(4);
263    SkipArray(4);
264  }
265  else if(strncmp((char *)fBuffer, "2.00", 4) == 0)
266  {
267    SkipBytes(12);
268
269    SkipArray(4);
270    SkipArray(4);
271    SkipArray(4);
272    SkipArray(4);
273    SkipArray(8);
274  }
275}
276
277//---------------------------------------------------------------------------
278
279void DelphesSTDHEPReader::ReadEventHeader()
280{
281  bool skipNTuples = false;
282  int skipSize = 4;
283
284  // version
285  fReader[0].ReadString(fBuffer, 100);
286  if(strncmp((char *)fBuffer, "2.00", 4) == 0)
287  {
288    skipNTuples = true;
289  }
290  else if(strncmp((char *)fBuffer, "3.00", 4) == 0)
291  {
292    skipNTuples = true;
293    skipSize = 8;
294  }
295
296  SkipBytes(20);
297
298  uint32_t dimBlocks = 0;
299  fReader[0].ReadValue(&dimBlocks, 4);
300
301  uint32_t dimNTuples = 0;
302  if(skipNTuples)
303  {
304    SkipBytes(4);
305    fReader[0].ReadValue(&dimNTuples, 4);
306  }
307
308  // Processing blocks extraction
309  if(dimBlocks > 0)
310  {
311    SkipArray(4);
312    SkipArray(skipSize);
313  }
314
315  // Processing blocks extraction
316  if(skipNTuples && dimNTuples > 0)
317  {
318    SkipArray(4);
319    SkipArray(skipSize);
320  }
321}
322
323//---------------------------------------------------------------------------
324
325void DelphesSTDHEPReader::ReadSTDCM1()
326{
327  // version
328  fReader[0].ReadString(fBuffer, 100);
329
330  // skip 5*4 + 2*8 = 36 bytes
331  SkipBytes(36);
332
333  if((strncmp((char *)fBuffer, "1.", 2) == 0) || (strncmp((char *)fBuffer, "2.", 2) == 0) ||
334     (strncmp((char *)fBuffer, "3.", 2) == 0) || (strncmp((char *)fBuffer, "4.", 2) == 0) ||
335     (strncmp((char *)fBuffer, "5.00", 4) == 0))
336  {
337    return;
338  }
339
340  SkipArray(1);
341  SkipArray(1);
342
343  if(strncmp((char *)fBuffer, "5.01", 4) == 0)
344  {
345    return;
346  }
347
348  SkipBytes(4);
349}
350
351//---------------------------------------------------------------------------
352
353void DelphesSTDHEPReader::ReadSTDHEP()
354{
355  uint32_t idhepSize, isthepSize, jmohepSize, jdahepSize, phepSize, vhepSize;
356
357  // version
358  fReader[0].ReadString(fBuffer, 100);
359
360  // Extracting the event number
361  fReader[0].ReadValue(&fEventNumber, 4);
362
363  // Extracting the number of particles
364  fReader[0].ReadValue(&fEventSize, 4);
365
366  if(fEventSize >= kBufferSize)
367  {
368    throw runtime_error("too many particles in event");
369  }
370
371  // 4*n + 4*n + 8*n + 8*n + 40*n + 32*n +
372  // 4 + 4 + 4 + 4 + 4 + 4 = 96*n + 24
373
374  fReader[0].ReadRaw(fBuffer, 96 * fEventSize + 24);
375
376  fReader[1].SetBuffer(fBuffer);
377  fReader[2].SetBuffer(fBuffer + 4 * 1 + 4 * 1 * fEventSize);
378  fReader[3].SetBuffer(fBuffer + 4 * 2 + 4 * 2 * fEventSize);
379  fReader[4].SetBuffer(fBuffer + 4 * 3 + 4 * 4 * fEventSize);
380  fReader[5].SetBuffer(fBuffer + 4 * 4 + 4 * 6 * fEventSize);
381  fReader[6].SetBuffer(fBuffer + 4 * 5 + 4 * 16 * fEventSize);
382
383  fReader[1].ReadValue(&idhepSize, 4);
384  fReader[2].ReadValue(&isthepSize, 4);
385  fReader[3].ReadValue(&jmohepSize, 4);
386  fReader[4].ReadValue(&jdahepSize, 4);
387  fReader[5].ReadValue(&phepSize, 4);
388  fReader[6].ReadValue(&vhepSize, 4);
389
390  if(fEventSize < 0 ||
391     fEventSize != (int)idhepSize      || fEventSize != (int)isthepSize     ||
392     (2*fEventSize) != (int)jmohepSize || (2*fEventSize) != (int)jdahepSize ||
393     (5*fEventSize) != (int)phepSize   || (4*fEventSize) != (int)vhepSize)
394  {
395    throw runtime_error("Inconsistent size of arrays. File is probably corrupted.");
396  }
397
398  fWeight = 1.0;
399  fAlphaQED = 0.0;
400  fAlphaQCD = 0.0;
401  fScaleSize = 0;
402  memset(fScale, 0, 10 * sizeof(double));
403}
404
405//---------------------------------------------------------------------------
406
407void DelphesSTDHEPReader::ReadSTDHEP4()
408{
409  uint32_t number;
410
411  // Extracting the event weight
412  fReader[0].ReadValue(&fWeight, 8);
413
414  // Extracting alpha QED
415  fReader[0].ReadValue(&fAlphaQED, 8);
416
417  // Extracting alpha QCD
418  fReader[0].ReadValue(&fAlphaQCD, 8);
419
420  // Extracting the event scale
421  fReader[0].ReadValue(&fScaleSize, 4);
422  for(number = 0; number < fScaleSize; ++number)
423  {
424    fReader[0].ReadValue(&fScale[number], 8);
425  }
426
427  SkipArray(8);
428  SkipArray(4);
429
430  SkipBytes(4);
431}
432
433//---------------------------------------------------------------------------
434
435void DelphesSTDHEPReader::AnalyzeEvent(ExRootTreeBranch *branch, long long eventNumber,
436  TStopwatch *readStopWatch, TStopwatch *procStopWatch)
437{
438  LHEFEvent *element;
439
440  element = static_cast<LHEFEvent *>(branch->NewEntry());
441
442  element->Number = fEventNumber;
443
444  element->ProcessID = 0;
445
446  element->Weight = fWeight;
447  element->ScalePDF = fScale[0];
448  element->AlphaQED = fAlphaQED;
449  element->AlphaQCD = fAlphaQCD;
450
451  element->ReadTime = readStopWatch->RealTime();
452  element->ProcTime = procStopWatch->RealTime();
453}
454
455//---------------------------------------------------------------------------
456
457void DelphesSTDHEPReader::AnalyzeParticles(DelphesFactory *factory,
458  TObjArray *allParticleOutputArray,
459  TObjArray *stableParticleOutputArray,
460  TObjArray *partonOutputArray)
461{
462  Candidate *candidate;
463  TParticlePDG *pdgParticle;
464  int pdgCode;
465
466  int number;
467  int32_t pid, status, m1, m2, d1, d2;
468  double px, py, pz, e, mass;
469  double x, y, z, t;
470
471  for(number = 0; number < fEventSize; ++number)
472  {
473    fReader[1].ReadValue(&status, 4);
474    fReader[2].ReadValue(&pid, 4);
475    fReader[3].ReadValue(&m1, 4);
476    fReader[3].ReadValue(&m2, 4);
477    fReader[4].ReadValue(&d1, 4);
478    fReader[4].ReadValue(&d2, 4);
479
480    fReader[5].ReadValue(&px, 8);
481    fReader[5].ReadValue(&py, 8);
482    fReader[5].ReadValue(&pz, 8);
483    fReader[5].ReadValue(&e, 8);
484    fReader[5].ReadValue(&mass, 8);
485
486    fReader[6].ReadValue(&x, 8);
487    fReader[6].ReadValue(&y, 8);
488    fReader[6].ReadValue(&z, 8);
489    fReader[6].ReadValue(&t, 8);
490
491    candidate = factory->NewCandidate();
492
493    candidate->PID = pid;
494    pdgCode = TMath::Abs(candidate->PID);
495
496    candidate->Status = status;
497
498    candidate->M1 = m1 - 1;
499    candidate->M2 = m2 - 1;
500
501    candidate->D1 = d1 - 1;
502    candidate->D2 = d2 - 1;
503
504    pdgParticle = fPDG->GetParticle(pid);
505    candidate->Charge = pdgParticle ? int(pdgParticle->Charge() / 3.0) : -999;
506    candidate->Mass = mass;
507
508    candidate->Momentum.SetPxPyPzE(px, py, pz, e);
509
510    candidate->Position.SetXYZT(x, y, z, t);
511
512    allParticleOutputArray->Add(candidate);
513
514    if(!pdgParticle) continue;
515
516    if(status == 1)
517    {
518      stableParticleOutputArray->Add(candidate);
519    }
520    else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
521    {
522      partonOutputArray->Add(candidate);
523    }
524  }
525}
526
527//---------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.