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
0168       if (OpenNextFile())
0169       {
0170         std::cout << "Fun4AllHepMCInputManager::run::" << Name() << ": No Input file from filelist opened" << std::endl;
0171         return -1;
0172       }
0173     }
0174
0175     if (m_EventPushedBackFlag)  // if an event was pushed back, reuse save copy
0176     {
0177       HepMC::IO_GenEvent ascii_tmp_in(m_HepMCTmpFile, std::ios::in);
0178       ascii_tmp_in >> evt;
0179       m_EventPushedBackFlag = 0;
0180       remove(m_HepMCTmpFile.c_str());
0181     }
0182     else
0183     {
0184       if (m_ReadOscarFlag)
0185       {
0186         evt = ConvertFromOscar();
0187       }
0188       else
0189       {
0190         evt = ascii_in->read_next_event();
0191       }
0192     }
0193
0194     if (!evt)
0195     {
0196       if (Verbosity() > 1)
0197       {
0198         std::cout << "Fun4AllHepMCInputManager::run::" << Name()
0199                   << ": error type: " << ascii_in->error_type()
0200                   << ", rdstate: " << ascii_in->rdstate() << std::endl;
0201       }
0202       fileclose();
0203     }
0204     else
0205     {
0206       MySyncManager()->CurrentEvent(evt->event_number());
0207       if (Verbosity() > 0)
0208       {
0209         std::cout << "Fun4AllHepMCInputManager::run::" << Name()
0210                   << ": hepmc evt no: " << evt->event_number() << std::endl;
0211       }
0212       m_MyEvent.push_back(evt->event_number());
0213       PHHepMCGenEventMap::Iter ievt = PHHepMCGenHelper::get_geneventmap()->find(PHHepMCGenHelper::get_embedding_id());
0214       if (ievt != PHHepMCGenHelper::get_geneventmap()->end())
0215       {
0216         // override existing event
0217         ievt->second->addEvent(evt);
0218       }
0219       else
0220       {
0221         PHHepMCGenHelper::insert_event(evt);
0222       }
0223
0224       events_total++;
0225       events_thisfile++;
0226
0227       // check if the local SubsysReco discards this event
0228       if (RejectEvent() != Fun4AllReturnCodes::EVENT_OK)
0229       {
0230         // if this event is discarded we only need to remove the event from the list of event numbers
0231         // the new event will overwrite the event on the node tree without issues
0232         m_MyEvent.pop_back();
0233       }
0234       else
0235       {
0236         break;  // have a good event, move on
0237       }
0238     }
0239   }  // attempt to retrieve a valid event from inputs
0240   return Fun4AllReturnCodes::EVENT_OK;
0241 }
0242
0243 int Fun4AllHepMCInputManager::fileclose()
0244 {
0245   if (!IsOpen())
0246   {
0247     std::cout << Name() << ": fileclose: No Input file open" << std::endl;
0248     return -1;
0249   }
0250   if (m_ReadOscarFlag)
0251   {
0252     theOscarFile.close();
0253   }
0254   else
0255   {
0256     delete ascii_in;
0257     ascii_in = nullptr;
0258   }
0259   IsOpen(0);
0260   // if we have a file list, move next entry to top of the list
0261   // or repeat the same entry again
0262   UpdateFileList();
0263   return 0;
0264 }
0265
0266 void Fun4AllHepMCInputManager::Print(const std::string &what) const
0267 {
0268   Fun4AllInputManager::Print(what);
0269   std::cout << Name() << " Vertex Settings: " << std::endl;
0270   PHHepMCGenHelper::Print(what);
0271   return;
0272 }
0273
0274 int Fun4AllHepMCInputManager::PushBackEvents(const int i)
0275 {
0276   // PushBackEvents is supposedly pushing events back on the stack which works
0277   // easily with root trees (just grab a different entry) but hard in these HepMC ASCII files.
0278   // A special case is when the synchronization fails and we need to only push back a single
0279   // event. In this case we save the evt in a temporary file which is read back in the run method
0280   // instead of getting the next event.
0281   if (i > 0)
0282   {
0283     if (i == 1 && evt)  // check on evt pointer makes sure it is not done from the cmd line
0284     {
0285       // root barfs when writing the node to the output.
0286       // Saving the pointer - even using a deep copy and reusing it did not work
0287       // The hackaround which works is to write this event into a temporary file and read it back
0288       if (m_HepMCTmpFile.empty())
0289       {
0290         // we need to create this filename just once, we reuse it. Do it only if we need it
0291         m_HepMCTmpFile = "/tmp/HepMCTmpEvent-" + Name() + "-" + std::to_string(getpid()) + ".hepmc";
0292       }
0293       HepMC::IO_GenEvent ascii_io(m_HepMCTmpFile, std::ios::out);
0294       ascii_io << evt;
0295       m_EventPushedBackFlag = 1;
0296       m_MyEvent.pop_back();
0297
0298       if (Verbosity() > 3)
0299       {
0300         std::cout << Name() << ": pushing back evt no " << evt->event_number() << std::endl;
0301       }
0302       return 0;
0303     }
0304     std::cout << PHWHERE << Name()
0305               << " Fun4AllHepMCInputManager cannot pop back more than 1 event"
0306               << std::endl;
0307     return -1;
0308   }
0309   if (!IsOpen())
0310   {
0311     std::cout << PHWHERE << Name()
0312               << " no file opened yet" << std::endl;
0313     return -1;
0314   }
0315   // Skipping events is implemented as
0316   // pushing a negative number of events on the stack, so in order to implement
0317   // the skipping of events we read -i events.
0318   int nevents = -i;  // negative number of events to push back -> skip num events
0319   int errorflag = 0;
0320   while (nevents > 0 && !errorflag)
0321   {
0322     evt = ascii_in->read_next_event();
0323     if (!evt)
0324     {
0325       std::cout << "Error after skipping " << i - nevents << std::endl;
0326       std::cout << "error type: " << ascii_in->error_type()
0327                 << ", rdstate: " << ascii_in->rdstate() << std::endl;
0328       errorflag = -1;
0329       fileclose();
0330     }
0331     else
0332     {
0333       m_MyEvent.push_back(evt->event_number());
0334       if (Verbosity() > 3)
0335       {
0336         std::cout << "Skipping evt no: " << evt->event_number() << std::endl;
0337       }
0338     }
0339     delete evt;
0340     nevents--;
0341   }
0342   return errorflag;
0343 }
0344
0345 HepMC::GenEvent *
0346 Fun4AllHepMCInputManager::ConvertFromOscar()
0347 {
0348   if (theOscarFile.eof())  // if the file is exhausted bail out during this next read
0349   {
0350     return nullptr;
0351   }
0352
0353   delete evt;
0354   // use PHENIX unit
0355   evt = new HepMC::GenEvent(HepMC::Units::GEV, HepMC::Units::CM);
0356
0357   if (Verbosity() > 1)
0358   {
0359     std::cout << "Reading Oscar Event " << events_total << std::endl;
0360   }
0361   // Grab New Event From Oscar
0362   std::string theLine;
0363   std::vector<std::vector<double> > theEventVec;
0364   std::vector<HepMC::FourVector> theVtxVec;
0365   while (getline(theOscarFile, theLine))
0366   {
0367     if (theLine.compare(0, 1, "#") == 0)
0368     {
0369       continue;
0370     }
0371     std::vector<double> theInfo;  // format: N,pid,px,py,pz,E,mass,xvtx,yvtx,zvtx,?
0372     double number = NAN;
0373     for (std::istringstream numbers_iss(theLine); numbers_iss >> number;)
0374     {
0375       theInfo.push_back(number);
0376     }
0377
0378     if (theInfo.size() == 2 && theInfo[0] == 0 && theInfo[1] == 0)
0379     {
0380       break;
0381     }
0382     if (theInfo.size() == 2 && theInfo[0] == 0 && theInfo[1] > 0)
0383     {
0384       continue;
0385     }
0386
0387     theEventVec.push_back(theInfo);
0388     HepMC::FourVector vert(theInfo[8] * toMM, theInfo[9] * toMM, theInfo[10] * toMM, theInfo[11]);
0389     theVtxVec.push_back(vert);
0390
0391   }  // while(getline)
0392
0393   // Set Event Number
0394   evt->set_event_number(events_total);
0395
0396   // Loop Over One Event, Fill HepMC
0397   for (unsigned int i = 0; i < theEventVec.size(); i++)
0398   {
0399     // int N = (int)theEventVec[i][0];
0400     int pid = (int) theEventVec[i][1];
0401     double px = theEventVec[i][3];
0402     double py = theEventVec[i][4];
0403     double pz = theEventVec[i][5];
0404     double E = theEventVec[i][6];
0405     double m = theEventVec[i][7];
0406     int status = 1;  // oscar only writes final state particles
0407
0408     HepMC::GenVertex *v = new HepMC::GenVertex(theVtxVec[i]);
0409     evt->add_vertex(v);
0410
0411     HepMC::GenParticle *p = new HepMC::GenParticle(HepMC::FourVector(px, py, pz, E), pid, status);
0412     p->setGeneratedMass(m);
0413     p->suggest_barcode(i + 1);
0414     v->add_particle_out(p);
0415   }
0416   if (Verbosity() > 3)
0417   {
0418     evt->print();
0419   }
0420   return evt;
0421 }
0422
0423 int Fun4AllHepMCInputManager::ResetEvent()
0424 {
0425   m_MyEvent.clear();
0426   return 0;
0427 }
0428
0429 int Fun4AllHepMCInputManager::MyCurrentEvent(const unsigned int index) const
0430 {
0431   if (m_MyEvent.empty())
0432   {
0433     return 0;
0434   }
0435   return m_MyEvent.at(index);
0436 }