Fork me on GitHub

source: git/classes/DelphesSTDHEPReader.cc@ c04eb8a

ImprovedOutputFile Timing dual_readout llp
Last change on this file since c04eb8a was 1fa50c2, checked in by Pavel Demin <pavel.demin@…>, 10 years ago

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