Herwig7Interface/src/Herwig7HepMC3Interface.cc

0001 /** \class Herwig7Interface
0002  *
0003  *  Marco A. Harrendorf marco.harrendorf@cern.ch
0004  *  Dominik Beutel dominik.beutel@cern.ch
0005  */
0006
0007 #include <cmath>
0008 #include <cstdlib>
0009 #include <fstream>
0010 #include <iostream>
0011 #include <memory>
0012 #include <sstream>
0013
0014 #include <algorithm>
0015
0016 #include <HepMC3/GenEvent.h>
0017
0018 #include <Herwig/API/HerwigAPI.h>
0019
0020 #include <ThePEG/Utilities/DynamicLoader.h>
0021 #include <ThePEG/Repository/Repository.h>
0022 #include <ThePEG/Handlers/EventHandler.h>
0023 #include <ThePEG/Handlers/XComb.h>
0024 #include <ThePEG/EventRecord/Event.h>
0025 #include <ThePEG/EventRecord/Particle.h>
0026 #include <ThePEG/EventRecord/Collision.h>
0027 #include <ThePEG/EventRecord/TmpTransform.h>
0028 #include <ThePEG/Config/ThePEG.h>
0029 #include <ThePEG/PDF/PartonExtractor.h>
0030 #include <ThePEG/PDF/PDFBase.h>
0031 #include <ThePEG/Utilities/UtilityBase.h>
0032
0033 #include <GeneratorInterface/Herwig7Interface/interface/HepMC3Converter.h>
0034
0035 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0036
0037 #include "GeneratorInterface/Core/interface/ParameterCollector.h"
0038
0039 #include "GeneratorInterface/Herwig7Interface/interface/Proxy.h"
0040 #include "GeneratorInterface/Herwig7Interface/interface/RandomEngineGlue.h"
0041 #include "GeneratorInterface/Herwig7Interface/interface/Herwig7HepMC3Interface.h"
0042
0043 #include "CLHEP/Random/RandomEngine.h"
0044
0045 using namespace std;
0046 using namespace gen;
0047
0048 Herwig7HepMC3Interface::Herwig7HepMC3Interface(const edm::ParameterSet &pset)
0049     : randomEngineGlueProxy_(ThePEG::RandomEngineGlue::Proxy::create()),
0050       dataLocation_(ParameterCollector::resolve(pset.getParameter<string>("dataLocation"))),
0051       generator_(pset.getParameter<string>("generatorModule")),
0052       run_(pset.getParameter<string>("run")),
0053       dumpConfig_(pset.getUntrackedParameter<string>("dumpConfig", "HerwigConfig.in")),
0054       skipEvents_(pset.getUntrackedParameter<unsigned int>("skipEvents", 0)) {
0055   // Write events in hepmc ascii format for debugging purposes
0056 #if 0
0057   string dumpEvents = pset.getUntrackedParameter<string>("dumpEvents", "");
0058   if (!dumpEvents.empty()) {
0059     iobc_ = std::make_unique<HepMC::IO_GenEvent>(dumpEvents, ios::out);
0060     edm::LogInfo("ThePEGSource") << "Event logging switched on (=> " << dumpEvents << ")";
0061   }
0062 #endif
0063   // Clear dumpConfig target
0064   if (!dumpConfig_.empty())
0065     ofstream cfgDump(dumpConfig_.c_str(), ios_base::trunc);
0066 }
0067
0068 Herwig7HepMC3Interface::~Herwig7HepMC3Interface() noexcept {
0069   if (eg_)
0070     eg_->finalize();
0071   edm::LogInfo("Herwig7HepMC3Interface") << "Event generator finalized";
0072 }
0073
0074 void Herwig7HepMC3Interface::setPEGRandomEngine(CLHEP::HepRandomEngine *v) {
0075   randomEngineGlueProxy_->setRandomEngine(v);
0076   randomEngine = v;
0077   ThePEG::RandomEngineGlue *rnd = randomEngineGlueProxy_->getInstance();
0078   if (rnd) {
0079     rnd->setRandomEngine(v);
0080   }
0081 }
0082
0083 void Herwig7HepMC3Interface::initRepository(const edm::ParameterSet &pset) {
0084   std::string runModeTemp = pset.getUntrackedParameter<string>("runModeList", "read,run");
0085
0086   // To Lower
0087   std::transform(runModeTemp.begin(), runModeTemp.end(), runModeTemp.begin(), ::tolower);
0088
0089   while (!runModeTemp.empty()) {
0090     // Split first part of List
0091     std::string choice;
0092     size_t pos = runModeTemp.find(',');
0093     if (std::string::npos == pos)
0094       choice = runModeTemp;
0095     else
0096       choice = runModeTemp.substr(0, pos);
0097
0098     if (pos == std::string::npos)
0099       runModeTemp.erase();
0100     else
0101       runModeTemp.erase(0, pos + 1);
0102
0103     HwUI_ = std::make_shared<Herwig::HerwigUIProvider>(pset, dumpConfig_, Herwig::RunMode::READ);
0104     edm::LogInfo("Herwig7HepMC3Interface")
0105         << "HerwigUIProvider object with run mode " << HwUI_->runMode() << " created.\n";
0106
0107     // Chose run mode
0108     if (choice == "read") {
0109       createInputFile(pset);
0110       HwUI_->setRunMode(Herwig::RunMode::READ, pset, dumpConfig_);
0111       edm::LogInfo("Herwig7HepMC3Interface")
0112           << "Input file " << dumpConfig_ << " will be passed to Herwig for the read step.\n";
0113       callHerwigGenerator();
0114     } else if (choice == "build") {
0115       createInputFile(pset);
0116       HwUI_->setRunMode(Herwig::RunMode::BUILD, pset, dumpConfig_);
0117       edm::LogInfo("Herwig7HepMC3Interface")
0118           << "Input file " << dumpConfig_ << " will be passed to Herwig for the build step.\n";
0119       callHerwigGenerator();
0120
0121     } else if (choice == "integrate") {
0122       std::string runFileName = run_ + ".run";
0123       edm::LogInfo("Herwig7HepMC3Interface")
0124           << "Run file " << runFileName << " will be passed to Herwig for the integrate step.\n";
0125       HwUI_->setRunMode(Herwig::RunMode::INTEGRATE, pset, runFileName);
0126       callHerwigGenerator();
0127
0128     } else if (choice == "run") {
0129       std::string runFileName = run_ + ".run";
0130       edm::LogInfo("Herwig7HepMC3Interface")
0131           << "Run file " << runFileName << " will be passed to Herwig for the run step.\n";
0132       HwUI_->setRunMode(Herwig::RunMode::RUN, pset, runFileName);
0133     } else {
0134       edm::LogInfo("Herwig7HepMC3Interface") << "Cannot recognize \"" << choice << "\".\n"
0135                                              << "Trying to skip step.\n";
0136       continue;
0137     }
0138   }
0139 }
0140
0141 void Herwig7HepMC3Interface::callHerwigGenerator() {
0142   try {
0143     edm::LogInfo("Herwig7HepMC3Interface")
0144         << "callHerwigGenerator function invoked with run mode " << HwUI_->runMode() << ".\n";
0145
0146     // Call program switches according to runMode
0147     switch (HwUI_->runMode()) {
0148       case Herwig::RunMode::INIT:
0149         Herwig::API::init(*HwUI_);
0150         break;
0151       case Herwig::RunMode::READ:
0152         Herwig::API::read(*HwUI_);
0153         break;
0154       case Herwig::RunMode::BUILD:
0155         Herwig::API::build(*HwUI_);
0156         break;
0157       case Herwig::RunMode::INTEGRATE:
0158         Herwig::API::integrate(*HwUI_);
0159         break;
0160       case Herwig::RunMode::MERGEGRIDS:
0161         Herwig::API::mergegrids(*HwUI_);
0162         break;
0163       case Herwig::RunMode::RUN: {
0164         HwUI_->setSeed(randomEngine->getSeed());
0165         eg_ = Herwig::API::prepareRun(*HwUI_);
0166         break;
0167       }
0168       case Herwig::RunMode::ERROR:
0169         edm::LogError("Herwig7HepMC3Interface") << "Error during read in of command line parameters.\n"
0170                                                 << "Program execution will stop now.";
0171         return;
0172       default:
0173         HwUI_->quitWithHelp();
0174     }
0175
0176     return;
0177
0178   } catch (ThePEG::Exception &e) {
0179     edm::LogError("Herwig7HepMC3Interface") << ": ThePEG::Exception caught.\n"
0180                                             << e.what() << '\n'
0181                                             << "See logfile for details.\n";
0182     return;
0183   } catch (std::exception &e) {
0184     edm::LogError("Herwig7HepMC3Interface") << ": " << e.what() << '\n';
0185     return;
0186   } catch (const char *what) {
0187     edm::LogError("Herwig7HepMC3Interface") << ": caught exception: " << what << "\n";
0188     return;
0189   }
0190 }
0191
0192 bool Herwig7HepMC3Interface::initGenerator() {
0193   if (HwUI_->runMode() == Herwig::RunMode::RUN) {
0194     edm::LogInfo("Herwig7HepMC3Interface") << "Starting EventGenerator initialization";
0195     callHerwigGenerator();
0196     edm::LogInfo("Herwig7HepMC3Interface") << "EventGenerator initialized";
0197
0198     // Skip events
0199     for (unsigned int i = 0; i < skipEvents_; i++) {
0200       flushRandomNumberGenerator();
0201       eg_->shoot();
0202       edm::LogInfo("Herwig7HepMC3Interface") << "Event discarded";
0203     }
0204
0205     return true;
0206
0207   } else {
0208     edm::LogInfo("Herwig7HepMC3Interface") << "Stopped EventGenerator due to missing run mode.";
0209     return false;
0210     /*
0211         throw cms::Exception("Herwig7HepMC3Interface")
0212             << "EventGenerator could not be initialized due to wrong run mode!" << endl;
0213 */
0214   }
0215 }
0216
0217 void Herwig7HepMC3Interface::flushRandomNumberGenerator() {
0218   /*ThePEG::RandomEngineGlue *rnd = randomEngineGlueProxy_->getInstance();
0219
0220     if (!rnd)
0221         edm::LogWarning("ProxyMissing")
0222             << "ThePEG not initialised with RandomEngineGlue.";
0223     else
0224         rnd->flush();
0225       */
0226 }
0227
0228 unique_ptr<HepMC3::GenEvent> Herwig7HepMC3Interface::convert(const ThePEG::EventPtr &event) {
0229   return std::unique_ptr<HepMC3::GenEvent>(ThePEG::HepMCConverter<HepMC3::GenEvent>::convert(*event));
0230 }
0231
0232 double Herwig7HepMC3Interface::pthat(const ThePEG::EventPtr &event) {
0233   using namespace ThePEG;
0234
0235   if (!event->primaryCollision())
0236     return -1.0;
0237
0238   tSubProPtr sub = event->primaryCollision()->primarySubProcess();
0239   TmpTransform<tSubProPtr> tmp(sub, Utilities::getBoostToCM(sub->incoming()));
0240
0241   double pthat = (*sub->outgoing().begin())->momentum().perp() / ThePEG::GeV;
0242   for (PVector::const_iterator it = sub->outgoing().begin(); it != sub->outgoing().end(); ++it)
0243     pthat = std::min<double>(pthat, (*it)->momentum().perp() / ThePEG::GeV);
0244
0245   return pthat;
0246 }
0247
0248 void Herwig7HepMC3Interface::createInputFile(const edm::ParameterSet &pset) {
0249   /* Initialize the input config for Herwig from
0250      * 1. the Herwig7 config files
0251      * 2. the CMSSW config blocks
0252      * Writes them to an output file which is read by Herwig
0253     */
0254
0255   stringstream logstream;
0256
0257   // Contains input config passed to Herwig
0258   stringstream herwiginputconfig;
0259
0260   // Define output file to which input config is written, too, if dumpConfig parameter is set.
0261   // Otherwise use default file HerwigConfig.in which is read in by Herwig
0262   ofstream cfgDump;
0263   cfgDump.open(dumpConfig_.c_str(), ios_base::app);
0264
0265   // Read Herwig config files as input
0266   vector<string> configFiles = pset.getParameter<vector<string> >("configFiles");
0267   // Loop over the config files
0268   for (const auto &iter : configFiles) {
0269     // Open external config file
0270     ifstream externalConfigFile(iter);
0271     if (externalConfigFile.is_open()) {
0272       edm::LogInfo("Herwig7HepMC3Interface") << "Reading config file (" << iter << ")" << endl;
0273       stringstream configFileStream;
0274       configFileStream << externalConfigFile.rdbuf();
0275       string configFileContent = configFileStream.str();
0276
0277       // Comment out occurence of saverun in config file since it is set later considering run and generator option
0278       string searchKeyword("saverun");
0279       if (configFileContent.find(searchKeyword) != std::string::npos) {
0280         edm::LogInfo("Herwig7HepMC3Interface")
0281             << "Commented out saverun command in external input config file(" << iter << ")" << endl;
0282         configFileContent.insert(configFileContent.find(searchKeyword), "#");
0283       }
0284       herwiginputconfig << "# Begin Config file input" << endl
0285                         << configFileContent << endl
0286                         << "# End Config file input";
0287       edm::LogInfo("Herwig7HepMC3Interface") << "Finished reading config file (" << iter << ")" << endl;
0288     } else {
0289       edm::LogWarning("Herwig7HepMC3Interface") << "Could not read config file (" << iter << ")" << endl;
0290     }
0291   }
0292
0293   edm::LogInfo("Herwig7HepMC3Interface") << "Start with processing CMSSW config" << endl;
0294   // Read CMSSW config file parameter sets starting from "parameterSets"
0295   ParameterCollector collector(pset);
0296   ParameterCollector::const_iterator iter;
0297   iter = collector.begin();
0298   herwiginputconfig << endl << "# Begin Parameter set input\n" << endl;
0299   for (; iter != collector.end(); ++iter) {
0300     herwiginputconfig << *iter << endl;
0301   }
0302
0303   // Add some additional necessary lines to the Herwig input config
0304   herwiginputconfig << "saverun " << run_ << " " << generator_ << endl;
0305   // write the ProxyID for the RandomEngineGlue to fill its pointer in
0306   ostringstream ss;
0307   ss << randomEngineGlueProxy_->getID();
0308   //herwiginputconfig << "set " << generator_ << ":RandomNumberGenerator:ProxyID " << ss.str() << endl;
0309
0310   // Dump Herwig input config to file, so that it can be read by Herwig
0311   cfgDump << herwiginputconfig.str() << endl;
0312   cfgDump.close();
0313 }