Fork me on GitHub

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

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

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