generators/phhepmc/Fun4AllHepMCInputManager.cc

0001 #include "Fun4AllHepMCInputManager.h"
0002
0003 #include "PHHepMCGenEvent.h"
0004 #include "PHHepMCGenEventMap.h"
0005
0006 #include <frog/FROG.h>
0007
0008 #include <fun4all/Fun4AllInputManager.h>  // for Fun4AllInpu...
0009 #include <fun4all/Fun4AllReturnCodes.h>
0010 #include <fun4all/Fun4AllServer.h>
0011 #include <fun4all/Fun4AllSyncManager.h>
0012
0013 #include <phool/PHCompositeNode.h>
0014 #include <phool/PHIODataNode.h>    // for PHIODataNode
0015 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0016 #include <phool/PHObject.h>        // for PHObject
0017 #include <phool/getClass.h>
0018 #include <phool/phool.h>  // for PHWHERE
0019 #include <phool/recoConsts.h>
0020
0021 #include <HepMC/GenEvent.h>
0022 #include <HepMC/GenParticle.h>  // for GenParticle
0023 #include <HepMC/GenVertex.h>    // for GenVertex
0024 #include <HepMC/IO_GenEvent.h>
0025 #include <HepMC/SimpleVector.h>  // for FourVector
0026 #include <HepMC/Units.h>         // for CM, GEV
0027
0028 #include <TPRegexp.h>
0029 #include <TString.h>
0030
0031 #include <boost/iostreams/filter/bzip2.hpp>
0032 #include <boost/iostreams/filter/gzip.hpp>
0033
0034 #include <cmath>
0035 #include <cstdlib>
0036 #include <fstream>
0037 #include <iostream>
0038 #include <map>  // for _Rb_tree_it...
0039 #include <sstream>
0040 #include <vector>  // for vector
0041
0042 static const double toMM = 1.e-12;
0043
0044 Fun4AllHepMCInputManager::Fun4AllHepMCInputManager(const std::string &name, const std::string &nodename, const std::string &topnodename)
0045   : Fun4AllInputManager(name, nodename, topnodename)
0046   , topNodeName(topnodename)
0047 {
0048   Fun4AllServer *se = Fun4AllServer::instance();
0049   topNode = se->topNode(topNodeName);
0050   PHNodeIterator iter(topNode);
0051   PHCompositeNode *dstNode = se->getNode(InputNode(), topNodeName);
0052
0053   PHHepMCGenEventMap *geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0054   if (!geneventmap)
0055   {
0056     geneventmap = new PHHepMCGenEventMap();
0057     PHIODataNode<PHObject> *newmapnode = new PHIODataNode<PHObject>(geneventmap, "PHHepMCGenEventMap", "PHObject");
0058     dstNode->addNode(newmapnode);
0059   }
0060
0061   PHHepMCGenHelper::set_geneventmap(geneventmap);
0062
0063   return;
0064 }
0065
0066 Fun4AllHepMCInputManager::~Fun4AllHepMCInputManager()
0067 {
0068   fileclose();
0069   if (!m_HepMCTmpFile.empty())
0070   {
0071     // okay if the file does not exist
0072     remove(m_HepMCTmpFile.c_str());
0073   }
0074   delete ascii_in;
0075   delete filestream;
0076   delete unzipstream;
0077 }
0078
0079 int Fun4AllHepMCInputManager::fileopen(const std::string &filenam)
0080 {
0081   if (!MySyncManager())
0082   {
0083     std::cout << "Call fileopen only after you registered your Input Manager " << Name() << " with the Fun4AllServer" << std::endl;
0084     exit(1);
0085   }
0086   if (IsOpen())
0087   {
0088     std::cout << "Closing currently open file "
0089               << filename
0090               << " and opening " << filenam << std::endl;
0091     fileclose();
0092   }
0093   filename = filenam;
0094   FROG frog;
0095   std::string fname(frog.location(filename));
0096   if (Verbosity() > 0)
0097   {
0098     std::cout << Name() << ": opening file " << fname << std::endl;
0099   }
0100
0101   if (m_ReadOscarFlag)
0102   {
0103     theOscarFile.open(fname);
0104   }
0105   else
0106   {
0107     TString tstr(fname);
0108     TPRegexp bzip_ext(".bz2$");
0109     TPRegexp gzip_ext(".gz$");
0110     if (tstr.Contains(bzip_ext))
0111     {
0112       // use boost iosteam library to decompress bz2 on the fly
0113       filestream = new std::ifstream(fname, std::ios::in | std::ios::binary);
0114       zinbuffer.push(boost::iostreams::bzip2_decompressor());
0115       zinbuffer.push(*filestream);
0116       unzipstream = new std::istream(&zinbuffer);
0117       ascii_in = new HepMC::IO_GenEvent(*unzipstream);
0118     }
0119     else if (tstr.Contains(gzip_ext))
0120     {
0121       // use boost iosream to decompress the gzip file on the fly
0122       filestream = new std::ifstream(fname, std::ios::in | std::ios::binary);
0123       zinbuffer.push(boost::iostreams::gzip_decompressor());
0124       zinbuffer.push(*filestream);
0125       unzipstream = new std::istream(&zinbuffer);
0126       ascii_in = new HepMC::IO_GenEvent(*unzipstream);
0127     }
0128     else
0129     {
0130       // expects normal ascii hepmc file
0131       ascii_in = new HepMC::IO_GenEvent(fname, std::ios::in);
0132     }
0133   }
0134
0135   recoConsts *rc = recoConsts::instance();
0136   static bool run_number_forced = rc->FlagExist("RUNNUMBER");
0137   if (run_number_forced)
0138   {
0139     MySyncManager()->CurrentRun(rc->get_IntFlag("RUNNUMBER"));
0140   }
0141   else
0142   {
0143     MySyncManager()->CurrentRun(-1);
0144   }
0145   events_thisfile = 0;
0146   IsOpen(1);
0147   AddToFileOpened(fname);  // add file to the list of files which were opened
0148   return 0;
0149 }
0150
0151 int Fun4AllHepMCInputManager::run(const int /*nevents*/)
0152 {
0153   // attempt to retrieve a valid event from inputs
0154   while (true)
0155   {
0156     if (!IsOpen())
0157     {
0158       if (FileListEmpty())
0159
0160       {
0161         if (Verbosity() > 0)
0162         {
0163           std::cout << "Fun4AllHepMCInputManager::run::" << Name() << ": No Input file open" << std::endl;
0164         }
0165         return -1;
0166       }
0167       else
0168       {
0169         if (OpenNextFile())
0170         {
0171           std::cout << "Fun4AllHepMCInputManager::run::" << Name() << ": No Input file from filelist opened" << std::endl;
0172           return -1;
0173         }
0174       }
0175     }
0176
0177     if (m_EventPushedBackFlag)  // if an event was pushed back, reuse save copy
0178     {
0179       HepMC::IO_GenEvent ascii_tmp_in(m_HepMCTmpFile, std::ios::in);
0180       ascii_tmp_in >> evt;
0181       m_EventPushedBackFlag = 0;
0182       remove(m_HepMCTmpFile.c_str());
0183     }
0184     else
0185     {
0186       if (m_ReadOscarFlag)
0187       {
0188         evt = ConvertFromOscar();
0189       }
0190       else
0191       {
0192         evt = ascii_in->read_next_event();
0193       }
0194     }
0195
0196     if (!evt)
0197     {
0198       if (Verbosity() > 1)
0199       {
0200         std::cout << "Fun4AllHepMCInputManager::run::" << Name()
0201                   << ": error type: " << ascii_in->error_type()
0202                   << ", rdstate: " << ascii_in->rdstate() << std::endl;
0203       }
0204       fileclose();
0205     }
0206     else
0207     {
0208       MySyncManager()->CurrentEvent(evt->event_number());
0209       if (Verbosity() > 0)
0210       {
0211         std::cout << "Fun4AllHepMCInputManager::run::" << Name()
0212                   << ": hepmc evt no: " << evt->event_number() << std::endl;
0213       }
0214       m_MyEvent.push_back(evt->event_number());
0215       PHHepMCGenEventMap::Iter ievt = PHHepMCGenHelper::get_geneventmap()->find(PHHepMCGenHelper::get_embedding_id());
0216       if (ievt != PHHepMCGenHelper::get_geneventmap()->end())
0217       {
0218         // override existing event
0219         ievt->second->addEvent(evt);
0220       }
0221       else
0222       {
0223         PHHepMCGenHelper::insert_event(evt);
0224       }
0225
0226       events_total++;
0227       events_thisfile++;
0228
0229       // check if the local SubsysReco discards this event
0230       if (RejectEvent() != Fun4AllReturnCodes::EVENT_OK)
0231       {
0232         // if this event is discarded we only need to remove the event from the list of event numbers
0233         // the new event will overwrite the event on the node tree without issues
0234         m_MyEvent.pop_back();
0235       }
0236       else
0237       {
0238         break;  // have a good event, move on
0239       }
0240     }
0241   }  // attempt to retrieve a valid event from inputs
0242   return Fun4AllReturnCodes::EVENT_OK;
0243 }
0244
0245 int Fun4AllHepMCInputManager::fileclose()
0246 {
0247   if (!IsOpen())
0248   {
0249     std::cout << Name() << ": fileclose: No Input file open" << std::endl;
0250     return -1;
0251   }
0252   if (m_ReadOscarFlag)
0253   {
0254     theOscarFile.close();
0255   }
0256   else
0257   {
0258     delete ascii_in;
0259     ascii_in = nullptr;
0260   }
0261   IsOpen(0);
0262   // if we have a file list, move next entry to top of the list
0263   // or repeat the same entry again
0264   UpdateFileList();
0265   return 0;
0266 }
0267
0268 void Fun4AllHepMCInputManager::Print(const std::string &what) const
0269 {
0270   Fun4AllInputManager::Print(what);
0271   std::cout << Name() << " Vertex Settings: " << std::endl;
0272   PHHepMCGenHelper::Print(what);
0273   return;
0274 }
0275
0276 int Fun4AllHepMCInputManager::PushBackEvents(const int i)
0277 {
0278   // PushBackEvents is supposedly pushing events back on the stack which works
0279   // easily with root trees (just grab a different entry) but hard in these HepMC ASCII files.
0280   // A special case is when the synchronization fails and we need to only push back a single
0281   // event. In this case we save the evt in a temporary file which is read back in the run method
0282   // instead of getting the next event.
0283   if (i > 0)
0284   {
0285     if (i == 1 && evt)  // check on evt pointer makes sure it is not done from the cmd line
0286     {
0287       // root barfs when writing the node to the output.
0288       // Saving the pointer - even using a deep copy and reusing it did not work
0289       // The hackaround which works is to write this event into a temporary file and read it back
0290       if (m_HepMCTmpFile.empty())
0291       {
0292         // we need to create this filename just once, we reuse it. Do it only if we need it
0293         m_HepMCTmpFile = "/tmp/HepMCTmpEvent-" + Name() + "-" + std::to_string(getpid()) + ".hepmc";
0294       }
0295       HepMC::IO_GenEvent ascii_io(m_HepMCTmpFile, std::ios::out);
0296       ascii_io << evt;
0297       m_EventPushedBackFlag = 1;
0298       m_MyEvent.pop_back();
0299
0300       if (Verbosity() > 3)
0301       {
0302         std::cout << Name() << ": pushing back evt no " << evt->event_number() << std::endl;
0303       }
0304       return 0;
0305     }
0306     std::cout << PHWHERE << Name()
0307               << " Fun4AllHepMCInputManager cannot pop back more than 1 event"
0308               << std::endl;
0309     return -1;
0310   }
0311   if (!IsOpen())
0312   {
0313     std::cout << PHWHERE << Name()
0314               << " no file opened yet" << std::endl;
0315     return -1;
0316   }
0317   // Skipping events is implemented as
0318   // pushing a negative number of events on the stack, so in order to implement
0319   // the skipping of events we read -i events.
0320   int nevents = -i;  // negative number of events to push back -> skip num events
0321   int errorflag = 0;
0322   while (nevents > 0 && !errorflag)
0323   {
0324     evt = ascii_in->read_next_event();
0325     if (!evt)
0326     {
0327       std::cout << "Error after skipping " << i - nevents << std::endl;
0328       std::cout << "error type: " << ascii_in->error_type()
0329                 << ", rdstate: " << ascii_in->rdstate() << std::endl;
0330       errorflag = -1;
0331       fileclose();
0332     }
0333     else
0334     {
0335       m_MyEvent.push_back(evt->event_number());
0336       if (Verbosity() > 3)
0337       {
0338         std::cout << "Skipping evt no: " << evt->event_number() << std::endl;
0339       }
0340     }
0341     delete evt;
0342     nevents--;
0343   }
0344   return errorflag;
0345 }
0346
0347 HepMC::GenEvent *
0348 Fun4AllHepMCInputManager::ConvertFromOscar()
0349 {
0350   if (theOscarFile.eof())  // if the file is exhausted bail out during this next read
0351   {
0352     return nullptr;
0353   }
0354
0355   delete evt;
0356   //use PHENIX unit
0357   evt = new HepMC::GenEvent(HepMC::Units::GEV, HepMC::Units::CM);
0358
0359   if (Verbosity() > 1) std::cout << "Reading Oscar Event " << events_total << std::endl;
0360   //Grab New Event From Oscar
0361   std::string theLine;
0362   std::vector<std::vector<double> > theEventVec;
0363   std::vector<HepMC::FourVector> theVtxVec;
0364   while (getline(theOscarFile, theLine))
0365   {
0366     if (theLine.compare(0, 1, "#") == 0) continue;
0367     std::vector<double> theInfo;  //format: N,pid,px,py,pz,E,mass,xvtx,yvtx,zvtx,?
0368     double number = NAN;
0369     for (std::istringstream numbers_iss(theLine); numbers_iss >> number;)
0370     {
0371       theInfo.push_back(number);
0372     }
0373
0374     if (theInfo.size() == 2 && theInfo[0] == 0 && theInfo[1] == 0)
0375     {
0376       break;
0377     }
0378     else if (theInfo.size() == 2 && theInfo[0] == 0 && theInfo[1] > 0)
0379     {
0380       continue;
0381     }
0382     else
0383     {
0384       theEventVec.push_back(theInfo);
0385       HepMC::FourVector vert(theInfo[8] * toMM, theInfo[9] * toMM, theInfo[10] * toMM, theInfo[11]);
0386       theVtxVec.push_back(vert);
0387     }
0388
0389   }  //while(getline)
0390
0391   //Set Event Number
0392   evt->set_event_number(events_total);
0393
0394   //Loop Over One Event, Fill HepMC
0395   for (unsigned int i = 0; i < theEventVec.size(); i++)
0396   {
0397     //int N = (int)theEventVec[i][0];
0398     int pid = (int) theEventVec[i][1];
0399     double px = theEventVec[i][3];
0400     double py = theEventVec[i][4];
0401     double pz = theEventVec[i][5];
0402     double E = theEventVec[i][6];
0403     double m = theEventVec[i][7];
0404     int status = 1;  //oscar only writes final state particles
0405
0406     HepMC::GenVertex *v = new HepMC::GenVertex(theVtxVec[i]);
0407     evt->add_vertex(v);
0408
0409     HepMC::GenParticle *p = new HepMC::GenParticle(HepMC::FourVector(px, py, pz, E), pid, status);
0410     p->setGeneratedMass(m);
0411     p->suggest_barcode(i + 1);
0412     v->add_particle_out(p);
0413   }
0414   if (Verbosity() > 3)
0415   {
0416     evt->print();
0417   }
0418   return evt;
0419 }
0420
0421 int Fun4AllHepMCInputManager::ResetEvent()
0422 {
0423   m_MyEvent.clear();
0424   return 0;
0425 }
0426
0427 int Fun4AllHepMCInputManager::MyCurrentEvent(const unsigned int index) const
0428 {
0429   if (m_MyEvent.empty())
0430   {
0431     return 0;
0432   }
0433   return m_MyEvent.at(index);
0434 }
0435