generators/phhepmc/HepMCFlowAfterBurner.cc

0001 //
0002 // Inspired by code from ATLAS.  Thanks!
0003 //
0004 #include "HepMCFlowAfterBurner.h"
0005
0006 #include "PHHepMCGenEvent.h"
0007 #include "PHHepMCGenEventMap.h"
0008
0009 #include <flowafterburner/flowAfterburner.h>
0010
0011 #include <fun4all/Fun4AllReturnCodes.h>
0012 #include <fun4all/SubsysReco.h>               // for SubsysReco
0013
0014 #include <phool/PHRandomSeed.h>
0015 #include <phool/getClass.h>
0016 #include <phool/phool.h>
0017
0018 #include <CLHEP/Random/MTwistEngine.h>
0019 #include <CLHEP/Random/RandomEngine.h>
0020
0021 #include <iostream>
0022 #include <iterator>                           // for operator!=, reverse_ite...
0023 #include <set>                                // for set, _Rb_tree_const_ite...
0024 #include <string>
0025 #include <utility>                            // for pair
0026
0027 using namespace std;
0028
0029 class PHCompositeNode;
0030 namespace HepMC { class GenEvent; }
0031
0032 CLHEP::HepRandomEngine *engine = nullptr;
0033
0034 set<string> algoset = {"MINBIAS", "MINBIAS_V2_ONLY", "CUSTOM"};
0035
0036 // we want to keep the default eta range identical between here and
0037 // the flowAfterburner executable. If you change the default eta range here
0038 // please apply the same change to generators/flowAfterburner/main.cc
0039 HepMCFlowAfterBurner::HepMCFlowAfterBurner(const std::string &name)
0040   : SubsysReco(name)
0041 {
0042 }
0043
0044 int HepMCFlowAfterBurner::Init(PHCompositeNode */*topNode*/)
0045 {
0046   if (seedset)
0047   {
0048     randomSeed = seed;
0049   }
0050   else
0051   {
0052     randomSeed = PHRandomSeed();
0053   }
0054
0055   engine = new CLHEP::MTwistEngine(randomSeed);
0056
0057   return 0;
0058 }
0059
0060 int HepMCFlowAfterBurner::process_event(PHCompositeNode *topNode)
0061 {
0062   PHHepMCGenEventMap *genevtmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0063   for (PHHepMCGenEventMap::ReverseIter iter = genevtmap->rbegin(); iter != genevtmap->rend(); ++iter)
0064   {
0065     PHHepMCGenEvent *genevt = iter->second;
0066
0067     HepMC::GenEvent *evt = genevt->getEvent();
0068     if (!evt)
0069     {
0070       cout << PHWHERE << " no evt pointer under HEPMC Node found" << endl;
0071       return Fun4AllReturnCodes::ABORTEVENT;
0072     }
0073     if (Verbosity() > 0)
0074     {
0075       cout << "calling flowAfterburner with algorithm "
0076            << algorithmName << ", mineta " << mineta
0077            << ", maxeta: " << maxeta << ", minpt: " << minpt
0078            << ", maxpt: " << maxpt << endl;
0079     }
0080     flowAfterburner(evt, engine, algorithmName, mineta, maxeta, minpt, maxpt);
0081   }
0082   return Fun4AllReturnCodes::EVENT_OK;
0083 }
0084
0085 void HepMCFlowAfterBurner::setSeed(const long il)
0086 {
0087   seedset = 1;
0088   seed = il;
0089   randomSeed = seed;
0090   // just in case we are already running, kill the engine and make
0091   // a new one using the selected seed
0092   if (engine)
0093   {
0094     delete engine;
0095     engine = new CLHEP::MTwistEngine(randomSeed);
0096   }
0097   return;
0098 }
0099
0100 void HepMCFlowAfterBurner::SaveRandomState(const string &savefile)
0101 {
0102   if (engine)
0103   {
0104     engine->saveStatus(savefile.c_str());
0105     return;
0106   }
0107   cout << PHWHERE << " Random engine not started yet" << endl;
0108 }
0109
0110 void HepMCFlowAfterBurner::RestoreRandomState(const string &savefile)
0111 {
0112   if (engine)
0113   {
0114     engine->restoreStatus(savefile.c_str());
0115     return;
0116   }
0117   cout << PHWHERE << " Random engine not started yet" << endl;
0118 }
0119
0120 void HepMCFlowAfterBurner::Print(const string &/*what*/) const
0121 {
0122   cout << "FlowAfterBurner parameters:" << endl;
0123   cout << "algorithm: " << algorithmName << endl;
0124   cout << "mineta: " << mineta << ", maxeta: " << maxeta << endl;
0125   cout << "minpt: " << minpt << ", maxpt: " << maxpt << endl;
0126   cout << "Implemented algorithms: MINBIAS (default), MINBIAS_V2_ONLY, CUSTOM"
0127        << endl;
0128   return;
0129 }
0130
0131 void HepMCFlowAfterBurner::setAlgorithmName(const std::string &name)
0132 {
0133   auto it = algoset.find(name);
0134   if (it != algoset.end())
0135   {
0136     algorithmName = *it;
0137   }
0138   else
0139   {
0140     cout << "algorithm " << name << " not in list of possible algorithms" << endl;
0141     cout << "possible algorithms are" << endl;
0142     for (auto &al : algoset)
0143     {
0144       cout << al << endl;
0145     }
0146   }
0147   return;
0148 }