Fork me on GitHub

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

Last change on this file since 1096 was 1053, checked in by Pavel Demin, 12 years ago

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