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/AfterburnerAlgo.h>
0010 #include <flowafterburner/flowAfterburner.h>  // Afterburner class
0011
0012 #include <fun4all/Fun4AllReturnCodes.h>
0013 #include <fun4all/SubsysReco.h>  // for SubsysReco
0014
0015 #include <phool/PHRandomSeed.h>
0016 #include <phool/getClass.h>
0017 #include <phool/phool.h>
0018
0019 #include <CLHEP/Random/MTwistEngine.h>
0020 #include <CLHEP/Random/RandomEngine.h>
0021
0022 #include <iostream>
0023 #include <iterator>  // for operator!=, reverse_ite...
0024 #include <ranges>
0025 #include <set>  // for set, _Rb_tree_const_ite...
0026 #include <string>
0027 #include <utility>  // for pair
0028
0029 class PHCompositeNode;
0030 namespace HepMC
0031 {
0032   class GenEvent;
0033 }
0034
0035 // we want to keep the default eta range identical between here and
0036 // the flowAfterburner executable. If you change the default eta range here
0037 // please apply the same change to generators/flowAfterburner/main.cc
0038 HepMCFlowAfterBurner::HepMCFlowAfterBurner(const std::string &name)
0039   : SubsysReco(name)
0040 {
0041 }
0042
0043 int HepMCFlowAfterBurner::Init(PHCompositeNode * /*topNode*/)
0044 {
0045   if (seedset)
0046   {
0047     randomSeed = seed;
0048   }
0049   else
0050   {
0051     randomSeed = PHRandomSeed();
0052   }
0053
0054   m_engine = new CLHEP::MTwistEngine(randomSeed);
0055   m_afterburner = new Afterburner(algorithmName, m_engine, mineta, maxeta, minpt, maxpt);
0056   m_flowalgo = m_afterburner->getAlgo();
0057   // you can set other algo parameters here if needed
0058   if (enableFlucuations)
0059   {
0060     m_flowalgo->enable_fluctuations();
0061   }
0062
0063   for (unsigned int i = 1; i <= 6; ++i)
0064   {  // apply the scale factors to the flow harmonics if any
0065     if (flowScales[i - 1] != 1.0F)
0066     {
0067       m_flowalgo->set_single_scale_N(i, flowScales[i - 1]);
0068     }
0069   }
0070
0071   return 0;
0072 }
0073
0074 int HepMCFlowAfterBurner::process_event(PHCompositeNode *topNode)
0075 {
0076   PHHepMCGenEventMap *genevtmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0077   for (auto &iter : std::ranges::reverse_view(*genevtmap))
0078   {
0079     PHHepMCGenEvent *genevt = iter.second;
0080
0081     HepMC::GenEvent *evt = genevt->getEvent();
0082     if (!evt)
0083     {
0084       std::cout << PHWHERE << " no evt pointer under HEPMC Node found" << std::endl;
0085       return Fun4AllReturnCodes::ABORTEVENT;
0086     }
0087     if (Verbosity() > 0)
0088     {
0089       std::cout << "calling flowAfterburner with algorithm "
0090                 << algorithmName << ", mineta " << mineta
0091                 << ", maxeta: " << maxeta << ", minpt: " << minpt
0092                 << ", maxpt: " << maxpt << std::endl;
0093     }
0094
0095     m_afterburner->flowAfterburner(evt);
0096
0097     for (unsigned int i = 1; i <= 6; ++i)
0098     {
0099       genevt->set_flow_psi(i, m_afterburner->getPsiN(i));
0100       if (Verbosity() > 1)
0101       {
0102         std::cout << "  set reaction plane angle psi_" << i << " = " << genevt->get_flow_psi(i) << std::endl;
0103       }
0104     }
0105
0106     if (Verbosity() > 1)
0107     {
0108       m_afterburner->getAlgo()->print();
0109     }
0110
0111     // flowAfterburner(evt, engine, algorithmName, mineta, maxeta, minpt, maxpt);
0112   }
0113   return Fun4AllReturnCodes::EVENT_OK;
0114 }
0115
0116 void HepMCFlowAfterBurner::setSeed(const long il)
0117 {
0118   seedset = 1;
0119   seed = il;
0120   randomSeed = seed;
0121   // just in case we are already running, kill the engine and make
0122   // a new one using the selected seed
0123   if (m_engine)
0124   {
0125     delete m_engine;
0126     m_engine = new CLHEP::MTwistEngine(randomSeed);
0127   }
0128   return;
0129 }
0130
0131 void HepMCFlowAfterBurner::SaveRandomState(const std::string &savefile)
0132 {
0133   if (m_engine)
0134   {
0135     m_engine->saveStatus(savefile.c_str());
0136     return;
0137   }
0138   std::cout << PHWHERE << " Random engine not started yet" << std::endl;
0139 }
0140
0141 void HepMCFlowAfterBurner::RestoreRandomState(const std::string &savefile)
0142 {
0143   if (m_engine)
0144   {
0145     m_engine->restoreStatus(savefile.c_str());
0146     return;
0147   }
0148   std::cout << PHWHERE << " Random engine not started yet" << std::endl;
0149 }
0150
0151 void HepMCFlowAfterBurner::Print(const std::string & /*what*/) const
0152 {
0153   std::cout << "FlowAfterBurner parameters:" << std::endl;
0154   std::cout << "algorithm: " << algorithmName << std::endl;
0155   std::cout << "mineta: " << mineta << ", maxeta: " << maxeta << std::endl;
0156   std::cout << "minpt: " << minpt << ", maxpt: " << maxpt << std::endl;
0157   std::cout << "Implemented algorithms: MINBIAS (default), MINBIAS_V2_ONLY, CUSTOM"
0158             << std::endl;
0159   return;
0160 }
0161
0162 void HepMCFlowAfterBurner::setAlgorithmName(const std::string &name)
0163 {
0164   m_flowalgorithm = AfterburnerAlgo::getAlgoFromName(name);       // will print error message if algo name is unknown
0165   algorithmName = AfterburnerAlgo::getAlgoName(m_flowalgorithm);  // make sure the name is consistent
0166   return;
0167 }
0168
0169 void HepMCFlowAfterBurner::scaleFlow(const float scale, const unsigned int n)
0170 {
0171   if (n == 0)
0172   {  // set all scales
0173     for (unsigned int i = 0; i < 6; ++i)
0174     {
0175       flowScales[i] = scale;
0176     }
0177   }
0178   else if (n > 0 && n <= 6)
0179   {  // set specific harmonic
0180     flowScales[n - 1] = scale;
0181   }
0182   else
0183   {  // out of range
0184     std::cout << "HepMCFlowAfterBurner::scaleFlow - ERROR: n = " << n << " is out of range.  Must be between 0 (all) or 1,..,6." << std::endl;
0185   }
0186 }