source: HiSusy/trunk/Delphes/Delphes-3.0.9/classes/DelphesSTDHEPReader.cc @ 5

Last change on this file since 5 was 5, checked in by zerwas, 11 years ago

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