SherpaInterface/src/SherpaHadronizer.cc

0001 #include <cassert>
0002 #include <iostream>
0003 #include <sstream>
0004 #include <string>
0005 #include <memory>
0006 #include <cstdint>
0007 #include <vector>
0008
0009 #include "SHERPA/Main/Sherpa.H"
0010 #include "ATOOLS/Math/Random.H"
0011
0012 #include "ATOOLS/Org/Run_Parameter.H"
0013 #include "ATOOLS/Org/MyStrStream.H"
0014 #include "ATOOLS/Org/CXXFLAGS.H"
0015 #include "ATOOLS/Org/CXXFLAGS_PACKAGES.H"
0016 #include "ATOOLS/Org/My_MPI.H"
0017
0018 #include "GeneratorInterface/Core/interface/ParameterCollector.h"
0019 #include "GeneratorInterface/Core/interface/BaseHadronizer.h"
0020 #include "GeneratorInterface/Core/interface/GeneratorFilter.h"
0021 #include "GeneratorInterface/Core/interface/HadronizerFilter.h"
0022 #include "GeneratorInterface/SherpaInterface/interface/SherpackFetcher.h"
0023
0024 #include "CLHEP/Random/RandomEngine.h"
0025
0026 #include "FWCore/Utilities/interface/Exception.h"
0027
0028 //This unnamed namespace is used (instead of static variables) to pass the
0029 //randomEngine passed to doSetRandomEngine to the External Random
0030 //Number Generator CMS_SHERPA_RNG of sherpa
0031 //The advantage of the unnamed namespace over static variables is
0032 //that it is only accessible in this file
0033
0034 namespace {
0035   CLHEP::HepRandomEngine *ExternalEngine = nullptr;
0036   CLHEP::HepRandomEngine *GetExternalEngine() { return ExternalEngine; }
0037   void SetExternalEngine(CLHEP::HepRandomEngine *v) { ExternalEngine = v; }
0038 }  // namespace
0039
0040 class SherpaHadronizer : public gen::BaseHadronizer {
0041 public:
0042   SherpaHadronizer(const edm::ParameterSet &params);
0043   ~SherpaHadronizer() override;
0044
0045   bool readSettings(int) { return true; }
0046   bool initializeForInternalPartons();
0047   bool declareStableParticles(const std::vector<int> &pdgIds);
0048   bool declareSpecialSettings(const std::vector<std::string> &) { return true; }
0049   void statistics();
0050   bool generatePartonsAndHadronize();
0051   bool decay();
0052   bool rearrangeWeights;
0053   bool residualDecay();
0054   void finalizeEvent();
0055   std::unique_ptr<GenLumiInfoHeader> getGenLumiInfoHeader() const override;
0056   const char *classname() const { return "SherpaHadronizer"; }
0057
0058 private:
0059   void doSetRandomEngine(CLHEP::HepRandomEngine *v) override;
0060
0061   std::string SherpaProcess;
0062   std::string SherpaChecksum;
0063   std::string SherpaPath;
0064   std::string SherpaPathPiece;
0065   std::string SherpaResultDir;
0066   double SherpaDefaultWeight;
0067   edm::ParameterSet SherpaParameterSet;
0068   unsigned int maxEventsToPrint;
0069   std::vector<std::string> arguments;
0070   SHERPA::Sherpa *Generator = new SHERPA::Sherpa();
0071   bool isInitialized;
0072   bool isRNGinitialized;
0073   std::vector<std::string> weightlist;
0074   std::vector<std::string> variationweightlist;
0075 };
0076
0077 class CMS_SHERPA_RNG : public ATOOLS::External_RNG {
0078 public:
0079   CMS_SHERPA_RNG() : randomEngine(nullptr) {
0080     edm::LogVerbatim("SherpaHadronizer") << "Use stored reference for the external RNG";
0081     setRandomEngine(GetExternalEngine());
0082   }
0083   void setRandomEngine(CLHEP::HepRandomEngine *v) { randomEngine = v; }
0084
0085 private:
0086   double Get() override;
0087   CLHEP::HepRandomEngine *randomEngine;
0088 };
0089
0090 void SherpaHadronizer::doSetRandomEngine(CLHEP::HepRandomEngine *v) {
0091   CMS_SHERPA_RNG *cmsSherpaRng = dynamic_cast<CMS_SHERPA_RNG *>(ATOOLS::ran->GetExternalRng());
0092   if (cmsSherpaRng == nullptr) {
0093     //First time call to this function makes the interface store the reference in the unnamed namespace
0094     if (!isRNGinitialized) {
0095       isRNGinitialized = true;
0096       edm::LogVerbatim("SherpaHadronizer") << "Store assigned reference of the randomEngine";
0097       SetExternalEngine(v);
0098       // Throw exception if there is no reference to an external RNG and it is not the first call!
0099     } else {
0100       if (isInitialized and v != nullptr) {
0101         throw edm::Exception(edm::errors::LogicError)
0102             << "The Sherpa interface got a randomEngine reference but there is "
0103                "no reference to the external RNG to hand it over to\n";
0104       }
0105     }
0106   } else {
0107     cmsSherpaRng->setRandomEngine(v);
0108   }
0109 }
0110
0111 SherpaHadronizer::SherpaHadronizer(const edm::ParameterSet &params)
0112     : BaseHadronizer(params),
0113       SherpaParameterSet(params.getParameter<edm::ParameterSet>("SherpaParameters")),
0114       isRNGinitialized(false) {
0115   if (!params.exists("SherpaProcess"))
0116     SherpaProcess = "";
0117   else
0118     SherpaProcess = params.getParameter<std::string>("SherpaProcess");
0119   if (!params.exists("SherpaPath"))
0120     SherpaPath = "";
0121   else
0122     SherpaPath = params.getParameter<std::string>("SherpaPath");
0123   if (!params.exists("SherpaPathPiece"))
0124     SherpaPathPiece = "";
0125   else
0126     SherpaPathPiece = params.getParameter<std::string>("SherpaPathPiece");
0127   if (!params.exists("SherpaResultDir"))
0128     SherpaResultDir = "Result";
0129   else
0130     SherpaResultDir = params.getParameter<std::string>("SherpaResultDir");
0131   if (!params.exists("SherpaDefaultWeight"))
0132     SherpaDefaultWeight = 1.;
0133   else
0134     SherpaDefaultWeight = params.getParameter<double>("SherpaDefaultWeight");
0135   if (!params.exists("maxEventsToPrint"))
0136     maxEventsToPrint = 0;
0137   else
0138     maxEventsToPrint = params.getParameter<int>("maxEventsToPrint");
0139   // if hepmcextendedweights is used the event weights have to be reordered ( unordered list can be accessed via event->weights().write() )
0140   // two lists have to be provided:
0141   // 1) SherpaWeights
0142   // - containing nominal event weight, combined matrix element and phase space weight, event normalization, and possibly other sherpa weights
0143   // 2) SherpaVariationsWeights
0144   // - containing weights from scale and PDF variations ( have to be defined in the runcard )
0145   // - in case of unweighted events these weights are also divided by the event normalization (see list 1 )
0146   // Sherpa Documentation: http://sherpa.hepforge.org/doc/SHERPA-MC-2.2.0.html#Scale-and-PDF-variations
0147   if (!params.exists("SherpaWeightsBlock")) {
0148     rearrangeWeights = false;
0149   } else {
0150     rearrangeWeights = true;
0151     edm::ParameterSet WeightsBlock = params.getParameter<edm::ParameterSet>("SherpaWeightsBlock");
0152     if (WeightsBlock.exists("SherpaWeights"))
0153       weightlist = WeightsBlock.getParameter<std::vector<std::string> >("SherpaWeights");
0154     else
0155       throw cms::Exception("SherpaInterface") << "SherpaWeights does not exists in SherpaWeightsBlock" << std::endl;
0156     if (WeightsBlock.exists("SherpaVariationWeights"))
0157       variationweightlist = WeightsBlock.getParameter<std::vector<std::string> >("SherpaVariationWeights");
0158     else
0159       throw cms::Exception("SherpaInterface")
0160           << "SherpaVariationWeights does not exists in SherpaWeightsBlock" << std::endl;
0161     edm::LogVerbatim("SherpaHadronizer") << "SherpaHadronizer will try rearrange the event weights according to "
0162                                             "SherpaWeights and SherpaVariationWeights";
0163   }
0164
0165   spf::SherpackFetcher Fetcher(params);
0166   int retval = Fetcher.Fetch();
0167   if (retval != 0) {
0168     throw cms::Exception("SherpaInterface") << "SherpaHadronizer: Preparation of Sherpack failed ... ";
0169   }
0170   // The ids (names) of parameter sets to be read (Analysis,Run) to create Analysis.dat, Run.dat
0171   //They are given as a vstring.
0172   std::vector<std::string> setNames = SherpaParameterSet.getParameter<std::vector<std::string> >("parameterSets");
0173   //Loop all set names...
0174   for (unsigned i = 0; i < setNames.size(); ++i) {
0175     // ...and read the parameters for each set given in vstrings
0176     std::vector<std::string> pars = SherpaParameterSet.getParameter<std::vector<std::string> >(setNames[i]);
0177     edm::LogVerbatim("SherpaHadronizer") << "Write Sherpa parameter set " << setNames[i] << " to " << setNames[i]
0178                                          << ".dat ";
0179     std::string datfile = SherpaPath + "/" + setNames[i] + ".dat";
0180     std::ofstream os(datfile.c_str());
0181     // Loop over all strings and write the according *.dat
0182     for (std::vector<std::string>::const_iterator itPar = pars.begin(); itPar != pars.end(); ++itPar) {
0183       os << (*itPar) << std::endl;
0184     }
0185   }
0186
0187   //To be conform to the default Sherpa usage create a command line:
0188   //name of executable  (only for demonstration, could also be empty)
0189   std::string shRun = "./Sherpa";
0190   //Path where the Sherpa libraries are stored
0191   std::string shPath = "PATH=" + SherpaPath;
0192   // new for Sherpa 1.3.0, suggested by authors
0193   std::string shPathPiece = "PATH_PIECE=" + SherpaPathPiece;
0194   //Path where results are stored
0195   std::string shRes = "RESULT_DIRECTORY=" + SherpaResultDir;  // from Sherpa 1.2.0 on
0196   //Name of the external random number class
0197   std::string shRng = "EXTERNAL_RNG=CMS_SHERPA_RNG";
0198
0199   //create the command line
0200   arguments.push_back(shRun);
0201   arguments.push_back(shPath);
0202   arguments.push_back(shPathPiece);
0203   arguments.push_back(shRes);
0204   arguments.push_back(shRng);
0205   isInitialized = false;
0206   //initialization of Sherpa moved to initializeForInternalPartons
0207 #ifdef USING__MPI
0208   MPI::Init();
0209 #endif
0210 }
0211
0212 SherpaHadronizer::~SherpaHadronizer() {
0213   Generator->~Sherpa();
0214 #ifdef USING__MPI
0215   MPI::Finalize();
0216 #endif
0217 }
0218
0219 bool SherpaHadronizer::initializeForInternalPartons() {
0220   //initialize Sherpa but only once
0221   if (!isInitialized) {
0222     int argc = arguments.size();
0223     char *argv[argc];
0224     for (int l = 0; l < argc; l++)
0225       argv[l] = (char *)arguments[l].c_str();
0226     Generator->InitializeTheRun(argc, argv);
0227     Generator->InitializeTheEventHandler();
0228     isInitialized = true;
0229   }
0230   return true;
0231 }
0232
0233 bool SherpaHadronizer::declareStableParticles(const std::vector<int> &pdgIds) { return false; }
0234
0235 void SherpaHadronizer::statistics() {
0236   //calculate statistics
0237   Generator->SummarizeRun();
0238
0239   //get the xsec & err
0240   double xsec_val = Generator->TotalXS();
0241   double xsec_err = Generator->TotalErr();
0242
0243   //set the internal cross section in pb in GenRunInfoProduct
0244   runInfo().setInternalXSec(GenRunInfoProduct::XSec(xsec_val, xsec_err));
0245 }
0246
0247 bool SherpaHadronizer::generatePartonsAndHadronize() {
0248   //get the next event and check if it produced
0249   bool rc = false;
0250   int itry = 0;
0251   bool gen_event = true;
0252   while ((itry < 3) && gen_event) {
0253     try {
0254       rc = Generator->GenerateOneEvent();
0255       gen_event = false;
0256     } catch (...) {
0257       ++itry;
0258       std::cerr << "Exception from Generator->GenerateOneEvent() catch. Call # " << itry << " for this event\n";
0259     }
0260   }
0261   if (rc) {
0262     //convert it to HepMC2
0263     auto evt = std::make_unique<HepMC::GenEvent>();
0264     Generator->FillHepMCEvent(*evt);
0265
0266     // in case of unweighted events sherpa puts the max weight as event weight.
0267     // this is not optimal, we want 1 for unweighted events, so we check
0268     // whether we are producing unweighted events ("EVENT_GENERATION_MODE" == "1")
0269     // the information about the weights to the HepMC weight vector:
0270     // [0] event weight
0271     // [1] combined matrix element and phase space weight (missing only PDF information, thus directly suitable for PDF reweighting)
0272     // [2] event weight normalisation (in case of unweighted events event weights of ~ +/-1 can be obtained by (event weight)/(event weight normalisation))
0273     // [3] number of trials.
0274     // see also: https://sherpa.hepforge.org/doc/SHERPA-MC-2.1.0.html#Event-output-formats
0275     bool unweighted = false;
0276     double weight_normalization = -1;
0277     if (ATOOLS::ToType<int>(ATOOLS::rpa->gen.Variable("EVENT_GENERATION_MODE")) == 1) {
0278       if (evt->weights().size() > 2) {
0279         unweighted = true;
0280         weight_normalization = evt->weights()[2];
0281       } else {
0282         throw cms::Exception("SherpaInterface")
0283             << "Requested unweighted production. Missing normalization weight." << std::endl;
0284       }
0285     }
0286
0287     // vector to fill new weights in correct order
0288     std::vector<double> newWeights;
0289     if (rearrangeWeights) {
0290       for (auto &i : weightlist) {
0291         if (evt->weights().has_key(i)) {
0292           newWeights.push_back(evt->weights()[i]);
0293         } else {
0294           throw cms::Exception("SherpaInterface")
0295               << "Missing weights! Key " << i << " not found, please check the weight definition!" << std::endl;
0296         }
0297       }
0298       for (auto &i : variationweightlist) {
0299         if (evt->weights().has_key(i)) {
0300           if (unweighted) {
0301             newWeights.push_back(evt->weights()[i] / weight_normalization);
0302           } else {
0303             newWeights.push_back(evt->weights()[i]);
0304           }
0305         } else {
0306           throw cms::Exception("SherpaInterface")
0307               << "Missing weights! Key " << i << " not found, please check the weight definition!" << std::endl;
0308         }
0309       }
0310
0311       //Change original weights for reordered ones
0312       evt->weights().clear();
0313       for (auto &elem : newWeights) {
0314         evt->weights().push_back(elem);
0315       }
0316     }
0317
0318     if (unweighted) {
0319       evt->weights()[0] /= weight_normalization;
0320     }
0321     resetEvent(std::move(evt));
0322     return true;
0323   } else {
0324     return false;
0325   }
0326 }
0327
0328 bool SherpaHadronizer::decay() { return true; }
0329
0330 bool SherpaHadronizer::residualDecay() { return true; }
0331
0332 void SherpaHadronizer::finalizeEvent() {
0333   //******** Verbosity *******
0334   if (maxEventsToPrint > 0) {
0335     maxEventsToPrint--;
0336     event()->print();
0337   }
0338 }
0339
0340 //GETTER for the external random numbers
0341 DECLARE_GETTER(CMS_SHERPA_RNG, "CMS_SHERPA_RNG", ATOOLS::External_RNG, ATOOLS::RNG_Key);
0342
0343 ATOOLS::External_RNG *ATOOLS::Getter<ATOOLS::External_RNG, ATOOLS::RNG_Key, CMS_SHERPA_RNG>::operator()(
0344     const ATOOLS::RNG_Key &) const {
0345   return new CMS_SHERPA_RNG();
0346 }
0347
0348 void ATOOLS::Getter<ATOOLS::External_RNG, ATOOLS::RNG_Key, CMS_SHERPA_RNG>::PrintInfo(std::ostream &str,
0349                                                                                       const size_t) const {
0350   str << "CMS_SHERPA_RNG interface";
0351 }
0352
0353 double CMS_SHERPA_RNG::Get() {
0354   if (randomEngine == nullptr) {
0355     throw edm::Exception(edm::errors::LogicError) << "The Sherpa code attempted to a generate random number while\n"
0356                                                   << "the engine pointer was null. This might mean that the code\n"
0357                                                   << "was modified to generate a random number outside the event and\n"
0358                                                   << "beginLuminosityBlock methods, which is not allowed.\n";
0359   }
0360   return randomEngine->flat();
0361 }
0362 std::unique_ptr<GenLumiInfoHeader> SherpaHadronizer::getGenLumiInfoHeader() const {
0363   auto genLumiInfoHeader = BaseHadronizer::getGenLumiInfoHeader();
0364
0365   if (rearrangeWeights) {
0366     edm::LogPrint("SherpaHadronizer") << "The order of event weights was changed!";
0367     for (auto &i : weightlist) {
0368       genLumiInfoHeader->weightNames().push_back(i);
0369       edm::LogVerbatim("SherpaHadronizer") << i;
0370     }
0371     for (auto &i : variationweightlist) {
0372       genLumiInfoHeader->weightNames().push_back(i);
0373       edm::LogVerbatim("SherpaHadronizer") << i;
0374     }
0375   }
0376   return genLumiInfoHeader;
0377 }
0378 #include "GeneratorInterface/ExternalDecays/interface/ExternalDecayDriver.h"
0379
0380 typedef edm::GeneratorFilter<SherpaHadronizer, gen::ExternalDecayDriver> SherpaGeneratorFilter;
0381 DEFINE_FWK_MODULE(SherpaGeneratorFilter);