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