Herwig7Interface/plugins/Herwig7HepMC3Hadronizer.cc

0001 #include <memory>
0002 #include <sstream>
0003 #include <fstream>
0004
0005 #include <HepMC3/GenEvent.h>
0006 //#include <HepMC/IO_BaseClass.h>
0007 #include "HepMC3/Print.h"
0008
0009 #include <ThePEG/Repository/Repository.h>
0010 #include <ThePEG/EventRecord/Event.h>
0011 #include <ThePEG/Config/ThePEG.h>
0012 #include <ThePEG/LesHouches/LesHouchesReader.h>
0013
0014 #include "FWCore/Framework/interface/LuminosityBlock.h"
0015 #include "FWCore/Framework/interface/EventSetup.h"
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017
0018 #include "SimDataFormats/GeneratorProducts/interface/HepMC3Product.h"
0019 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0020 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct3.h"
0021 #include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
0022
0023 #include "GeneratorInterface/Core/interface/BaseHadronizer.h"
0024 #include "GeneratorInterface/Core/interface/GeneratorFilter.h"
0025 #include "GeneratorInterface/Core/interface/HadronizerFilter.h"
0026
0027 #include "GeneratorInterface/LHEInterface/interface/LHEProxy.h"
0028
0029 #include "GeneratorInterface/Herwig7Interface/interface/Herwig7HepMC3Interface.h"
0030
0031 #include <Herwig/API/HerwigAPI.h>
0032 #include "CLHEP/Random/RandomEngine.h"
0033
0034 namespace CLHEP {
0035   class HepRandomEngine;
0036 }
0037
0038 class Herwig7HepMC3Hadronizer : public Herwig7HepMC3Interface, public gen::BaseHadronizer {
0039 public:
0040   Herwig7HepMC3Hadronizer(const edm::ParameterSet& params);
0041   ~Herwig7HepMC3Hadronizer() override;
0042
0043   bool readSettings(int) { return true; }
0044   bool initializeForInternalPartons();
0045   bool initializeForExternalPartons();
0046   bool declareStableParticles(const std::vector<int>& pdgIds);
0047   bool declareSpecialSettings(const std::vector<std::string>) { return true; }
0048
0049   void statistics();
0050
0051   bool generatePartonsAndHadronize();
0052   bool hadronize();
0053   bool decay();
0054   bool residualDecay();
0055   void finalizeEvent();
0056
0057   const char* classname() const { return "Herwig7HepMC3Hadronizer"; }
0058   std::unique_ptr<GenLumiInfoHeader> getGenLumiInfoHeader() const override;
0059   void beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&);
0060   void randomizeIndex(edm::LuminosityBlock const& lumi, CLHEP::HepRandomEngine* rengine);
0061
0062 private:
0063   void doSetRandomEngine(CLHEP::HepRandomEngine* v) override { setPEGRandomEngine(v); }
0064
0065   unsigned int eventsToPrint;
0066
0067   ThePEG::EventPtr thepegEvent;
0068   bool haveEvt = false;
0069
0070   std::shared_ptr<lhef::LHEProxy> proxy_;
0071   const std::string handlerDirectory_;
0072   edm::ParameterSet paramSettings;
0073   const std::string runFileName;
0074
0075   unsigned int firstLumiBlock = 0;
0076   unsigned int currentLumiBlock = 0;
0077 };
0078
0079 Herwig7HepMC3Hadronizer::Herwig7HepMC3Hadronizer(const edm::ParameterSet& pset)
0080     : Herwig7HepMC3Interface(pset),
0081       BaseHadronizer(pset),
0082       eventsToPrint(pset.getUntrackedParameter<unsigned int>("eventsToPrint", 0)),
0083       handlerDirectory_(pset.getParameter<std::string>("eventHandlers")),
0084       runFileName(pset.getParameter<std::string>("run")) {
0085   ivhepmc = 3;
0086   initRepository(pset);
0087   paramSettings = pset;
0088 }
0089
0090 Herwig7HepMC3Hadronizer::~Herwig7HepMC3Hadronizer() {}
0091
0092 bool Herwig7HepMC3Hadronizer::initializeForInternalPartons() {
0093   if (currentLumiBlock == firstLumiBlock) {
0094     std::ifstream runFile(runFileName + ".run");
0095     if (runFile.fail())  //required for showering of LHE files
0096     {
0097       initRepository(paramSettings);
0098     }
0099     if (!initGenerator()) {
0100       edm::LogInfo("Generator|Herwig7HepMC3Hadronizer") << "No run step for Herwig chosen. Program will be aborted.";
0101       exit(0);
0102     }
0103   }
0104   return true;
0105 }
0106
0107 bool Herwig7HepMC3Hadronizer::initializeForExternalPartons() {
0108   if (currentLumiBlock == firstLumiBlock) {
0109     std::ifstream runFile(runFileName + ".run");
0110     if (runFile.fail())  //required for showering of LHE files
0111     {
0112       initRepository(paramSettings);
0113     }
0114     if (!initGenerator()) {
0115       edm::LogInfo("Generator|Herwig7HepMC3Hadronizer") << "No run step for Herwig chosen. Program will be aborted.";
0116       exit(0);
0117     }
0118   }
0119   return true;
0120 }
0121
0122 bool Herwig7HepMC3Hadronizer::declareStableParticles(const std::vector<int>& pdgIds) { return false; }
0123
0124 void Herwig7HepMC3Hadronizer::statistics() {
0125   if (eg_) {
0126     runInfo().setInternalXSec(
0127         GenRunInfoProduct::XSec(eg_->integratedXSec() / ThePEG::picobarn, eg_->integratedXSecErr() / ThePEG::picobarn));
0128   }
0129 }
0130
0131 bool Herwig7HepMC3Hadronizer::generatePartonsAndHadronize() {
0132   edm::LogInfo("Generator|Herwig7HepMC3Hadronizer") << "Start production";
0133
0134   try {
0135     thepegEvent = eg_->shoot();
0136   } catch (std::exception& exc) {
0137     edm::LogWarning("Generator|Herwig7HepMC3Hadronizer")
0138         << "EGPtr::shoot() thrown an exception, event skipped: " << exc.what();
0139     return false;
0140   }
0141
0142   if (!thepegEvent) {
0143     edm::LogWarning("Generator|Herwig7HepMC3Hadronizer") << "thepegEvent not initialized";
0144     return false;
0145   }
0146
0147   event3() = convert(thepegEvent);
0148   if (!event3().get()) {
0149     edm::LogWarning("Generator|Herwig7HepMC3Hadronizer") << "genEvent not initialized";
0150     return false;
0151   }
0152
0153   return true;
0154 }
0155
0156 bool Herwig7HepMC3Hadronizer::hadronize() {
0157   if (!haveEvt) {
0158     try {
0159       thepegEvent = eg_->shoot();
0160       haveEvt = true;
0161     } catch (std::exception& exc) {
0162       edm::LogWarning("Generator|Herwig7HepMC3Hadronizer")
0163           << "EGPtr::shoot() thrown an exception, event skipped: " << exc.what();
0164       return false;
0165     }
0166   }
0167   int evtnum = lheEvent()->evtnum();
0168   if (evtnum == -1) {
0169     edm::LogError("Generator|Herwig7HepMC3Hadronizer")
0170         << "Event number not set in lhe file, needed for correctly aligning Herwig and LHE events!";
0171     return false;
0172   }
0173   if (thepegEvent->number() < evtnum) {
0174     edm::LogError("Herwig7 interface") << "Herwig does not seem to be generating events in order, did you set "
0175                                           "/Herwig/EventHandlers/FxFxLHReader:AllowedToReOpen Yes?";
0176     return false;
0177   } else if (thepegEvent->number() == evtnum) {
0178     haveEvt = false;
0179     if (!thepegEvent) {
0180       edm::LogWarning("Generator|Herwig7HepMC3Hadronizer") << "thepegEvent not initialized";
0181       return false;
0182     }
0183
0184     event3() = convert(thepegEvent);
0185     if (!event3().get()) {
0186       edm::LogWarning("Generator|Herwig7HepMC3Hadronizer") << "genEvent not initialized";
0187       return false;
0188     }
0189     return true;
0190   }
0191   edm::LogWarning("Generator|Herwig7HepMC3Hadronizer")
0192       << "Event " << evtnum << " not generated (likely skipped in merging)";
0193   return false;
0194 }
0195
0196 void Herwig7HepMC3Hadronizer::finalizeEvent() {
0197   eventInfo3() = std::make_unique<GenEventInfoProduct3>(event3().get());
0198   eventInfo3()->setBinningValues(std::vector<double>(1, pthat(thepegEvent)));
0199
0200   if (eventsToPrint) {
0201     eventsToPrint--;
0202     //event3()->print();
0203     HepMC3::Print::listing(*(event3().get()));
0204   }
0205
0206   //  if (iobc_.get())
0207   //    iobc_->write_event(event().get());
0208
0209   edm::LogInfo("Generator|Herwig7HepMC3Hadronizer") << "Event produced";
0210 }
0211
0212 bool Herwig7HepMC3Hadronizer::decay() { return true; }
0213
0214 bool Herwig7HepMC3Hadronizer::residualDecay() { return true; }
0215
0216 std::unique_ptr<GenLumiInfoHeader> Herwig7HepMC3Hadronizer::getGenLumiInfoHeader() const {
0217   auto genLumiInfoHeader = BaseHadronizer::getGenLumiInfoHeader();
0218
0219   if (thepegEvent) {
0220     int weights_number = thepegEvent->optionalWeights().size();
0221
0222     if (weights_number > 1) {
0223       genLumiInfoHeader->weightNames().reserve(weights_number + 1);
0224       genLumiInfoHeader->weightNames().push_back("nominal");
0225       std::map<std::string, double> weights_map = thepegEvent->optionalWeights();
0226       for (std::map<std::string, double>::iterator it = weights_map.begin(); it != weights_map.end(); it++) {
0227         genLumiInfoHeader->weightNames().push_back(it->first);
0228       }
0229     }
0230   }
0231
0232   return genLumiInfoHeader;
0233 }
0234
0235 void Herwig7HepMC3Hadronizer::randomizeIndex(edm::LuminosityBlock const& lumi, CLHEP::HepRandomEngine* rengine) {
0236   BaseHadronizer::randomizeIndex(lumi, rengine);
0237
0238   if (firstLumiBlock == 0) {
0239     firstLumiBlock = lumi.id().luminosityBlock();
0240   }
0241   currentLumiBlock = lumi.id().luminosityBlock();
0242 }
0243
0244 #include "GeneratorInterface/ExternalDecays/interface/ExternalDecayDriver.h"
0245
0246 typedef edm::GeneratorFilter<Herwig7HepMC3Hadronizer, gen::ExternalDecayDriver> Herwig7HepMC3GeneratorFilter;
0247 DEFINE_FWK_MODULE(Herwig7HepMC3GeneratorFilter);
0248
0249 typedef edm::HadronizerFilter<Herwig7HepMC3Hadronizer, gen::ExternalDecayDriver> Herwig7HepMC3HadronizerFilter;
0250 DEFINE_FWK_MODULE(Herwig7HepMC3HadronizerFilter);