Herwig7Interface/plugins/Herwig7Hadronizer.cc

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