generators/phhepmc/Fun4AllHepMCPileupInputManager.cc

0001 #include "Fun4AllHepMCPileupInputManager.h"
0002
0003 #include "PHHepMCGenEvent.h"
0004 #include "PHHepMCGenEventMap.h"
0005 #include "PHHepMCGenHelper.h"  // for PHHepMCGenHelper, PHHepMCGen...
0006
0007 #include <fun4all/Fun4AllBase.h>  // for Fun4AllBase::VERBOSITY_SOME
0008 #include <fun4all/Fun4AllReturnCodes.h>
0009
0010 #include <phool/PHRandomSeed.h>
0011
0012 #include <HepMC/GenEvent.h>
0013 #include <HepMC/IO_GenEvent.h>
0014
0015 #include <gsl/gsl_randist.h>
0016 #include <gsl/gsl_rng.h>
0017
0018 #include <cassert>  // for assert
0019 #include <fstream>
0020 #include <iostream>
0021
0022 Fun4AllHepMCPileupInputManager::Fun4AllHepMCPileupInputManager(
0023     const std::string &name, const std::string &nodename, const std::string &topnodename)
0024   : Fun4AllHepMCInputManager(name, nodename, topnodename)
0025   , RandomGenerator(gsl_rng_alloc(gsl_rng_mt19937))
0026 {
0027   //! repeatedly read the input file
0028   Repeat(1);
0029
0030   //! If set_embedding_id(i) with a negative number or 0, the pile up event will be inserted with increasing positive embedding_id. This is the default operation mode.
0031   PHHepMCGenHelper::set_embedding_id(-1);
0032
0033   unsigned int seed = PHRandomSeed();  // fixed seed is handled in this funtcion
0034   gsl_rng_set(RandomGenerator, seed);
0035
0036   return;
0037 }
0038
0039 Fun4AllHepMCPileupInputManager::~Fun4AllHepMCPileupInputManager()
0040 {
0041   gsl_rng_free(RandomGenerator);
0042 }
0043
0044 int Fun4AllHepMCPileupInputManager::SkipForThisManager(const int nevents)
0045 {
0046   for (int i = 0; i < nevents; ++i)
0047   {
0048     if (m_SignalInputManager)
0049     {
0050       m_SignalEventNumber = m_SignalInputManager->MyCurrentEvent(i);
0051     }
0052     int iret = run(1, true);
0053     if (iret)
0054     {
0055       return iret;
0056     }
0057   }
0058   return 0;
0059 }
0060
0061 int Fun4AllHepMCPileupInputManager::run(const int /*nevents*/, const bool skip)
0062 {
0063   if (_first_run)
0064   {
0065     _first_run = false;
0066
0067     _ave_coll_per_crossing = _collision_rate * _time_between_crossings * 1e-9;
0068     _min_crossing = _min_integration_time / _time_between_crossings;
0069     _max_crossing = _max_integration_time / _time_between_crossings;
0070
0071     if (Verbosity() >= VERBOSITY_SOME)
0072     {
0073       std::cout << "Fun4AllHepMCPileupInputManager::run:first event:  - ";
0074       std::cout << " _ave_coll_per_crossing = " << _ave_coll_per_crossing;
0075       std::cout << " _min_crossing = " << _min_crossing;
0076       std::cout << " _max_crossing = " << _max_crossing;
0077       std::cout << ". Start first event." << std::endl;
0078     }
0079   }
0080   if (m_SignalInputManager && !skip)
0081   {
0082     m_SignalEventNumber = m_SignalInputManager->MyCurrentEvent();
0083   }
0084   // if we have pushed back events we read them here and return
0085   if (m_EventPushedBackFlag)
0086   {
0087     // if m_EventPushedBackFlag = -1 we have no events but still created the temporary HepMC file
0088     // so we still need to remove it
0089     // if m_EventPushedBackFlag > 0 we have events pushed back which need to be recovered
0090     if (m_EventPushedBackFlag > 0)
0091     {
0092       HepMC::IO_GenEvent ascii_tmp_in(m_HepMCTmpFile, std::ios::in);
0093       HepMC::GenEvent *evttmp = ascii_tmp_in.read_next_event();
0094       while (ascii_tmp_in.rdstate() == 0)
0095       {
0096         if (evttmp->event_number() != m_SignalEventNumber)
0097         {
0098           double crossing_time = m_EventNumberMap.find(evttmp->event_number())->second;
0099           InsertEvent(evttmp, crossing_time);
0100         }
0101         evttmp = ascii_tmp_in.read_next_event();
0102       }
0103     }
0104     m_EventNumberMap.clear();
0105     m_EventPushedBackFlag = 0;
0106     remove(m_HepMCTmpFile.c_str());
0107     return 0;
0108   }
0109   // toss multiple crossings all the way back
0110   for (int icrossing = _min_crossing; icrossing <= _max_crossing; ++icrossing)
0111   {
0112     double crossing_time = _time_between_crossings * icrossing;
0113
0114     int ncollisions = gsl_ran_poisson(RandomGenerator, _ave_coll_per_crossing);
0115     //    if (icrossing == 0) --ncollisions; // Jin: this is not necessary.
0116     // Triggering an interesting crossing does not change the distribution of the number of background collisions in that crossing
0117
0118     for (int icollision = 0; icollision < ncollisions; ++icollision)
0119     {
0120       // loop until retrieve a valid event
0121       while (true)
0122       {
0123         if (!IsOpen())
0124         {
0125           if (FileListEmpty())
0126           {
0127             if (Verbosity() > 0)
0128             {
0129               std::cout << Name() << ": No Input file open" << std::endl;
0130             }
0131             return -1;
0132           }
0133
0134           if (OpenNextFile())
0135           {
0136             std::cout << Name() << ": No Input file from filelist opened" << std::endl;
0137             return -1;
0138           }
0139         }
0140         {
0141           if (ReadOscar())
0142           {
0143             evt = ConvertFromOscar();
0144             if (evt && m_SignalEventNumber == evt->event_number())
0145             {
0146               delete evt;
0147               evt = ConvertFromOscar();
0148             }
0149           }
0150           else
0151           {
0152             evt = ascii_in->read_next_event();
0153             if (evt && m_SignalEventNumber == evt->event_number())
0154             {
0155               delete evt;
0156               evt = ascii_in->read_next_event();
0157             }
0158           }
0159         }
0160
0161         if (!evt)
0162         {
0163           if (Verbosity() > 1)
0164           {
0165             std::cout << "error type: " << ascii_in->error_type()
0166                       << ", rdstate: " << ascii_in->rdstate() << std::endl;
0167           }
0168           fileclose();
0169         }
0170         else
0171         {
0172           if (Verbosity() > 0)
0173           {
0174             std::cout << "Fun4AllHepMCPileupInputManager::run::" << Name();
0175             if (skip)
0176             {
0177               std::cout << " skip";
0178             }
0179             std::cout << " hepmc evt no: " << evt->event_number() << std::endl;
0180           }
0181           events_total++;
0182           events_thisfile++;
0183           // check if the local SubsysReco discards this event
0184           if (RejectEvent() != Fun4AllReturnCodes::EVENT_OK)
0185           {
0186             delete evt;
0187             evt = nullptr;
0188             ResetEvent();
0189           }
0190           else
0191           {
0192             m_EventNumberMap.insert(std::make_pair(evt->event_number(), crossing_time));
0193             break;  // got the evt, move on
0194           }
0195         }
0196       }  // loop until retrieve a valid event
0197       if (!skip)
0198       {
0199         InsertEvent(evt, crossing_time);
0200       }
0201       else
0202       {
0203         delete evt;
0204         evt = nullptr;
0205       }
0206     }  //    for (int icollision = 0; icollision < ncollisions; ++icollision)
0207
0208   }  //  for (int icrossing = _min_crossing; icrossing <= _max_crossing; ++icrossing)
0209   return 0;
0210 }
0211
0212 int Fun4AllHepMCPileupInputManager::ResetEvent()
0213 {
0214   m_EventNumberMap.clear();
0215   m_SignalEventNumber = 0;
0216   return 0;
0217 }
0218
0219 int Fun4AllHepMCPileupInputManager::PushBackEvents(const int i)
0220 {
0221   if (i == 1)
0222   {
0223     if (m_HepMCTmpFile.empty())
0224     {
0225       // we need to create this filename just once, we reuse it. Do it only if we need it
0226       m_HepMCTmpFile = "/tmp/HepMCTmpPileEvent-" + Name() + "-" + std::to_string(getpid()) + ".hepmc";
0227     }
0228     m_EventPushedBackFlag = -1;
0229     HepMC::IO_GenEvent ascii_io(m_HepMCTmpFile, std::ios::out);
0230     PHHepMCGenEventMap *geneventmap = PHHepMCGenHelper::get_geneventmap();
0231     for (auto &iter : *geneventmap)
0232     {
0233       if (m_EventNumberMap.contains((iter.second)->getEvent()->event_number()))
0234       {
0235         m_EventPushedBackFlag = 1;
0236         HepMC::GenEvent *evttmp = (iter.second)->getEvent();
0237         ascii_io << evttmp;
0238       }
0239     }
0240     return 0;
0241   }
0242   return -1;
0243 }
0244 int Fun4AllHepMCPileupInputManager::InsertEvent(HepMC::GenEvent *evt, const double crossing_time)
0245 {
0246   PHHepMCGenEventMap *geneventmap = PHHepMCGenHelper::get_geneventmap();
0247   PHHepMCGenEvent *genevent = nullptr;
0248   if (PHHepMCGenHelper::get_embedding_id() > 0)
0249   {
0250     //! If set_embedding_id(i) with a positive number, the pile up event will be inserted with increasing positive embedding_id. This would be a strange way to use pile up.
0251
0252     genevent = geneventmap->insert_active_event(get_PHHepMCGenEvent_template());
0253   }
0254   else
0255   {
0256     //! If set_embedding_id(i) with a negative number or 0, the pile up event will be inserted with increasing positive embedding_id. This is the default operation mode.
0257
0258     genevent = geneventmap->insert_background_event(get_PHHepMCGenEvent_template());
0259   }
0260   assert(genevent);
0261   assert(evt);
0262   genevent->addEvent(evt);
0263   PHHepMCGenHelper::HepMC2Lab_boost_rotation_translation(genevent);
0264   // place to the crossing center in time
0265   genevent->moveVertex(0, 0, 0, crossing_time);
0266   return 0;
0267 }