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