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