Herwig7Hadronizer.cc

Go to the documentation of this file.

#include <memory>
#include <sstream>
#include <fstream>

#include <HepMC/GenEvent.h>
#include <HepMC/IO_BaseClass.h>

#include <ThePEG/Repository/Repository.h>
#include <ThePEG/EventRecord/Event.h>
#include <ThePEG/Config/ThePEG.h>
#include <ThePEG/LesHouches/LesHouchesReader.h>

#include "FWCore/Framework/interface/LuminosityBlock.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
#include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"

#include "GeneratorInterface/Core/interface/BaseHadronizer.h"
#include "GeneratorInterface/Core/interface/GeneratorFilter.h"
#include "GeneratorInterface/Core/interface/HadronizerFilter.h"

#include "GeneratorInterface/LHEInterface/interface/LHEProxy.h"

#include "GeneratorInterface/Herwig7Interface/interface/Herwig7Interface.h"

#include <Herwig/API/HerwigAPI.h>
#include "CLHEP/Random/RandomEngine.h"

namespace CLHEP {
  class HepRandomEngine;
}

class Herwig7Hadronizer : public Herwig7Interface, public gen::BaseHadronizer {
public:
  Herwig7Hadronizer(const edm::ParameterSet& params);
  ~Herwig7Hadronizer() override;

  bool readSettings(int) { return true; }
  bool initializeForInternalPartons();
  bool initializeForExternalPartons();
  bool declareStableParticles(const std::vector<int>& pdgIds);
  bool declareSpecialSettings(const std::vector<std::string>) { return true; }

  void statistics();

  bool generatePartonsAndHadronize();
  bool hadronize();
  bool decay();
  bool residualDecay();
  void finalizeEvent();

  const char* classname() const { return "Herwig7Hadronizer"; }
  std::unique_ptr<GenLumiInfoHeader> getGenLumiInfoHeader() const override;
  void beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&);
  void randomizeIndex(edm::LuminosityBlock const& lumi, CLHEP::HepRandomEngine* rengine);

private:
  void doSetRandomEngine(CLHEP::HepRandomEngine* v) override { setPEGRandomEngine(v); }

  unsigned int eventsToPrint;

  ThePEG::EventPtr thepegEvent;

  std::shared_ptr<lhef::LHEProxy> proxy_;
  const std::string handlerDirectory_;
  edm::ParameterSet paramSettings;
  const std::string runFileName;

  unsigned int firstLumiBlock = 0;
  unsigned int currentLumiBlock = 0;
};

Herwig7Hadronizer::Herwig7Hadronizer(const edm::ParameterSet& pset)
    : Herwig7Interface(pset),
      BaseHadronizer(pset),
      eventsToPrint(pset.getUntrackedParameter<unsigned int>("eventsToPrint", 0)),
      handlerDirectory_(pset.getParameter<std::string>("eventHandlers")),
      runFileName(pset.getParameter<std::string>("run")) {
  initRepository(pset);
  paramSettings = pset;
}

Herwig7Hadronizer::~Herwig7Hadronizer() {}

bool Herwig7Hadronizer::initializeForInternalPartons() {
  if (currentLumiBlock == firstLumiBlock) {
    std::ifstream runFile(runFileName + ".run");
    if (runFile.fail())  //required for showering of LHE files
    {
      initRepository(paramSettings);
    }
    if (!initGenerator()) {
      edm::LogInfo("Generator|Herwig7Hadronizer") << "No run step for Herwig chosen. Program will be aborted.";
      exit(0);
    }
  }
  return true;
}

bool Herwig7Hadronizer::initializeForExternalPartons() {
  edm::LogError("Herwig7 interface")
      << "Read in of LHE files is not supported in this way. You can read them manually if necessary.";
  return false;
}

bool Herwig7Hadronizer::declareStableParticles(const std::vector<int>& pdgIds) { return false; }

void Herwig7Hadronizer::statistics() {
  if (eg_) {
    runInfo().setInternalXSec(
        GenRunInfoProduct::XSec(eg_->integratedXSec() / ThePEG::picobarn, eg_->integratedXSecErr() / ThePEG::picobarn));
  }
}

bool Herwig7Hadronizer::generatePartonsAndHadronize() {
  edm::LogInfo("Generator|Herwig7Hadronizer") << "Start production";

  try {
    thepegEvent = eg_->shoot();
  } catch (std::exception& exc) {
    edm::LogWarning("Generator|Herwig7Hadronizer")
        << "EGPtr::shoot() thrown an exception, event skipped: " << exc.what();
    return false;
  }

  if (!thepegEvent) {
    edm::LogWarning("Generator|Herwig7Hadronizer") << "thepegEvent not initialized";
    return false;
  }

  event() = convert(thepegEvent);
  if (!event().get()) {
    edm::LogWarning("Generator|Herwig7Hadronizer") << "genEvent not initialized";
    return false;
  }

  return true;
}

bool Herwig7Hadronizer::hadronize() {
  edm::LogError("Herwig7 interface")
      << "Read in of LHE files is not supported in this way. You can read them manually if necessary.";
  return false;
}

void Herwig7Hadronizer::finalizeEvent() {
  eventInfo() = std::make_unique<GenEventInfoProduct>(event().get());
  eventInfo()->setBinningValues(std::vector<double>(1, pthat(thepegEvent)));

  if (eventsToPrint) {
    eventsToPrint--;
    event()->print();
  }

  if (iobc_.get())
    iobc_->write_event(event().get());

  edm::LogInfo("Generator|Herwig7Hadronizer") << "Event produced";
}

bool Herwig7Hadronizer::decay() { return true; }

bool Herwig7Hadronizer::residualDecay() { return true; }

std::unique_ptr<GenLumiInfoHeader> Herwig7Hadronizer::getGenLumiInfoHeader() const {
  auto genLumiInfoHeader = BaseHadronizer::getGenLumiInfoHeader();

  if (thepegEvent) {
    int weights_number = thepegEvent->optionalWeights().size();

    if (weights_number > 1) {
      genLumiInfoHeader->weightNames().reserve(weights_number + 1);
      genLumiInfoHeader->weightNames().push_back("nominal");
      std::map<std::string, double> weights_map = thepegEvent->optionalWeights();
      for (std::map<std::string, double>::iterator it = weights_map.begin(); it != weights_map.end(); it++) {
        genLumiInfoHeader->weightNames().push_back(it->first);
      }
    }
  }

  return genLumiInfoHeader;
}

void Herwig7Hadronizer::randomizeIndex(edm::LuminosityBlock const& lumi, CLHEP::HepRandomEngine* rengine) {
  BaseHadronizer::randomizeIndex(lumi, rengine);

  if (firstLumiBlock == 0) {
    firstLumiBlock = lumi.id().luminosityBlock();
  }
  currentLumiBlock = lumi.id().luminosityBlock();
}

#include "GeneratorInterface/ExternalDecays/interface/ExternalDecayDriver.h"

typedef edm::GeneratorFilter<Herwig7Hadronizer, gen::ExternalDecayDriver> Herwig7GeneratorFilter;
DEFINE_FWK_MODULE(Herwig7GeneratorFilter);

typedef edm::HadronizerFilter<Herwig7Hadronizer, gen::ExternalDecayDriver> Herwig7HadronizerFilter;
DEFINE_FWK_MODULE(Herwig7HadronizerFilter);