ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Fun4AllHepMCInputManager.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Fun4AllHepMCInputManager.cc
2 
3 #include "PHHepMCGenEvent.h"
4 #include "PHHepMCGenEventMap.h"
5 
6 #include <frog/FROG.h>
7 
8 #include <fun4all/Fun4AllInputManager.h> // for Fun4AllInpu...
10 #include <fun4all/Fun4AllServer.h>
12 
13 #include <phool/PHCompositeNode.h>
14 #include <phool/PHIODataNode.h> // for PHIODataNode
15 #include <phool/PHNodeIterator.h> // for PHNodeIterator
16 #include <phool/PHObject.h> // for PHObject
17 #include <phool/getClass.h>
18 #include <phool/phool.h> // for PHWHERE
19 #include <phool/recoConsts.h>
20 
21 #include <HepMC/GenEvent.h>
22 #include <HepMC/GenParticle.h> // for GenParticle
23 #include <HepMC/GenVertex.h> // for GenVertex
24 #include <HepMC/IO_GenEvent.h>
25 #include <HepMC/SimpleVector.h> // for FourVector
26 #include <HepMC/Units.h> // for CM, GEV
27 
28 #include <TDirectory.h>
29 #include <TPRegexp.h>
30 #include <TROOT.h>
31 #include <TString.h>
32 #include <TSystem.h>
33 
34 #include <boost/iostreams/filter/bzip2.hpp>
35 #include <boost/iostreams/filter/gzip.hpp>
36 
37 #include <cmath>
38 #include <cstdlib>
39 #include <fstream>
40 #include <iostream>
41 #include <map> // for _Rb_tree_it...
42 #include <sstream>
43 #include <vector> // for vector
44 
45 static const double toMM = 1.e-12;
46 
47 Fun4AllHepMCInputManager::Fun4AllHepMCInputManager(const std::string &name, const std::string &nodename, const std::string &topnodename)
48  : Fun4AllInputManager(name, nodename, topnodename)
49  , topNodeName(topnodename)
50 {
51  set_embedding_id(0); // default embedding ID. Welcome to change via macro
52 
55  PHNodeIterator iter(topNode);
56  PHCompositeNode *dstNode = se->getNode(InputNode(), topNodeName);
57 
58  PHHepMCGenEventMap *geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
59  if (!geneventmap)
60  {
61  geneventmap = new PHHepMCGenEventMap();
62  PHIODataNode<PHObject> *newmapnode = new PHIODataNode<PHObject>(geneventmap, "PHHepMCGenEventMap", "PHObject");
63  dstNode->addNode(newmapnode);
64  }
65 
67 
68  return;
69 }
70 
72 {
73  fileclose();
74  if (!m_HepMCTmpFile.empty())
75  {
76  // okay if the file does not exist
77  remove(m_HepMCTmpFile.c_str());
78  }
79  delete ascii_in;
80  delete filestream;
81  delete unzipstream;
82 }
83 
84 int Fun4AllHepMCInputManager::fileopen(const std::string &filenam)
85 {
86  if (!MySyncManager())
87  {
88  std::cout << "Call fileopen only after you registered your Input Manager " << Name() << " with the Fun4AllServer" << std::endl;
89  exit(1);
90  }
91  if (IsOpen())
92  {
93  std::cout << "Closing currently open file "
94  << filename
95  << " and opening " << filenam << std::endl;
96  fileclose();
97  }
98  filename = filenam;
99  FROG frog;
100  std::string fname(frog.location(filename));
101  if (Verbosity() > 0)
102  {
103  std::cout << Name() << ": opening file " << fname << std::endl;
104  }
105 
106  if (m_ReadOscarFlag)
107  {
108  theOscarFile.open(fname);
109  }
110  else
111  {
112  TString tstr(fname);
113  TPRegexp bzip_ext(".bz2$");
114  TPRegexp gzip_ext(".gz$");
115  if (tstr.Contains(bzip_ext))
116  {
117  // use boost iosteam library to decompress bz2 on the fly
118  filestream = new std::ifstream(fname, std::ios::in | std::ios::binary);
119  zinbuffer.push(boost::iostreams::bzip2_decompressor());
120  zinbuffer.push(*filestream);
121  unzipstream = new std::istream(&zinbuffer);
122  ascii_in = new HepMC::IO_GenEvent(*unzipstream);
123  }
124  else if (tstr.Contains(gzip_ext))
125  {
126  // use boost iosream to decompress the gzip file on the fly
127  filestream = new std::ifstream(fname, std::ios::in | std::ios::binary);
128  zinbuffer.push(boost::iostreams::gzip_decompressor());
129  zinbuffer.push(*filestream);
130  unzipstream = new std::istream(&zinbuffer);
131  ascii_in = new HepMC::IO_GenEvent(*unzipstream);
132  }
133  else
134  {
135  // expects normal ascii hepmc file
136  ascii_in = new HepMC::IO_GenEvent(fname, std::ios::in);
137  }
138  }
139 
141  static bool run_number_forced = rc->FlagExist("RUNNUMBER");
142  if (run_number_forced)
143  {
144  MySyncManager()->CurrentRun(rc->get_IntFlag("RUNNUMBER"));
145  }
146  else
147  {
148  MySyncManager()->CurrentRun(-1);
149  }
150  events_thisfile = 0;
151  IsOpen(1);
152  AddToFileOpened(fname); // add file to the list of files which were opened
153  return 0;
154 }
155 
156 int Fun4AllHepMCInputManager::run(const int /*nevents*/)
157 {
158  // attempt to retrieve a valid event from inputs
159  while (true)
160  {
161  if (!IsOpen())
162  {
163  if (FileListEmpty())
164 
165  {
166  if (Verbosity() > 0)
167  {
168  std::cout << "Fun4AllHepMCInputManager::run::" << Name() << ": No Input file open" << std::endl;
169  }
170  return -1;
171  }
172  else
173  {
174  if (OpenNextFile())
175  {
176  std::cout << "Fun4AllHepMCInputManager::run::" << Name() << ": No Input file from filelist opened" << std::endl;
177  return -1;
178  }
179  }
180  }
181 
182  if (m_EventPushedBackFlag) // if an event was pushed back, reuse save copy
183  {
184  HepMC::IO_GenEvent ascii_tmp_in(m_HepMCTmpFile, std::ios::in);
185  ascii_tmp_in >> evt;
187  remove(m_HepMCTmpFile.c_str());
188  }
189  else
190  {
191  if (m_ReadOscarFlag)
192  {
193  evt = ConvertFromOscar();
194  }
195  else
196  {
197  evt = ascii_in->read_next_event();
198  }
199  }
200 
201  if (!evt)
202  {
203  if (Verbosity() > 1)
204  {
205  std::cout << "Fun4AllHepMCInputManager::run::" << Name()
206  << ": error type: " << ascii_in->error_type()
207  << ", rdstate: " << ascii_in->rdstate() << std::endl;
208  }
209  fileclose();
210  }
211  else
212  {
213  MySyncManager()->CurrentEvent(evt->event_number());
214  if (Verbosity() > 0)
215  {
216  std::cout << "Fun4AllHepMCInputManager::run::" << Name()
217  << ": hepmc evt no: " << evt->event_number() << std::endl;
218  }
219  m_MyEvent.push_back(evt->event_number());
221  if (ievt != PHHepMCGenHelper::get_geneventmap()->end())
222  {
223  // override existing event
224  ievt->second->addEvent(evt);
225  }
226  else
227  {
229  }
230 
231  events_total++;
232  events_thisfile++;
233 
234  // check if the local SubsysReco discards this event
236  {
237  // if this event is discarded we only need to remove the event from the list of event numbers
238  // the new event will overwrite the event on the node tree without issues
239  m_MyEvent.pop_back();
240  }
241  else
242  {
243  break; // have a good event, move on
244  }
245  }
246  } // attempt to retrieve a valid event from inputs
248 }
249 
251 {
252  if (!IsOpen())
253  {
254  std::cout << Name() << ": fileclose: No Input file open" << std::endl;
255  return -1;
256  }
257  if (m_ReadOscarFlag)
258  {
259  theOscarFile.close();
260  }
261  else
262  {
263  delete ascii_in;
264  ascii_in = nullptr;
265  }
266  IsOpen(0);
267  // if we have a file list, move next entry to top of the list
268  // or repeat the same entry again
269  UpdateFileList();
270  return 0;
271 }
272 
273 void Fun4AllHepMCInputManager::Print(const std::string &what) const
274 {
276  std::cout << Name() << " Vertex Settings: " << std::endl;
278  return;
279 }
280 
282 {
283  // PushBackEvents is supposedly pushing events back on the stack which works
284  // easily with root trees (just grab a different entry) but hard in these HepMC ASCII files.
285  // A special case is when the synchronization fails and we need to only push back a single
286  // event. In this case we save the evt in a temporary file which is read back in the run method
287  // instead of getting the next event.
288  if (i > 0)
289  {
290  if (i == 1 && evt) // check on evt pointer makes sure it is not done from the cmd line
291  {
292  // root barfs when writing the node to the output.
293  // Saving the pointer - even using a deep copy and reusing it did not work
294  // The hackaround which works is to write this event into a temporary file and read it back
295  if (m_HepMCTmpFile.empty())
296  {
297  // we need to create this filename just once, we reuse it. Do it only if we need it
298  m_HepMCTmpFile = "/tmp/HepMCTmpEvent-" + Name() + "-" + std::to_string(getpid()) + ".hepmc";
299  }
300  HepMC::IO_GenEvent ascii_io(m_HepMCTmpFile, std::ios::out);
301  ascii_io << evt;
303  m_MyEvent.pop_back();
304 
305  if (Verbosity() > 3)
306  {
307  std::cout << Name() << ": pushing back evt no " << evt->event_number() << std::endl;
308  }
309  return 0;
310  }
311  std::cout << PHWHERE << Name()
312  << " Fun4AllHepMCInputManager cannot pop back more than 1 event"
313  << std::endl;
314  return -1;
315  }
316  if (!IsOpen())
317  {
318  std::cout << PHWHERE << Name()
319  << " no file opened yet" << std::endl;
320  return -1;
321  }
322  // Skipping events is implemented as
323  // pushing a negative number of events on the stack, so in order to implement
324  // the skipping of events we read -i events.
325  int nevents = -i; // negative number of events to push back -> skip num events
326  int errorflag = 0;
327  while (nevents > 0 && !errorflag)
328  {
329  evt = ascii_in->read_next_event();
330  if (!evt)
331  {
332  std::cout << "Error after skipping " << i - nevents << std::endl;
333  std::cout << "error type: " << ascii_in->error_type()
334  << ", rdstate: " << ascii_in->rdstate() << std::endl;
335  errorflag = -1;
336  fileclose();
337  }
338  else
339  {
340  m_MyEvent.push_back(evt->event_number());
341  if (Verbosity() > 3)
342  {
343  std::cout << "Skipping evt no: " << evt->event_number() << std::endl;
344  }
345  }
346  delete evt;
347  nevents--;
348  }
349  return errorflag;
350 }
351 
352 HepMC::GenEvent *
354 {
355  if (theOscarFile.eof()) // if the file is exhausted bail out during this next read
356  {
357  return nullptr;
358  }
359 
360  delete evt;
361  //use PHENIX unit
362  evt = new HepMC::GenEvent(HepMC::Units::GEV, HepMC::Units::CM);
363 
364  if (Verbosity() > 1) std::cout << "Reading Oscar Event " << events_total << std::endl;
365  //Grab New Event From Oscar
366  std::string theLine;
367  std::vector<std::vector<double> > theEventVec;
368  std::vector<HepMC::FourVector> theVtxVec;
369  while (getline(theOscarFile, theLine))
370  {
371  if (theLine.compare(0, 1, "#") == 0) continue;
372  std::vector<double> theInfo; //format: N,pid,px,py,pz,E,mass,xvtx,yvtx,zvtx,?
373  double number = NAN;
374  for (std::istringstream numbers_iss(theLine); numbers_iss >> number;)
375  {
376  theInfo.push_back(number);
377  }
378 
379  if (theInfo.size() == 2 && theInfo[0] == 0 && theInfo[1] == 0)
380  {
381  break;
382  }
383  else if (theInfo.size() == 2 && theInfo[0] == 0 && theInfo[1] > 0)
384  {
385  continue;
386  }
387  else
388  {
389  theEventVec.push_back(theInfo);
390  HepMC::FourVector vert(theInfo[8] * toMM, theInfo[9] * toMM, theInfo[10] * toMM, theInfo[11]);
391  theVtxVec.push_back(vert);
392  }
393 
394  } //while(getline)
395 
396  //Set Event Number
397  evt->set_event_number(events_total);
398 
399  //Loop Over One Event, Fill HepMC
400  for (unsigned int i = 0; i < theEventVec.size(); i++)
401  {
402  //int N = (int)theEventVec[i][0];
403  int pid = (int) theEventVec[i][1];
404  double px = theEventVec[i][3];
405  double py = theEventVec[i][4];
406  double pz = theEventVec[i][5];
407  double E = theEventVec[i][6];
408  double m = theEventVec[i][7];
409  int status = 1; //oscar only writes final state particles
410 
411  HepMC::GenVertex *v = new HepMC::GenVertex(theVtxVec[i]);
412  evt->add_vertex(v);
413 
414  HepMC::GenParticle *p = new HepMC::GenParticle(HepMC::FourVector(px, py, pz, E), pid, status);
415  p->setGeneratedMass(m);
416  p->suggest_barcode(i + 1);
417  v->add_particle_out(p);
418  }
419  if (Verbosity() > 3)
420  {
421  evt->print();
422  }
423  return evt;
424 }
425 
427 {
428  m_MyEvent.clear();
429  return 0;
430 }
431 
432 int Fun4AllHepMCInputManager::MyCurrentEvent(const unsigned int index) const
433 {
434  if (m_MyEvent.empty())
435  {
436  return 0;
437  }
438  return m_MyEvent.at(index);
439 }
440