Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Fun4AllOscarInputManager.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Fun4AllOscarInputManager.cc
2 
3 #include "PHHepMCGenEvent.h" // for PHHepMCGenE...
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>
17 #include <phool/getClass.h>
18 #include <phool/recoConsts.h>
19 
20 #include <HepMC/GenEvent.h>
21 #include <HepMC/GenParticle.h> // for GenParticle
22 #include <HepMC/GenVertex.h> // for GenVertex
23 #include <HepMC/SimpleVector.h> // for FourVector
24 #include <HepMC/Units.h> // for GEV, MM
25 
26 #include <TPRegexp.h>
27 #include <TString.h>
28 
29 #include <boost/iostreams/filter/bzip2.hpp>
30 #include <boost/iostreams/filter/gzip.hpp>
31 #include <boost/iostreams/filtering_streambuf.hpp>
32 
33 #include <cstdlib>
34 #include <fstream>
35 #include <iostream>
36 #include <map> // for _Rb_tree_it...
37 #include <sstream>
38 #include <utility> // for swap, pair
39 #include <vector> // for vector
40 
41 using namespace std;
42 
43 static boost::iostreams::filtering_streambuf<boost::iostreams::input> zinbuffer;
44 static const double toMM = 1.e-12;
46 
47 Fun4AllOscarInputManager::Fun4AllOscarInputManager(const string &name, const string &topnodename)
48  : Fun4AllInputManager(name, "")
49  , events_total(0)
50  , events_thisfile(0)
51  , topNodeName(topnodename)
52  , evt(nullptr)
53  , skipEvents(0)
54  , skippedEvents(0)
55  , filestream(nullptr)
56  , unzipstream(nullptr)
57  , isCompressed(false)
58 {
60  topNode = se->topNode(topNodeName.c_str());
61  PHNodeIterator iter(topNode);
62  PHCompositeNode *dstNode = se->getNode(InputNode(), topNodeName);
63 
64  PHHepMCGenEventMap *geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
65  if (!geneventmap)
66  {
67  geneventmap = new PHHepMCGenEventMap();
68  PHIODataNode<PHObject> *newmapnode = new PHIODataNode<PHObject>(geneventmap, "PHHepMCGenEventMap", "PHObject");
69  dstNode->addNode(newmapnode);
70  }
71 
73 }
74 
76 {
77  fileclose();
78  delete filestream;
79  delete unzipstream;
80 }
81 
82 int Fun4AllOscarInputManager::fileopen(const string &filenam)
83 {
84  if (!MySyncManager())
85  {
86  cout << "Call fileopen only after you registered your Input Manager " << Name() << " with the Fun4AllServer" << endl;
87  exit(1);
88  }
89  if (IsOpen())
90  {
91  cout << "Closing currently open file "
92  << filename
93  << " and opening " << filenam << endl;
94  fileclose();
95  }
96  filename = filenam;
97  FROG frog;
98  string fname(frog.location(filename.c_str()));
99  if (Verbosity() > 0)
100  {
101  cout << Name() << ": opening file " << fname << endl;
102  }
103 
104  TString tstr(fname);
105  TPRegexp bzip_ext(".bz2$");
106  TPRegexp gzip_ext(".gz$");
107  if (tstr.Contains(bzip_ext))
108  {
109  // use boost iosteam library to decompress bz2 on the fly
110  filestream = new ifstream(fname.c_str(), std::ios::in | std::ios::binary);
111  zinbuffer.push(boost::iostreams::bzip2_decompressor());
112  zinbuffer.push(*filestream);
113  unzipstream = new istream(&zinbuffer);
114  isCompressed = true;
115  }
116  else if (tstr.Contains(gzip_ext))
117  {
118  // use boost iosream to decompress the gzip file on the fly
119  filestream = new ifstream(fname.c_str(), std::ios::in | std::ios::binary);
120  zinbuffer.push(boost::iostreams::gzip_decompressor());
121  zinbuffer.push(*filestream);
122  unzipstream = new istream(&zinbuffer);
123  isCompressed = true;
124  }
125  else
126  {
127  theOscarFile.open(fname.c_str());
128  }
129 
131  static bool run_number_forced = rc->FlagExist("RUNNUMBER");
132  if (run_number_forced)
133  {
134  MySyncManager()->CurrentRun(rc->get_IntFlag("RUNNUMBER"));
135  }
136  else
137  {
138  MySyncManager()->CurrentRun(-1);
139  }
140  events_thisfile = 0;
141  IsOpen(1);
142  AddToFileOpened(fname); // add file to the list of files which were opened
143  return 0;
144 }
145 
147 {
148 readagain:
149  if (!IsOpen())
150  {
151  if (FileListEmpty())
152  {
153  if (Verbosity() > 0)
154  {
155  cout << Name() << ": No Input file open" << endl;
156  }
157  return 0;
158  }
159  else
160  {
161  if (OpenNextFile())
162  {
163  cout << Name() << ": No Input file from filelist opened" << endl;
164  return 0;
165  }
166  }
167  }
168 
169  //Read oscar
170  evt = new HepMC::GenEvent(HepMC::Units::GEV, HepMC::Units::MM);
171  int code = ConvertFromOscar();
172  /*
173  if(skippedEvents < skipEvents || code == 3)
174  {
175  goto readagain;
176  }
177  */
178  if (code == 1 || evt == nullptr)
179  {
180  if (Verbosity() > 1) cout << "Finished file!" << endl;
181  fileclose();
182  goto readagain;
183  }
184 
185  // if(Verbosity() > 4) cout << "SIZE: " << phhepmcgenevt->size() << endl;
186  //MySyncManager()->CurrentEvent(evt->event_number());
187  events_total++;
188  events_thisfile++;
189 
190  // check if the local SubsysReco discards this event
192  {
193  //ResetEvent();
194  goto readagain;
195  }
196 
197  if (events_total < nevents)
198  {
199  //ResetEvent();
200  goto readagain;
201  }
202 
203  return 0;
204 }
205 
207 {
208  if (!IsOpen())
209  {
210  cout << Name() << ": fileclose: No Input file open" << endl;
211  return -1;
212  }
213  if (isCompressed)
214  {
215  filestream->close();
216  }
217  else
218  {
219  theOscarFile.close();
220  }
221  IsOpen(0);
222  // if we have a file list, move next entry to top of the list
223  // or repeat the same entry again
224  UpdateFileList();
225  return 0;
226 }
227 
228 void Fun4AllOscarInputManager::Print(const string &what) const
229 {
231  return;
232 }
233 
235 {
236  //delete evt;
237  //evt = nullptr;
238  return 0;
239 }
240 
242 {
243  int counter = 0;
244 
245  while (counter < i)
246  {
247  std::string theLine;
248  while (getline(theOscarFile, theLine))
249  {
250  if (theLine.compare(0, 1, "#") == 0) continue;
251  vector<double> theInfo;
252  double number = NAN;
253  for (istringstream numbers_iss(theLine); numbers_iss >> number;)
254  {
255  theInfo.push_back(number);
256  }
257 
258  if (theInfo.size() == 2 && theInfo[0] == 0 && theInfo[1] == 0)
259  {
260  counter++;
261  skippedEvents++;
262  break;
263  }
264  else if (theInfo.size() == 2 && theInfo[0] == 0 && theInfo[1] > 0)
265  {
266  continue;
267  }
268  }
269  }
270 
271  if (theOscarFile.eof()) return -1;
272 
273  return 0;
274 }
275 
277 {
278  delete evt;
279  evt = nullptr;
280 
281  if (theOscarFile.eof()) // if the file is exhausted bail out during this next read
282  {
283  cout << "Oscar EOF" << endl;
284  return 1;
285  }
286  evt = new HepMC::GenEvent(HepMC::Units::GEV, HepMC::Units::MM);
287 
288  if (Verbosity() > 1) cout << "Reading Oscar Event " << events_total + skippedEvents + 1 << endl;
289  //Grab New Event From Oscar
290  string theLine;
291  vector<vector<double> > theEventVec;
292  // vector<HepMC::FourVector> theVtxVec;
293  if (isCompressed)
294  {
295  // while(getline(unzipstream, theLine))
296  // {
297  // if(theLine.find("#") == 0) continue;
298  // vector<double> theInfo; //format: N,pid,px,py,pz,E,mass,xvtx,yvtx,zvtx,?
299  // double number;
300  // for(istringstream numbers_iss(theLine); numbers_iss >> number; )
301  // {
302  // theInfo.push_back(number);
303  // }
304 
305  // if(theInfo.size() == 2 && theInfo[0] == 0 && theInfo[1] == 0)
306  // {
307  // break;
308  // }
309  // else if (theInfo.size() == 2 && theInfo[0] == 0 && theInfo[1] > 0)
310  // {
311  // continue;
312  // }
313  // else
314  // {
315  // theEventVec.push_back(theInfo);
316  // HepMC::FourVector vert(theInfo[8]*toMM, theInfo[9]*toMM, theInfo[10]*toMM, theInfo[11]);
317  // theVtxVec.push_back(vert);
318  // }
319 
320  // }//while(getline)
321  }
322  else
323  {
324  while (getline(theOscarFile, theLine))
325  {
326  if (theLine.compare(0, 1, "#") == 0) continue;
327  vector<double> theInfo; //format: N,pid,px,py,pz,E,mass,xvtx,yvtx,zvtx,?
328  double number = NAN;
329  for (istringstream numbers_iss(theLine); numbers_iss >> number;)
330  {
331  theInfo.push_back(number);
332  }
333 
334  if (theInfo.size() == 2 && theInfo[0] == 0 && theInfo[1] == 0)
335  {
336  break;
337  }
338  else if (theInfo.size() == 2 && theInfo[0] == 0 && theInfo[1] > 0)
339  {
340  continue;
341  }
342  else
343  {
344  theEventVec.push_back(theInfo);
345  }
346 
347  } //while(getline)
348  }
349 
350  /*
351  if(skippedEvents < skipEvents)
352  {
353  skippedEvents++;
354  if (Verbosity() > 5) cout << "Skipping event " << skippedEvents << endl;
355  return 2;
356  }
357  */
358 
359  //Set Event Number
360  evt->set_event_number(events_total + 1);
361 
362  //Loop Over One Event, Fill particles
363  std::vector<HepMC::GenParticle *> hepevt_particles(theEventVec.size());
364  for (unsigned int i = 0; i < theEventVec.size(); i++)
365  {
366  //int N = (int)theEventVec[i][0];
367  int pid = (int) theEventVec[i][1];
368  double px = theEventVec[i][3];
369  double py = theEventVec[i][4];
370  double pz = theEventVec[i][5];
371  double E = theEventVec[i][6];
372  double m = theEventVec[i][7];
373  int status = 1; //oscar only writes final state particles
374 
375  hepevt_particles[i] = new HepMC::GenParticle(HepMC::FourVector(px, py, pz, E), pid, status);
376  hepevt_particles[i]->setGeneratedMass(m);
377  hepevt_particles[i]->suggest_barcode(i + 1);
378  }
379 
380  for (unsigned int i = 0; i < theEventVec.size(); i++)
381  {
382  HepMC::GenParticle *p = hepevt_particles[i];
383  HepMC::GenVertex *prod_vtx = p->production_vertex();
384  if (prod_vtx) prod_vtx->add_particle_out(p);
385 
386  bool found = false;
387  HepMC::FourVector prod_pos(theEventVec[i][8] * toMM, theEventVec[i][9] * toMM, theEventVec[i][10] * toMM, theEventVec[i][11]);
388  if (!prod_vtx)
389  {
390  //See if the vertex is already in the event
391  for (HepMC::GenEvent::vertex_iterator v = evt->vertices_begin(); v != evt->vertices_end(); ++v)
392  {
393  HepMC::GenVertex *theV = *v;
394  if (theV->position().x() != prod_pos.x()) continue;
395  if (theV->position().y() != prod_pos.y()) continue;
396  if (theV->position().z() != prod_pos.z()) continue;
397  found = true;
398  theV->add_particle_out(p);
399  }
400  if (!found)
401  {
402  //Didn't find vertex, add it
403  prod_vtx = new HepMC::GenVertex(prod_pos);
404  prod_vtx->add_particle_out(p);
405  evt->add_vertex(prod_vtx);
406  }
407  }
408 
409  // If prod_vtx doesn't already have position specified, fill it.
410  if (!found && prod_vtx && prod_vtx->position() == HepMC::FourVector()) prod_vtx->set_position(prod_pos);
411  }
412 
413  evt->print();
414  if (Verbosity() > 5) evt->print();
415  if (Verbosity() > 3) cout << "Adding Event to phhepmcgenevt" << endl;
416 
419  if (ievt != PHHepMCGenHelper::get_geneventmap()->end())
420  {
421  // override existing event
422  ievt->second->addEvent(evt);
423  }
424  else
426  return 0;
427 }