Fork me on GitHub

source: svn/trunk/classes/DelphesSTDHEPReader.cc@ 1125

Last change on this file since 1125 was 1120, checked in by Pavel Demin, 11 years ago

fix STDHEP reader
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id Revision Date
File size: 10.7 KB
Line 
1
2/** \class DelphesSTDHEPReader
3 *
4 * Reads STDHEP file
5 *
6 *
7 * $Date: 2013-05-21 20:03:37 +0000 (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.