Fork me on GitHub

source: git/classes/DelphesSTDHEPReader.cc

Last change on this file was 341014c, checked in by Pavel Demin <pavel-demin@…>, 6 years ago

apply .clang-format to all .h, .cc and .cpp files

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