g4simulation/g4main/Fun4AllDstPileupInputManager.cc

0001 /*!
0002  * \file Fun4AllDstPileupInputManager.cc
0003  * \author Hugo Pereira Da Costa <hugo.pereira-da-costa@cea.fr>
0004  */
0005
0006 #include "Fun4AllDstPileupInputManager.h"
0007 #include "Fun4AllDstPileupMerger.h"
0008
0009 #include <fun4all/Fun4AllInputManager.h>  // for Fun4AllInputManager
0010 #include <fun4all/Fun4AllReturnCodes.h>
0011 #include <fun4all/Fun4AllServer.h>
0012
0013 #include <ffaobjects/RunHeader.h>
0014
0015 #include <frog/FROG.h>
0016
0017 #include <phool/PHCompositeNode.h>
0018 #include <phool/PHNodeIOManager.h>
0019 #include <phool/PHNodeIntegrate.h>
0020 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0021 #include <phool/PHRandomSeed.h>
0022 #include <phool/getClass.h>
0023 #include <phool/phool.h>  // for PHWHERE, PHReadOnly, PHRunTree
0024
0025 #include <gsl/gsl_randist.h>
0026
0027 #include <cassert>
0028 #include <iostream>  // for operator<<, basic_ostream, endl
0029 #include <utility>   // for pair
0030
0031 //_____________________________________________________________________________
0032 Fun4AllDstPileupInputManager::Fun4AllDstPileupInputManager(const std::string &name, const std::string &nodename, const std::string &topnodename)
0033   : Fun4AllInputManager(name, nodename, topnodename)
0034 {
0035   // initialize random generator
0036   const uint seed = PHRandomSeed();
0037   m_rng.reset(gsl_rng_alloc(gsl_rng_mt19937));
0038   gsl_rng_set(m_rng.get(), seed);
0039 }
0040
0041 //_____________________________________________________________________________
0042 int Fun4AllDstPileupInputManager::fileopen(const std::string &filenam)
0043 {
0044   /*
0045   this is largely copied from fun4all/Fun4AllDstInputManager::fileopen
0046   with additional code to handle the background IManager
0047   */
0048
0049   auto se = Fun4AllServer::instance();
0050   if (IsOpen())
0051   {
0052     std::cout << "Closing currently open file "
0053               << FileName()
0054               << " and opening " << filenam << std::endl;
0055     fileclose();
0056   }
0057   FileName(filenam);
0058   FROG frog;
0059   m_fullfilename = frog.location(FileName());
0060   if (Verbosity() > 0)
0061   {
0062     std::cout << Name() << ": opening file " << m_fullfilename << std::endl;
0063   }
0064   // sanity check - the IManager must be nullptr when this method is executed
0065   // if not something is very very wrong and we must not continue
0066   assert(!m_IManager);
0067
0068   // first read the runnode if not disabled
0069   if (m_ReadRunTTree)
0070   {
0071     m_IManager.reset(new PHNodeIOManager(m_fullfilename, PHReadOnly, PHRunTree));
0072     if (m_IManager->isFunctional())
0073     {
0074       m_runNode = se->getNode(m_RunNode, TopNodeName());
0075       m_IManager->read(m_runNode);
0076
0077       // get the current run number
0078       auto runheader = findNode::getClass<RunHeader>(m_runNode, "RunHeader");
0079       if (runheader)
0080       {
0081         SetRunNumber(runheader->get_RunNumber());
0082       }
0083       // delete our internal copy of the runnode when opening subsequent files
0084       assert(!m_runNodeCopy);
0085       m_runNodeCopy.reset(new PHCompositeNode("RUNNODECOPY"));
0086       if (!m_runNodeSum)
0087       {
0088         m_runNodeSum.reset(new PHCompositeNode("RUNNODESUM"));
0089       }
0090
0091       {
0092         // read run node using temporary node iomanager
0093         PHNodeIOManager(m_fullfilename, PHReadOnly, PHRunTree).read(m_runNodeCopy.get());
0094       }
0095
0096       PHNodeIntegrate integrate;
0097       integrate.RunNode(m_runNode);
0098       integrate.RunSumNode(m_runNodeSum.get());
0099       // run recursively over internal run node copy and integrate objects
0100       PHNodeIterator mainIter(m_runNodeCopy.get());
0101       mainIter.forEach(integrate);
0102       // we do not need to keep the internal copy, keeping it would crate
0103       // problems in case a subsequent file does not contain all the
0104       // runwise objects from the previous file. Keeping this copy would then
0105       // integrate the missing object again with the old copy
0106       m_runNodeCopy.reset();
0107     }
0108     // DLW: move the delete outside the if block to cover the case where isFunctional() fails
0109     m_IManager.reset();
0110   }
0111
0112   // create internal dst node
0113   if (!m_dstNodeInternal)
0114   {
0115     m_dstNodeInternal.reset(new PHCompositeNode("DST_INTERNAL"));
0116   }
0117
0118   // update dst node from fun4all server
0119   m_dstNode = se->getNode(InputNode(), TopNodeName());
0120
0121   // open file in both active and background input manager
0122   m_IManager.reset(new PHNodeIOManager(m_fullfilename, PHReadOnly));
0123   if (m_IManager->isFunctional())
0124   {
0125     IsOpen(1);
0126     m_ievent_thisfile = 0;
0127     setBranches();                // set branch selections
0128     AddToFileOpened(FileName());  // add file to the list of files which were opened
0129     return 0;
0130   }
0131   else
0132   {
0133     std::cout << PHWHERE << ": " << Name() << " Could not open file " << FileName() << std::endl;
0134     m_IManager.reset();
0135     return -1;
0136   }
0137 }
0138
0139 //_____________________________________________________________________________
0140 int Fun4AllDstPileupInputManager::run(const int nevents)
0141 {
0142   if (nevents == 0)
0143   {
0144     return runOne(nevents);
0145   }
0146   else if (nevents > 1)
0147   {
0148     const auto result = runOne(nevents - 1);
0149     if (result != 0)
0150     {
0151       return result;
0152     }
0153   }
0154
0155   /*
0156    * assign/create relevant dst nodes if not already there
0157    * this normally happens in ::fileopen however, when the file is not oppened during first event, for instance because background rate is too low,
0158    * this can cause fun4all server to bark with "Someone changed the number of Output Nodes on the fly"
0159    */
0160   if (!m_dstNode)
0161   {
0162     auto se = Fun4AllServer::instance();
0163     m_dstNode = se->getNode(InputNode(), TopNodeName());
0164   }
0165
0166   if (!m_dstNodeInternal)
0167   {
0168     m_dstNodeInternal.reset(new PHCompositeNode("DST_INTERNAL"));
0169   }
0170
0171   // create merger node
0172   Fun4AllDstPileupMerger merger;
0173   merger.copyDetectorActiveCrossings(m_DetectorTiming);
0174   merger.load_nodes(m_dstNode);
0175
0176   // generate background collisions
0177   const double mu = m_collision_rate * m_time_between_crossings * 1e-9;
0178
0179   const int min_crossing = m_tmin / m_time_between_crossings;
0180   const int max_crossing = m_tmax / m_time_between_crossings;
0181   for (int icrossing = min_crossing; icrossing <= max_crossing; ++icrossing)
0182   {
0183     const double crossing_time = m_time_between_crossings * icrossing;
0184     const int ncollisions = gsl_ran_poisson(m_rng.get(), mu);
0185     for (int icollision = 0; icollision < ncollisions; ++icollision)
0186     {
0187       // read one event
0188       const auto result = runOne(1);
0189       if (result != 0)
0190       {
0191         return result;
0192       }
0193
0194       // merge
0195       if (Verbosity() > 0)
0196       {
0197         std::cout << "Fun4AllDstPileupInputManager::run - merged background event " << m_ievent_thisfile << " time: " << crossing_time << std::endl;
0198       }
0199       merger.copy_background_event(m_dstNodeInternal.get(), crossing_time);
0200     }
0201   }
0202
0203   return 0;
0204 }
0205
0206 //_____________________________________________________________________________
0207 int Fun4AllDstPileupInputManager::fileclose()
0208 {
0209   if (!IsOpen())
0210   {
0211     std::cout << Name() << ": fileclose: No Input file open" << std::endl;
0212     return -1;
0213   }
0214   m_IManager.reset();
0215   IsOpen(0);
0216   UpdateFileList();
0217   return 0;
0218 }
0219
0220 //_____________________________________________________________________________
0221 int Fun4AllDstPileupInputManager::BranchSelect(const std::string &branch, const int iflag)
0222 {
0223   int myflag = iflag;
0224   // if iflag > 0 the branch is set to read
0225   // if iflag = 0, the branch is set to NOT read
0226   // if iflag < 0 the branchname is erased from our internal branch read map
0227   // this does not have any effect on phool yet
0228   if (myflag < 0)
0229   {
0230     std::map<const std::string, int>::iterator branchiter;
0231     branchiter = m_branchread.find(branch);
0232     if (branchiter != m_branchread.end())
0233     {
0234       m_branchread.erase(branchiter);
0235     }
0236     return 0;
0237   }
0238
0239   if (myflag > 0)
0240   {
0241     if (Verbosity() > 1)
0242     {
0243       std::cout << "Setting Root Tree Branch: " << branch << " to read" << std::endl;
0244     }
0245     myflag = 1;
0246   }
0247   else
0248   {
0249     if (Verbosity() > 1)
0250     {
0251       std::cout << "Setting Root Tree Branch: " << branch << " to NOT read" << std::endl;
0252     }
0253   }
0254   m_branchread[branch] = myflag;
0255   return 0;
0256 }
0257
0258 //_____________________________________________________________________________
0259 int Fun4AllDstPileupInputManager::setBranches()
0260 {
0261   if (m_IManager)
0262   {
0263     if (!m_branchread.empty())
0264     {
0265       std::map<const std::string, int>::const_iterator branchiter;
0266       for (branchiter = m_branchread.begin(); branchiter != m_branchread.end(); ++branchiter)
0267       {
0268         m_IManager->selectObjectToRead(branchiter->first.c_str(), branchiter->second);
0269         if (Verbosity() > 0)
0270         {
0271           std::cout << branchiter->first << " set to " << branchiter->second << std::endl;
0272         }
0273       }
0274     }
0275   }
0276   else
0277   {
0278     std::cout << PHWHERE << " " << Name() << ": You can only call this function after a file has been opened" << std::endl;
0279     std::cout << "Do not worry, the branches will be set as soon as you open a file" << std::endl;
0280     return -1;
0281   }
0282   return 0;
0283 }
0284
0285 //_____________________________________________________________________________
0286 void Fun4AllDstPileupInputManager::Print(const std::string &what) const
0287 {
0288   if (what == "ALL" || what == "BRANCH")
0289   {
0290     // loop over the map and print out the content (name and location in memory)
0291     std::cout << "--------------------------------------" << std::endl
0292               << std::endl;
0293     std::cout << "List of selected branches in Fun4AllDstPileupInputManager " << Name() << ":" << std::endl;
0294
0295     std::map<const std::string, int>::const_iterator iter;
0296     for (iter = m_branchread.begin(); iter != m_branchread.end(); ++iter)
0297     {
0298       std::cout << iter->first << " is switched ";
0299       if (iter->second)
0300       {
0301         std::cout << "ON";
0302       }
0303       else
0304       {
0305         std::cout << "OFF";
0306       }
0307       std::cout << std::endl;
0308     }
0309   }
0310   if ((what == "ALL" || what == "PHOOL") && m_IManager)
0311   {
0312     // loop over the map and print out the content (name and location in memory)
0313     std::cout << "--------------------------------------" << std::endl
0314               << std::endl;
0315     std::cout << "PHNodeIOManager print in Fun4AllDstPileupInputManager " << Name() << ":" << std::endl;
0316     m_IManager->print();
0317   }
0318   Fun4AllInputManager::Print(what);
0319   return;
0320 }
0321
0322 //_____________________________________________________________________________
0323 int Fun4AllDstPileupInputManager::PushBackEvents(const int i)
0324 {
0325   if (m_IManager)
0326   {
0327     unsigned EventOnDst = m_IManager->getEventNumber();
0328     EventOnDst -= static_cast<unsigned>(i);
0329     m_ievent_thisfile -= i;
0330     m_ievent_total -= i;
0331     m_IManager->setEventNumber(EventOnDst);
0332     return 0;
0333   }
0334   std::cout << PHWHERE << Name() << ": could not push back events, Imanager is NULL"
0335             << " probably the dst is not open yet (you need to call fileopen or run 1 event for lists)" << std::endl;
0336   return -1;
0337 }
0338
0339 //_____________________________________________________________________________
0340 int Fun4AllDstPileupInputManager::runOne(const int nevents)
0341 {
0342   if (!IsOpen())
0343   {
0344     if (FileListEmpty())
0345     {
0346       if (Verbosity() > 0)
0347       {
0348         std::cout << Name() << ": No Input file open" << std::endl;
0349       }
0350       return -1;
0351     }
0352     else
0353     {
0354       if (OpenNextFile())
0355       {
0356         std::cout << Name() << ": No Input file from filelist opened" << std::endl;
0357         return -1;
0358       }
0359     }
0360   }
0361   if (Verbosity() > 3)
0362   {
0363     std::cout << "Getting Event from " << Name() << std::endl;
0364   }
0365
0366 readagain:
0367
0368   // read main event to dstNode
0369   auto dummy = m_IManager->read(m_dstNodeInternal.get());
0370   int ncount = 0;
0371   while (dummy)
0372   {
0373     ++ncount;
0374     if (nevents > 0 && ncount >= nevents)
0375     {
0376       break;
0377     }
0378     dummy = m_IManager->read(m_dstNodeInternal.get());
0379   }
0380   if (!dummy)
0381   {
0382     fileclose();
0383     if (!OpenNextFile())
0384     {
0385       // NOLINTNEXTLINE(hicpp-avoid-goto)
0386       goto readagain;
0387     }
0388     return -1;
0389   }
0390   m_ievent_total += ncount;
0391   m_ievent_thisfile += ncount;
0392   // check if the local SubsysReco discards this event
0393   if (RejectEvent() != Fun4AllReturnCodes::EVENT_OK)
0394   {
0395     // NOLINTNEXTLINE(hicpp-avoid-goto)
0396     goto readagain;
0397   }
0398   return 0;
0399 }
0400
0401 void Fun4AllDstPileupInputManager::setDetectorActiveCrossings(const std::string &name, const int nbcross)
0402 {
0403   setDetectorActiveCrossings(name, -nbcross, nbcross);
0404 }
0405
0406 void Fun4AllDstPileupInputManager::setDetectorActiveCrossings(const std::string &name, const int min, const int max)
0407 {
0408   std::string nodename = "G4HIT_" + name;
0409   // compensate that active for one bunch crossign means delta_t = 0
0410   m_DetectorTiming.insert(std::make_pair(nodename, std::make_pair(m_time_between_crossings * (min + 1), m_time_between_crossings * (max - 1))));
0411   return;
0412 }