Fork me on GitHub

source: git/classes/DelphesSTDHEPReader.cc@ f59b6d75

ImprovedOutputFile Timing dual_readout llp
Last change on this file since f59b6d75 was f59b6d75, checked in by pavel <pavel@…>, 11 years ago

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