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 {
0026   //! repeatedly read the input file
0027   Repeat(1);
0028
0029   //! 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.
0030   PHHepMCGenHelper::set_embedding_id(-1);
0031
0032   RandomGenerator = gsl_rng_alloc(gsl_rng_mt19937);
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           else
0134           {
0135             if (OpenNextFile())
0136             {
0137               std::cout << Name() << ": No Input file from filelist opened" << std::endl;
0138               return -1;
0139             }
0140           }
0141         }
0142         {
0143           if (ReadOscar())
0144           {
0145             evt = ConvertFromOscar();
0146             if (evt && m_SignalEventNumber == evt->event_number())
0147             {
0148               delete evt;
0149               evt = ConvertFromOscar();
0150             }
0151           }
0152           else
0153           {
0154             evt = ascii_in->read_next_event();
0155             if (evt && m_SignalEventNumber == evt->event_number())
0156             {
0157               delete evt;
0158               evt = ascii_in->read_next_event();
0159             }
0160           }
0161         }
0162
0163         if (!evt)
0164         {
0165           if (Verbosity() > 1)
0166           {
0167             std::cout << "error type: " << ascii_in->error_type()
0168                  << ", rdstate: " << ascii_in->rdstate() << std::endl;
0169           }
0170           fileclose();
0171         }
0172         else
0173         {
0174           if (Verbosity() > 0)
0175           {
0176             std::cout << "Fun4AllHepMCPileupInputManager::run::" << Name();
0177             if (skip) std::cout << " skip";
0178             std::cout << " hepmc evt no: " << evt->event_number() << std::endl;
0179           }
0180           events_total++;
0181           events_thisfile++;
0182           // check if the local SubsysReco discards this event
0183           if (RejectEvent() != Fun4AllReturnCodes::EVENT_OK)
0184           {
0185             delete evt;
0186             evt = nullptr;
0187             ResetEvent();
0188           }
0189           else
0190           {
0191             m_EventNumberMap.insert(std::make_pair(evt->event_number(), crossing_time));
0192             break;  // got the evt, move on
0193           }
0194         }
0195       }  // loop until retrieve a valid event
0196       if (!skip)
0197       {
0198         InsertEvent(evt, crossing_time);
0199       }
0200       else
0201       {
0202     delete evt;
0203     evt = nullptr;
0204       }
0205     }  //    for (int icollision = 0; icollision < ncollisions; ++icollision)
0206
0207   }  //  for (int icrossing = _min_crossing; icrossing <= _max_crossing; ++icrossing)
0208   return 0;
0209 }
0210
0211 int Fun4AllHepMCPileupInputManager::ResetEvent()
0212 {
0213   m_EventNumberMap.clear();
0214   m_SignalEventNumber = 0;
0215   return 0;
0216 }
0217
0218 int Fun4AllHepMCPileupInputManager::PushBackEvents(const int i)
0219 {
0220   if (i == 1)
0221   {
0222     if (m_HepMCTmpFile.empty())
0223     {
0224       // we need to create this filename just once, we reuse it. Do it only if we need it
0225       m_HepMCTmpFile = "/tmp/HepMCTmpPileEvent-" + Name() + "-" + std::to_string(getpid()) + ".hepmc";
0226     }
0227     m_EventPushedBackFlag = -1;
0228     HepMC::IO_GenEvent ascii_io(m_HepMCTmpFile, std::ios::out);
0229     PHHepMCGenEventMap *geneventmap = PHHepMCGenHelper::get_geneventmap();
0230     for (auto iter = geneventmap->begin(); iter != geneventmap->end(); ++iter)
0231     {
0232       if (m_EventNumberMap.find((iter->second)->getEvent()->event_number()) != m_EventNumberMap.end())
0233       {
0234         m_EventPushedBackFlag = 1;
0235         HepMC::GenEvent *evttmp = (iter->second)->getEvent();
0236         ascii_io << evttmp;
0237       }
0238     }
0239     return 0;
0240   }
0241   return -1;
0242 }
0243 int Fun4AllHepMCPileupInputManager::InsertEvent(HepMC::GenEvent *evt, const double crossing_time)
0244 {
0245   PHHepMCGenEventMap *geneventmap = PHHepMCGenHelper::get_geneventmap();
0246   PHHepMCGenEvent *genevent = nullptr;
0247   if (PHHepMCGenHelper::get_embedding_id() > 0)
0248   {
0249     //! 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.
0250
0251     genevent = geneventmap->insert_active_event(get_PHHepMCGenEvent_template() );
0252   }
0253   else
0254   {
0255     //! 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.
0256
0257     genevent = geneventmap->insert_background_event(get_PHHepMCGenEvent_template() );
0258   }
0259   assert(genevent);
0260   assert(evt);
0261   genevent->addEvent(evt);
0262   PHHepMCGenHelper::HepMC2Lab_boost_rotation_translation(genevent);
0263   // place to the crossing center in time
0264   genevent->moveVertex(0, 0, 0, crossing_time);
0265   return 0;
0266 }