generators/PHPythia8/PHPythia8.cc

0001 #include "PHPythia8.h"
0002
0003 #include "PHPy8GenTrigger.h"
0004
0005 #include <phhepmc/PHGenIntegral.h>  // for PHGenIntegral
0006 #include <phhepmc/PHGenIntegralv1.h>
0007 #include <phhepmc/PHHepMCGenHelper.h>  // for PHHepMCGenHelper
0008
0009 #include <fun4all/Fun4AllReturnCodes.h>
0010 #include <fun4all/SubsysReco.h>  // for SubsysReco
0011
0012 #include <phool/PHCompositeNode.h>
0013 #include <phool/PHIODataNode.h>
0014 #include <phool/PHNode.h>          // for PHNode
0015 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0016 #include <phool/PHObject.h>        // for PHObject
0017 #include <phool/PHRandomSeed.h>
0018 #include <phool/getClass.h>
0019 #include <phool/phool.h>  // for PHWHERE
0020
0021 #include <HepMC/GenEvent.h>
0022 #include <HepMC/Units.h>            // for GEV, MM
0023 #include <HepMC/WeightContainer.h>  // for WeightContainer
0024
0025 #include <Pythia8/Event.h>  // for Event
0026 #include <Pythia8/Info.h>   // for Info
0027 #include <Pythia8/Pythia.h>
0028 #include <Pythia8Plugins/HepMC2.h>
0029
0030 #include <cassert>
0031 #include <cstdlib>
0032 #include <format>
0033 #include <fstream>
0034 #include <iostream>  // for operator<<, endl
0035
0036 PHPythia8::PHPythia8(const std::string &name)
0037   : SubsysReco(name)
0038 {
0039   char *charPath = getenv("PYTHIA8");
0040   if (!charPath)
0041   {
0042     std::cout << "PHPythia8::Could not find $PYTHIA8 path!" << std::endl;
0043     return;
0044   }
0045
0046   std::string thePath(charPath);
0047   thePath += "/xmldoc/";
0048   m_Pythia8.reset(new Pythia8::Pythia(thePath));
0049
0050   m_Pythia8ToHepMC.reset(new HepMC::Pythia8ToHepMC());
0051   m_Pythia8ToHepMC->set_store_proc(true);
0052   m_Pythia8ToHepMC->set_store_pdf(true);
0053   m_Pythia8ToHepMC->set_store_xsec(true);
0054
0055   PHHepMCGenHelper::set_embedding_id(1);  // default embedding ID to 1
0056 }
0057
0058 int PHPythia8::Init(PHCompositeNode *topNode)
0059 {
0060   if (!m_ConfigFileName.empty())
0061   {
0062     read_config(m_ConfigFileName);
0063   }
0064   for (auto &m_Command : m_Commands)
0065   {
0066     m_Pythia8->readString(m_Command);
0067   }
0068
0069   create_node_tree(topNode);
0070
0071   // PYTHIA8 has very specific requires for its random number range
0072   // I map the designated unique seed from recoconst into something
0073   // acceptable for PYTHIA8
0074
0075   unsigned int seed = PHRandomSeed();
0076
0077   if (seed > 900000000)
0078   {
0079     seed = seed % 900000000;
0080   }
0081
0082   if ((seed > 0) && (seed <= 900000000))
0083   {
0084     m_Pythia8->readString("Random:setSeed = on");
0085     m_Pythia8->readString(std::format("Random:seed = {}", seed));
0086   }
0087   else
0088   {
0089     std::cout << PHWHERE << " ERROR: seed " << seed << " is not valid" << std::endl;
0090     exit(1);
0091   }
0092   // print out seed so we can make this is reproducible
0093   std::cout << "PHPythia8 random seed: " << seed << std::endl;
0094
0095   m_Pythia8->init();
0096
0097   return Fun4AllReturnCodes::EVENT_OK;
0098 }
0099
0100 int PHPythia8::End(PHCompositeNode * /*topNode*/)
0101 {
0102   if (Verbosity() >= VERBOSITY_MORE)
0103   {
0104     std::cout << "PHPythia8::End - I'm here!" << std::endl;
0105   }
0106
0107   if (Verbosity() >= VERBOSITY_SOME)
0108   {
0109     //-* dump out closing info (cross-sections, etc)
0110     m_Pythia8->stat();
0111
0112     // match pythia printout
0113     std::cout << " |                                                                "
0114               << "                                                 | " << std::endl;
0115     std::cout << "                         PHPythia8::End - " << m_EventCount
0116               << " events passed trigger" << std::endl;
0117     std::cout << "                         Fraction passed: " << m_EventCount
0118               << "/" << m_Pythia8->info.nAccepted()
0119               << " = " << m_EventCount / float(m_Pythia8->info.nAccepted()) << std::endl;
0120     std::cout << " *-------  End PYTHIA Trigger Statistics  ------------------------"
0121               << "-------------------------------------------------* " << std::endl;
0122
0123     if (m_IntegralNode)
0124     {
0125       std::cout << "Integral information on stored on node RUN/PHGenIntegral:" << std::endl;
0126       m_IntegralNode->identify();
0127       std::cout << " *-------  End PYTHIA Integral Node Print  ------------------------"
0128                 << "-------------------------------------------------* " << std::endl;
0129     }
0130   }
0131
0132   return Fun4AllReturnCodes::EVENT_OK;
0133 }
0134
0135 //__________________________________________________________
0136 int PHPythia8::read_config(const std::string &cfg_file)
0137 {
0138   m_ConfigFileName = cfg_file;
0139
0140   if (Verbosity() >= VERBOSITY_SOME)
0141   {
0142     std::cout << "PHPythia8::read_config - Reading " << m_ConfigFileName << std::endl;
0143   }
0144
0145   std::ifstream infile(m_ConfigFileName);
0146   if (infile.fail())
0147   {
0148     std::cout << "PHPythia8::read_config - Failed to open file " << m_ConfigFileName << std::endl;
0149     exit(2);
0150   }
0151
0152   m_Pythia8->readFile(m_ConfigFileName);
0153
0154   return Fun4AllReturnCodes::EVENT_OK;
0155 }
0156
0157 //-* print pythia config info
0158 void PHPythia8::print_config() const
0159 {
0160   m_Pythia8->info.list();
0161 }
0162
0163 int PHPythia8::process_event(PHCompositeNode * /*topNode*/)
0164 {
0165   if (Verbosity() >= VERBOSITY_MORE)
0166   {
0167     std::cout << "PHPythia8::process_event - event: " << m_EventCount << std::endl;
0168   }
0169
0170   bool passedGen = false;
0171   bool passedTrigger = false;
0172   //  int genCounter = 0;
0173
0174   while (!passedTrigger)
0175   {
0176     //    ++genCounter;
0177
0178     // generate another pythia event
0179     while (!passedGen)
0180     {
0181       passedGen = m_Pythia8->next();
0182     }
0183
0184     // test trigger logic
0185     bool andScoreKeeper = true;
0186     if (Verbosity() >= VERBOSITY_EVEN_MORE)
0187     {
0188       std::cout << "PHPythia8::process_event - triggersize: " << m_RegisteredTriggers.size() << std::endl;
0189     }
0190
0191     for (auto &m_RegisteredTrigger : m_RegisteredTriggers)
0192     {
0193       bool trigResult = m_RegisteredTrigger->Apply(m_Pythia8.get());
0194
0195       if (Verbosity() >= VERBOSITY_EVEN_MORE)
0196       {
0197         std::cout << "PHPythia8::process_event trigger: "
0198                   << m_RegisteredTrigger->GetName() << "  " << trigResult << std::endl;
0199       }
0200
0201       if (m_TriggersOR && trigResult)
0202       {
0203         passedTrigger = true;
0204         break;
0205       }
0206       if (m_TriggersAND)
0207       {
0208         andScoreKeeper &= trigResult;
0209       }
0210
0211       if (Verbosity() >= VERBOSITY_EVEN_MORE && !passedTrigger)
0212       {
0213         std::cout << "PHPythia8::process_event - failed trigger: "
0214                   << m_RegisteredTrigger->GetName() << std::endl;
0215       }
0216     }
0217
0218     if ((andScoreKeeper && m_TriggersAND) || (m_RegisteredTriggers.empty()))
0219     {
0220       passedTrigger = true;
0221       //      genCounter = 0;
0222     }
0223
0224     passedGen = false;
0225   }
0226
0227   // print
0228   if (Verbosity())
0229   {
0230     m_Pythia8->event.list();
0231   }
0232
0233   // fill HepMC object with event & pass to
0234
0235   auto *genevent = new HepMC::GenEvent(HepMC::Units::GEV, HepMC::Units::MM);
0236   m_Pythia8ToHepMC->fill_next_event(*m_Pythia8, genevent, m_EventCount);
0237   // Enable continuous reweighting by storing additional reweighting factor
0238   if (m_SaveEventWeightFlag)
0239   {
0240     genevent->weights().push_back(m_Pythia8->info.weight());
0241   }
0242
0243   /* pass HepMC to PHNode*/
0244   auto *success = PHHepMCGenHelper::insert_event(genevent);
0245   if (!success)
0246   {
0247     std::cout << "PHPythia8::process_event - Failed to add event to HepMC record!" << std::endl;
0248     return Fun4AllReturnCodes::ABORTRUN;
0249   }
0250
0251   // print outs
0252
0253   if (Verbosity() >= VERBOSITY_MORE)
0254   {
0255     std::cout << "PHPythia8::process_event - FINISHED WHOLE EVENT" << std::endl;
0256   }
0257   if (m_EventCount < 2 && Verbosity() >= VERBOSITY_SOME)
0258   {
0259     m_Pythia8->event.list();
0260   }
0261   if (m_EventCount >= 2 && Verbosity() >= VERBOSITY_A_LOT)
0262   {
0263     m_Pythia8->event.list();
0264   }
0265
0266   ++m_EventCount;
0267
0268   // save statistics
0269   if (m_IntegralNode)
0270   {
0271     m_IntegralNode->set_N_Generator_Accepted_Event(m_Pythia8->info.nAccepted());
0272     m_IntegralNode->set_N_Processed_Event(m_EventCount);
0273     m_IntegralNode->set_Sum_Of_Weight(m_Pythia8->info.weightSum());
0274     m_IntegralNode->set_Integrated_Lumi(m_Pythia8->info.nAccepted() / (m_Pythia8->info.sigmaGen() * 1e9));
0275   }
0276
0277   return Fun4AllReturnCodes::EVENT_OK;
0278 }
0279
0280 int PHPythia8::create_node_tree(PHCompositeNode *topNode)
0281 {
0282   // HepMC IO
0283   PHHepMCGenHelper::create_node_tree(topNode);
0284
0285   PHNodeIterator iter(topNode);
0286   PHCompositeNode *sumNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0287   if (!sumNode)
0288   {
0289     std::cout << PHWHERE << "RUN Node missing doing nothing" << std::endl;
0290     return Fun4AllReturnCodes::ABORTRUN;
0291   }
0292   if (m_SaveIntegratedLuminosityFlag)
0293   {
0294     m_IntegralNode = findNode::getClass<PHGenIntegral>(sumNode, "PHGenIntegral");
0295     if (!m_IntegralNode)
0296     {
0297       m_IntegralNode = new PHGenIntegralv1("PHPythia8 with embedding ID of " + std::to_string(PHHepMCGenHelper::get_embedding_id()));
0298       PHIODataNode<PHObject> *newmapnode = new PHIODataNode<PHObject>(m_IntegralNode, "PHGenIntegral", "PHObject");
0299       sumNode->addNode(newmapnode);
0300     }
0301     else
0302     {
0303       std::cout << "PHPythia8::create_node_tree - Fatal Error - "
0304                 << "RUN/PHGenIntegral node already exist. "
0305                 << "It is messy to overwrite integrated luminosities. Please turn off this function in the macro with " << std::endl;
0306       std::cout << "                              PHPythia8::save_integrated_luminosity(false);" << std::endl;
0307       std::cout << "The current RUN/PHGenIntegral node is ";
0308       m_IntegralNode->identify(std::cout);
0309
0310       exit(EXIT_FAILURE);
0311     }
0312     assert(m_IntegralNode);
0313   }
0314   return Fun4AllReturnCodes::EVENT_OK;
0315 }
0316
0317 void PHPythia8::register_trigger(PHPy8GenTrigger *theTrigger)
0318 {
0319   if (Verbosity() >= VERBOSITY_MORE)
0320   {
0321     std::cout << "PHPythia8::registerTrigger - trigger " << theTrigger->GetName() << " registered" << std::endl;
0322   }
0323   m_RegisteredTriggers.push_back(theTrigger);
0324 }