GEMC  2.2
Geant4 Monte-Carlo Framework
lStdHep.cc
Go to the documentation of this file.
1 //
3 // This class is based on the light-weight XDR class lXDR,
4 // and parses/writes StdHep files. It was mainly written,
5 // to provide a faster alternative to the more cumbersome
6 // methods using mcfio in StdHep.
7 //
8 // W.G.J. Langeveld, 24 May 2002
9 //
10 // Copied from the LCIO library.
11 //
13 #include "lStdHep.hh"
14 #include "string.h"
15 #include "stdlib.h"
16 
17 namespace UTIL{
18 
20 //
21 // The main lStdHep class.
22 //
24 
25 //
26 // Constructors, destructor
27 // ------------------------
28 // Constructor opens file, destructor closes file. Once opened for
29 // reading, the file cannot be written to, and v.v.
30 //
31 lStdHep::lStdHep(const char *filename, bool open_for_write) :
32  lXDR(filename, open_for_write),
33  ntot(0),version(0),title(0),comment(0),date(0),closingDate(0),numevts_expect(0),numevts(0),
34  firstTable(0),dimTable(0),nNTuples(0),nBlocks(0),blockIds(0),blockNames(0)
35 {
36  if (open_for_write) {
37  setError(LSH_NOTSUPPORTED);
38  }
39  else {
40  readFileHeader();
41  }
42  return;
43 }
44 
45 lStdHep::~lStdHep()
46 {
47  delete [] version;
48  delete [] date;
49  delete [] closingDate;
50  delete [] comment;
51  delete [] title;
52  delete [] blockIds;
53  if (blockNames) {
54  for (int i = 0; i < nBlocks; i++) {
55  delete [] blockNames[i];
56  }
57  delete [] blockNames;
58  }
59  return;
60 }
61 
62 void lStdHep::printFileHeader(FILE *fp)
63 {
64  if (fp == 0) fp = stdout;
65 
66  fprintf(fp, "====== File Header ===========\n");
67  fprintf(fp, " total blocks: %ld\n", ntot);
68  fprintf(fp, " version: %s\n", version);
69  fprintf(fp, " title: %s\n", title);
70  fprintf(fp, " comment: %s\n", comment);
71  fprintf(fp, " date: %s", date);
72  fprintf(fp, " closing date: %s", closingDate);
73  fprintf(fp, " expected events: %ld\n", numevts_expect);
74  fprintf(fp, " events: %ld\n", numevts);
75  fprintf(fp, " firstTable: %ld\n", firstTable);
76  fprintf(fp, " dimTable: %ld\n", dimTable);
77  fprintf(fp, " nNTuples: %ld\n", nNTuples);
78  fprintf(fp, " nBlocks: %ld\n", nBlocks);
79  if (nBlocks) fprintf(fp, " block names:\n");
80  for (int i = 0; i < nBlocks; i++) {
81  fprintf(fp, " : %s\n", blockNames[i]);
82  }
83  fprintf(fp, "=============================\n");
84  return;
85 }
86 
87 void lStdHep::printEventTable(FILE *fp)
88 {
89  if (fp == 0) fp = stdout;
90  eventTable.print(fp);
91  return;
92 }
93 
94 void lStdHep::printEventHeader(FILE *fp)
95 {
96  if (fp == 0) fp = stdout;
97  event.printHeader(fp);
98  return;
99 }
100 
101 void lStdHep::printEvent(FILE *fp)
102 {
103  if (fp == 0) fp = stdout;
104  event.print(fp);
105  return;
106 }
107 
108 void lStdHep::printTrack(int i, FILE *fp)
109 {
110  if (fp == 0) fp = stdout;
111  if (i < event.nhep) {
112  fprintf(fp, " Track: id: %ld, vtx: (%g, %g, %g, %g), mom: (%g, %g, %g, %g, %g)\n",
113  pid(i), X(i), Y(i), Z(i), T(i), Px(i), Py(i), Pz(i), E(i), M(i));
114  fprintf(fp, " Track: wgt: %g, alpha QED: %g, alpha QCD: %g, idrup: %ld\n",
115  eventweight(), alphaQED(), alphaQCD(), idrup());
116  }
117  return;
118 }
119 
120 void lStdHep::printBeginRunRecord(FILE *fp)
121 {
122  if (fp == 0) fp = stdout;
123  return;
124 }
125 
126 void lStdHep::printEndRunRecord(FILE *fp)
127 {
128  if (fp == 0) fp = stdout;
129  return;
130 }
131 
132 bool lStdHep::more(void)
133 {
134  return(getError() == LSH_SUCCESS);
135 }
136 
137 long lStdHep::readEvent(void)
138 {
139 //
140 // Look for an event or an event table
141 //
142  event.isEmpty = 1;
143  while (1) {
144  if (eventTable.ievt < eventTable.numEvts) {
145  if (filePosition(eventTable.ptrEvents[eventTable.ievt]) !=
146  eventTable.ptrEvents[eventTable.ievt]) return(getError());
147 
148  if (event.read(*this) != LSH_SUCCESS) return(getError());
149  eventTable.ievt++;
150 
151  if (event.isEmpty) continue;
152  return(getError());
153  }
154 
155  eventTable.isEmpty = 1;
156  while (eventTable.isEmpty) {
157  if (eventTable.nextlocator == -2) {
158 //
159 // This was the last event table, signal quitting. Not an error.
160 //
161  setError(LSH_ENDOFFILE);
162  return(getError());
163  }
164  else if (eventTable.nextlocator == -1) {
165  setError(LSH_EVTABLECORRUPT);
166  return(getError());
167  }
168 //
169 // Go to the next event table
170 //
171  if (filePosition(eventTable.nextlocator) != eventTable.nextlocator) return(getError());
172  if (eventTable.read(*this) != LSH_SUCCESS) return(getError());
173  }
174  }
175  return(getError());
176 }
177 
178 long lStdHep::getEvent(lStdEvent &lse) const
179 {
180  if (long status = getError() != LSH_SUCCESS) return(status);
181 
182  lse.evtNum = event.nevhep;
183 
184  lse.clear();
185  for (int i = 0; i < event.nhep; i++) {
186  lStdTrack lst;
187  lst.X = X(i);
188  lst.Y = Y(i);
189  lst.Z = Z(i);
190  lst.T = T(i);
191  lst.Px = Px(i);
192  lst.Py = Py(i);
193  lst.Pz = Pz(i);
194  lst.E = E(i);
195  lst.M = M(i);
196  lst.pid = pid(i);
197  lst.status = status(i);
198  lst.mother1 = mother1(i);
199  lst.mother2 = mother2(i);
200  lst.daughter1 = daughter1(i);
201  lst.daughter2 = daughter2(i);
202  lse.push_back(lst);
203  }
204  return(LSH_SUCCESS);
205 }
206 
207 long lStdHep::readEvent(lStdEvent &lse)
208 {
209  long status = readEvent();
210  if (status != LSH_SUCCESS) return(status);
211  return(getEvent(lse));
212 }
213 
214 long lStdHep::readFileHeader(void)
215 {
216  long len, blockid;
217 
218  blockid = readLong();
219  if (blockid != LSH_FILEHEADER) {
220  setError(LSH_BLOCKERROR);
221  return(getError());
222  }
223  ntot = readLong();
224  version = readString(len);
225 
226  title = readString(len);
227  comment = readString(len);
228  date = readString(len);
229  if ((strcmp(version, "2.00") == 0) || (strcmp(version, "1.00") == 0)) {
230  closingDate = new char[len + 1];
231  strcpy((char *) closingDate, date);
232  }
233  else {
234  closingDate = readString(len);
235  }
236 
237  numevts_expect = readLong();
238  numevts = readLong();
239  firstTable = readLong();
240  dimTable = readLong();
241  nBlocks = readLong();
242  if (*version != '2') {
243  nNTuples = 0;
244  }
245  else {
246  nNTuples = readLong();
247  }
248 
249  blockIds = readLongArray(nBlocks);
250  blockNames = new const char *[nBlocks];
251  for (int i = 0; i < nBlocks; i++) blockNames[i] = readString(len);
252  if (nNTuples > 0) setError(LSH_NOTSUPPORTED);
253 //
254 // Read the first event table
255 //
256  eventTable.read(*this);
257  return(getError());
258 }
259 
260 long lStdHep::writeEvent(void)
261 {
262  return(LSH_NOTSUPPORTED);
263 }
264 
265 long lStdHep::setEvent(const lStdEvent &lse)
266 {
267 // ***************set up event buffer!
268  setNTracks(lse.size());
269  setEvtNum(lse.evtNum);
270 
271  for (int i = 0; i < event.nhep; i++) {
272  setX (i, lse[i].X);
273  setY (i, lse[i].Y);
274  setZ (i, lse[i].Z);
275  setT (i, lse[i].T);
276  setPx (i, lse[i].Px);
277  setPy (i, lse[i].Py);
278  setPz (i, lse[i].Pz);
279  setE (i, lse[i].E);
280  setM (i, lse[i].M);
281  setPid (i, lse[i].pid);
282  setStatus (i, lse[i].status);
283  setMother1 (i, lse[i].mother1);
284  setMother2 (i, lse[i].mother2);
285  setDaughter1(i, lse[i].daughter1);
286  setDaughter2(i, lse[i].daughter2);
287  }
288  return(LSH_SUCCESS);
289 }
290 
291 long lStdHep::writeEvent(lStdEvent &lse)
292 {
293  long status = writeEvent();
294  if (status != LSH_SUCCESS) return(status);
295  return(setEvent(lse));
296 }
297 
298 lStdHep::EventTable::EventTable() :
299  isEmpty(1), ievt(0), blockid(0), ntot(0), version(0),
300  nextlocator(-3), numEvts(0), evtnums(0),
301  storenums(0), runnums(0), trigMasks(0), ptrEvents(0)
302 {
303  return;
304 }
305 
306 lStdHep::EventTable::~EventTable()
307 {
308  cleanup();
309  return;
310 }
311 
312 void lStdHep::EventTable::cleanup(void)
313 {
314  delete [] version; version = 0;
315  delete [] evtnums; evtnums = 0;
316  delete [] storenums; storenums = 0;
317  delete [] runnums; runnums = 0;
318  delete [] trigMasks; trigMasks = 0;
319  delete [] ptrEvents; ptrEvents = 0;
320  isEmpty = 1;
321  ievt = ntot = blockid = numEvts = 0; // leave nextlocator alone!
322  return;
323 }
324 
325 long lStdHep::EventTable::read(lStdHep &ls)
326 {
327  long len;
328 
329  cleanup();
330 
331  blockid = ls.readLong();
332  ntot = ls.readLong();
333  version = ls.readString(len);
334 
335  if (blockid != LSH_EVENTTABLE) {
336  ls.setError(LSH_NOEVENTTABLE);
337  return(ls.getError());
338  }
339  nextlocator = ls.readLong();
340  numEvts = ls.readLong();
341  evtnums = ls.readLongArray(len);
342  storenums = ls.readLongArray(len);
343  runnums = ls.readLongArray(len);
344  trigMasks = ls.readLongArray(len);
345  ptrEvents = ls.readLongArray(len);
346  if (numEvts > 0) isEmpty = 0;
347  return(ls.getError());
348 }
349 
350 long lStdHep::EventTable::print(FILE *fp)
351 {
352  fprintf(fp, " EventTable: blockid: %ld, ntot: %ld, version: %s\n", blockid, ntot, version);
353  fprintf(fp, " EventTable: nextlocator: %ld, numEvts: %ld\n", nextlocator, numEvts);
354  for (int i = 0; i < numEvts; i++) {
355  fprintf(fp, " EventTable: %d: evtnums %ld storenums %ld runnums %ld trigMasks %ld ptrEvents %ld\n",
356  i, evtnums[i], storenums[i], runnums[i], trigMasks[i], ptrEvents[i]);
357  if (i == 10) {
358  fprintf(fp, " EventTable: etc.\n");
359  break;
360  }
361  }
362  return(0);
363 }
364 
365 lStdHep::Event::Event() :
366  isEmpty(0), blockid(0),ntot(0),version(0),evtnum(0),storenum(0),runnum(0),
367  trigMask(0),nBlocks(0),dimBlocks(0),nNTuples(0),dimNTuples(0),blockIds(0),
368  ptrBlocks(0),nevhep(0),nhep(0),isthep(0),idhep(0),jmohep(0),jdahep(0),phep(0),
369  vhep(0),eventweight(0),alphaqed(0),alphaqcd(0),scale(0),spin(0),colorflow(0),idrup(0),
370  bnevtreq(0),bnevtgen(0),bnevtwrt(0),bstdecom(0),bstdxsec(0),bstdseed1(0),bstdseed2(0),
371  enevtreq(0),enevtgen(0),enevtwrt(0),estdecom(0),estdxsec(0),estdseed1(0),estdseed2(0)
372 
373 {
374  return;
375 }
376 
377 lStdHep::Event::~Event()
378 {
379  cleanup();
380 }
381 
382 void lStdHep::Event::cleanup(void)
383 {
384  delete [] version; version = 0;
385  delete [] ptrBlocks; ptrBlocks = 0;
386  delete [] blockIds; blockIds = 0;
387  delete [] isthep; isthep = 0;
388  delete [] idhep; idhep = 0;
389  delete [] jmohep; jmohep = 0;
390  delete [] jdahep; jdahep = 0;
391  delete [] phep; phep = 0;
392  delete [] vhep; vhep = 0;
393  delete [] scale; scale = 0;
394  delete [] spin; spin = 0;
395  delete [] colorflow; colorflow = 0;
396  blockid = ntot = nevhep = nhep = 0;
397  isEmpty = 1;
398  return;
399 }
400 
401 long lStdHep::Event::read(lStdHep &ls)
402 {
403 //
404 // Read event header
405 //
406  long len;
407 
408  cleanup();
409 
410  blockid = ls.readLong();
411  ntot = ls.readLong();
412  version = ls.readString(len);
413  if (blockid != LSH_EVENTHEADER) ls.setError(LSH_NOEVENT);
414 
415  evtnum = ls.readLong();
416  storenum = ls.readLong();
417  runnum = ls.readLong();
418  trigMask = ls.readLong();
419  nBlocks = ls.readLong();
420  dimBlocks = ls.readLong();
421 
422  if (*version == '2') {
423  nNTuples = ls.readLong();
424  dimNTuples = ls.readLong();
425  if (dimBlocks) {
426  blockIds = ls.readLongArray(len);
427  ptrBlocks = ls.readLongArray(len);
428  }
429  if (dimNTuples) {
430  ls.setError(LSH_NOTSUPPORTED);
431  return(ls.getError());
432  }
433  }
434  else {
435  nNTuples = 0;
436  dimNTuples = 0;
437  blockIds = ls.readLongArray(len);
438  ptrBlocks = ls.readLongArray(len);
439  }
440 //
441 // Read event
442 //
443  for (int i = 0; i < nBlocks; i++) {
444  blockid = ls.readLong();
445  ntot = ls.readLong();
446  if (version) delete [] version;
447  version = ls.readString(len);
448 
449  isEmpty = 0;
450  switch (blockIds[i]) {
451  case LSH_STDHEP : // 101
452  nevhep = ls.readLong();
453  nhep = ls.readLong();
454  if (isthep) delete [] isthep;
455  isthep = ls.readLongArray(len);
456  if (idhep) delete [] idhep;
457  idhep = ls.readLongArray(len);
458  if (jmohep) delete [] jmohep;
459  jmohep = ls.readLongArray(len);
460  if (jdahep) delete [] jdahep;
461  jdahep = ls.readLongArray(len);
462  if (phep) delete [] phep;
463  phep = ls.readDoubleArray(len);
464  if (vhep) delete [] vhep;
465  vhep = ls.readDoubleArray(len);
466  break;
467  case LSH_STDHEPEV4 : // 201
468  nevhep = ls.readLong();
469  nhep = ls.readLong();
470  if (isthep) delete [] isthep;
471  isthep = ls.readLongArray(len);
472  if (idhep) delete [] idhep;
473  idhep = ls.readLongArray(len);
474  if (jmohep) delete [] jmohep;
475  jmohep = ls.readLongArray(len);
476  if (jdahep) delete [] jdahep;
477  jdahep = ls.readLongArray(len);
478  if (phep) delete [] phep;
479  phep = ls.readDoubleArray(len);
480  if (vhep) delete [] vhep;
481  vhep = ls.readDoubleArray(len);
482 //
483 // New stuff for STDHEPEV4:
484 //
485  eventweight = ls.readDouble();
486  alphaqed = ls.readDouble();
487  alphaqcd = ls.readDouble();
488  if (scale) delete [] scale;
489  scale = ls.readDoubleArray(len);
490  if (spin) delete [] spin;
491  spin = ls.readDoubleArray(len);
492  if (colorflow) delete [] colorflow;
493  colorflow = ls.readLongArray(len);
494  idrup = ls.readLong();
495  break;
496  case LSH_OFFTRACKARRAYS : // 102
497  case LSH_OFFTRACKSTRUCT : // 103
498  case LSH_TRACEARRAYS : // 104
499  case LSH_STDHEPM : // 105
500  break;
501  case LSH_STDHEPBEG : // 106
502  bnevtreq = ls.readLong();
503  bnevtgen = ls.readLong();
504  bnevtwrt = ls.readLong();
505  bstdecom = ls.readFloat();
506  bstdxsec = ls.readFloat();
507  bstdseed1 = ls.readDouble();
508  bstdseed2 = ls.readDouble();
509  isEmpty = 1;
510  break;
511  case LSH_STDHEPEND : // 107
512  enevtreq = ls.readLong();
513  enevtgen = ls.readLong();
514  enevtwrt = ls.readLong();
515  estdecom = ls.readFloat();
516  estdxsec = ls.readFloat();
517  estdseed1 = ls.readDouble();
518  estdseed2 = ls.readDouble();
519  isEmpty = 1;
520  break;
521  case LSH_STDHEPCXX : // 108
522  break;
523  }
524  }
525  return(ls.getError());
526 }
527 
528 long lStdHep::Event::printHeader(FILE *fp)
529 {
530  fprintf(fp, " EventHeader: blockid: %ld, ntot: %ld, version: %s\n", blockid, ntot, version);
531  fprintf(fp, " : evtnum: %ld, storenum: %ld, runnum: %ld, trigMask: %ld, nBlocks: %ld, dimBlocks: %ld\n",
532  evtnum, storenum, runnum, trigMask, nBlocks, dimBlocks);
533  fprintf(fp, " : nNTuples: %ld, dimNTuples: %ld\n", nNTuples, dimNTuples);
534 
535  for (int i = 0; i < nBlocks; i++) {
536  const char *labels[10] = {"Event", "Off-track arrays", "Off-track struct", "Trace arrays",
537  "Event with multiple interactions", "Begin run", "End run", "StdHepCXX",
538  "EventV4", "Unknown" };
539  int j = (int)blockIds[i] - 101;
540  if (blockIds[i] == LSH_STDHEPEV4) j = 8;
541  if ((j < 0) || (j > 9)) j = 9;
542  fprintf(fp, " : %d: blockIds %ld (%s) ptrBlocks %ld\n",
543  i, blockIds[i], labels[j], ptrBlocks[i]);
544  }
545  return(0);
546 }
547 
548 long lStdHep::Event::print(FILE *fp)
549 {
550  fprintf(fp, " Event: nevhep: %ld, nhep: %ld\n", nevhep, nhep);
551  return(0);
552 }
553 
554 }
UTIL
Definition: lStdHep.cc:17
Generated on Wed Jul 6 2016 11:13:16 for GEMC by   doxygen 1.8.11