generators/PHHerwig7/PHHerwig7.cc

0001 #include "PHHerwig7.h"
0002 #include "PHHerwig7Core.h"
0003 #include "PHHerwig7GenTrigger.h"
0004
0005 #include <phhepmc/PHGenIntegral.h>
0006 #include <phhepmc/PHGenIntegralv1.h>
0007 #include <phhepmc/PHHepMCGenHelper.h>
0008
0009 #include <fun4all/Fun4AllReturnCodes.h>
0010
0011 #include <phool/PHCompositeNode.h>
0012 #include <phool/PHIODataNode.h>
0013 #include <phool/PHNodeIterator.h>
0014 #include <phool/PHObject.h>
0015 #include <phool/PHRandomSeed.h>
0016 #include <phool/getClass.h>
0017 #include <phool/phool.h>
0018
0019 // HepMC2 flavour of ThePEG's converter.  HepMCDefs.h must come first so the
0020 // traits see HEPMC_HAS_CROSS_SECTION and carry sigma into the record.
0021 #include <HepMC/GenEvent.h>
0022 #include <HepMC/HepMCDefs.h>
0023 #include <HepMC/Units.h>
0024 #include <HepMC/WeightContainer.h>
0025 #include <ThePEG/Config/HepMCHelper.h>
0026 #include <ThePEG/Vectors/HepMCConverter.h>
0027
0028 #include <cassert>
0029 #include <cstdlib>
0030 #include <iostream>
0031
0032 PHHerwig7::PHHerwig7(const std::string &name)
0033   : SubsysReco(name)
0034 {
0035   PHHepMCGenHelper::set_embedding_id(1);  // same default as PHPythia8
0036 }
0037
0038 int PHHerwig7::Init(PHCompositeNode *topNode)
0039 {
0040   if (m_RunFile.empty())
0041   {
0042     std::cout << PHWHERE << " no Herwig .run file set (PHHerwig7::set_run_file)" << std::endl;
0043     exit(EXIT_FAILURE);
0044   }
0045
0046   if (create_node_tree(topNode) != Fun4AllReturnCodes::EVENT_OK)
0047   {
0048     exit(EXIT_FAILURE);
0049   }
0050
0051   int seed = m_Seed;
0052   if (seed <= 0)
0053   {
0054     unsigned int s = PHRandomSeed();
0055     // ThePEG takes an int seed through the UI; keep it positive and in range
0056     seed = static_cast<int>(s % 2000000000U);
0057     if (seed <= 0)
0058     {
0059       seed = 1;
0060     }
0061   }
0062   std::cout << Name() << " Herwig random seed: " << seed << std::endl;
0063
0064   m_Core = std::make_unique<PHHerwig7Core>();
0065   m_Core->ui().m_RunFile = m_RunFile;
0066   m_Core->ui().m_SetupFile = m_SetupFile;
0067   m_Core->ui().m_Prepend = m_ReadDirs;
0068   m_Core->ui().m_Seed = seed;
0069
0070   // Herwig's cout formatting is not restored by ThePEG
0071   std::ios old_state(nullptr);
0072   old_state.copyfmt(std::cout);
0073   try
0074   {
0075     m_Core->init();
0076   }
0077   catch (const std::exception &e)
0078   {
0079     std::cout.copyfmt(old_state);
0080     std::cout << PHWHERE << " Herwig initialisation failed: " << e.what() << std::endl;
0081     // Fun4All drops a module whose Init fails and carries on writing empty
0082     // events; a generator that cannot start must kill the job (as PHPythia8
0083     // does on a bad seed).
0084     std::cout << PHWHERE << " Check the .run file and LHAPDF_DATA_PATH (hwrun.sh sets it to Source/lhapdf)" << std::endl;
0085     exit(EXIT_FAILURE);
0086   }
0087   std::cout.copyfmt(old_state);
0088
0089   if (Verbosity() > 0)
0090   {
0091     Print("CONFIG");
0092   }
0093   return Fun4AllReturnCodes::EVENT_OK;
0094 }
0095
0096 int PHHerwig7::process_event(PHCompositeNode * /*topNode*/)
0097 {
0098   std::ios old_state(nullptr);
0099   old_state.copyfmt(std::cout);
0100
0101   HepMC::GenEvent *genevent = nullptr;
0102   long attempts = 0;
0103   while (!genevent)
0104   {
0105     ThePEG::EventPtr ev;
0106     try
0107     {
0108       ev = m_Core->next();
0109     }
0110     catch (const std::exception &e)
0111     {
0112       std::cout.copyfmt(old_state);
0113       std::cout << PHWHERE << " Herwig threw during event generation: " << e.what() << std::endl;
0114       return Fun4AllReturnCodes::ABORTRUN;
0115     }
0116     if (!ev)
0117     {
0118       // shoot() returns null only when the generator has hit its own event
0119       // limit or given up; retrying would loop forever.
0120       std::cout.copyfmt(old_state);
0121       std::cout << PHWHERE << " Herwig returned no event after " << m_NGenerated
0122                 << " generated (generator event limit reached?), aborting run" << std::endl;
0123       return Fun4AllReturnCodes::ABORTRUN;
0124     }
0125     ++m_NGenerated;
0126     const double w = ev->weight();
0127     m_SumWAll += w;
0128
0129     HepMC::GenEvent *cand = ThePEG::HepMCConverter<HepMC::GenEvent>::convert(
0130         *ev, false, ThePEG::Units::GeV, ThePEG::Units::millimeter);
0131     if (!cand)
0132     {
0133       std::cout << PHWHERE << " HepMC conversion failed" << std::endl;
0134       continue;
0135     }
0136     if (m_SaveEventWeightFlag && cand->weights().empty())
0137     {
0138       cand->weights().push_back(w);
0139     }
0140
0141     // Apply the registered triggers to the event
0142     bool passedTrigger = false;
0143     bool andScoreKeeper = true;
0144
0145     for (auto *trigger : m_RegisteredTriggers)
0146     {
0147       const bool trigResult = trigger->Apply(cand);
0148
0149       if (Verbosity() > 1)
0150       {
0151         std::cout << Name() << ": trigger " << trigger->GetName() << " pass = " << trigResult << std::endl;
0152       }
0153
0154       if (m_TriggersOR && trigResult)
0155       {
0156         passedTrigger = true;
0157         break;
0158       }
0159
0160       if (m_TriggersAND)
0161       {
0162         andScoreKeeper = andScoreKeeper && trigResult;
0163       }
0164     }
0165
0166     if ((andScoreKeeper && m_TriggersAND) || m_RegisteredTriggers.empty())
0167     {
0168       passedTrigger = true;
0169     }
0170
0171     if (passedTrigger)
0172     {
0173       genevent = cand;
0174       m_SumWPass += w;
0175     }
0176     else
0177     {
0178       delete cand;
0179       ++attempts;
0180       if (m_MaxAttempts > 0 && attempts >= m_MaxAttempts)
0181       {
0182         std::cout.copyfmt(old_state);
0183         std::cout << PHWHERE << " failed to pass trigger after " << attempts
0184                   << " attempts, aborting run" << std::endl;
0185         return Fun4AllReturnCodes::ABORTRUN;
0186       }
0187     }
0188   }
0189
0190   PHHepMCGenEvent *success = PHHepMCGenHelper::insert_event(genevent);
0191   if (!success)
0192   {
0193     std::cout.copyfmt(old_state);
0194     std::cout << PHWHERE << " failed to add event to the HepMC record" << std::endl;
0195     return Fun4AllReturnCodes::ABORTRUN;
0196   }
0197   ++m_NPassed;
0198
0199   if (m_IntegralNode)
0200   {
0201     const double sigma_pb = m_Core->xsec_pb();
0202     m_IntegralNode->set_N_Generator_Accepted_Event(m_NGenerated);
0203     m_IntegralNode->set_N_Processed_Event(m_NPassed);
0204     m_IntegralNode->set_Sum_Of_Weight(m_SumWAll);
0205     // PHPythia8 convention: pb^-1
0206     m_IntegralNode->set_Integrated_Lumi(sigma_pb > 0 ? double(m_NGenerated) / sigma_pb : 0.0);
0207   }
0208
0209   std::cout.copyfmt(old_state);
0210   return Fun4AllReturnCodes::EVENT_OK;
0211 }
0212
0213 int PHHerwig7::End(PHCompositeNode * /*topNode*/)
0214 {
0215   std::ios old_state(nullptr);
0216   old_state.copyfmt(std::cout);
0217   if (m_Core)
0218   {
0219     m_Core->finish();
0220   }
0221   std::cout.copyfmt(old_state);
0222
0223   std::cout << Name() << ": generated " << m_NGenerated << ", passed trigger " << m_NPassed
0224             << " (eff " << (m_NGenerated ? double(m_NPassed) / double(m_NGenerated) : 0.0) << ")"
0225             << ", sumW all=" << m_SumWAll << " passed=" << m_SumWPass
0226             << ", sigma=" << (m_Core ? m_Core->xsec_pb() : 0.0) << " +- "
0227             << (m_Core ? m_Core->xsec_err_pb() : 0.0) << " pb" << std::endl;
0228   if (m_IntegralNode && Verbosity() > 0)
0229   {
0230     m_IntegralNode->identify();
0231   }
0232   return Fun4AllReturnCodes::EVENT_OK;
0233 }
0234
0235 int PHHerwig7::create_node_tree(PHCompositeNode *topNode)
0236 {
0237   PHHepMCGenHelper::create_node_tree(topNode);
0238
0239   PHNodeIterator iter(topNode);
0240   PHCompositeNode *sumNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0241   if (!sumNode)
0242   {
0243     std::cout << PHWHERE << " RUN node missing, doing nothing" << std::endl;
0244     return Fun4AllReturnCodes::ABORTRUN;
0245   }
0246   if (m_SaveIntegratedLuminosityFlag)
0247   {
0248     m_IntegralNode = findNode::getClass<PHGenIntegral>(sumNode, "PHGenIntegral");
0249     if (!m_IntegralNode)
0250     {
0251       m_IntegralNode = new PHGenIntegralv1("PHHerwig7 with embedding ID of " + std::to_string(PHHepMCGenHelper::get_embedding_id()));
0252       PHIODataNode<PHObject> *newnode = new PHIODataNode<PHObject>(m_IntegralNode, "PHGenIntegral", "PHObject");
0253       sumNode->addNode(newnode);
0254     }
0255     else
0256     {
0257       std::cout << PHWHERE << " RUN/PHGenIntegral already exists; refusing to overwrite. "
0258                 << "Call PHHerwig7::save_integrated_luminosity(false) if that is intended." << std::endl;
0259       return Fun4AllReturnCodes::ABORTRUN;
0260     }
0261   }
0262   return Fun4AllReturnCodes::EVENT_OK;
0263 }
0264
0265 void PHHerwig7::Print(const std::string & /*what*/) const
0266 {
0267   std::cout << Name() << " configuration:" << std::endl
0268             << "  run file       : " << m_RunFile << std::endl
0269             << "  setup file     : " << (m_SetupFile.empty() ? "(none)" : m_SetupFile) << std::endl
0270             << "  embedding id   : " << PHHepMCGenHelper::get_embedding_id() << std::endl
0271             << "  triggers       : " << m_RegisteredTriggers.size() << std::endl
0272             << "  combine        : " << (m_TriggersAND ? "AND" : "OR") << std::endl;
0273
0274   for (auto *trigger : m_RegisteredTriggers)
0275   {
0276     std::cout << "  trigger        : " << trigger->GetName() << std::endl;
0277   }
0278 }
0279
0280 // Same register_trigger() as PHPythia8
0281 void PHHerwig7::register_trigger(PHHerwig7GenTrigger *trigger)
0282 {
0283   if (Verbosity() > 0)
0284   {
0285     std::cout << Name() << ": trigger " << trigger->GetName() << " registered" << std::endl;
0286   }
0287
0288   m_RegisteredTriggers.push_back(trigger);
0289 }