diff --git a/GeneratorInterface/SherpaInterface/data/PrepareSherpaLibs.sh b/GeneratorInterface/SherpaInterface/data/PrepareSherpaLibs.sh index 732d33b43000c..47c9f713fcd47 100755 --- a/GeneratorInterface/SherpaInterface/data/PrepareSherpaLibs.sh +++ b/GeneratorInterface/SherpaInterface/data/PrepareSherpaLibs.sh @@ -31,6 +31,8 @@ function print_help() { echo " -L filename (optional) name of library file ( "${cflb}" )" && \ echo " -C filename (optional) name of cross section file ( "${cfcr}" )" && \ echo " -G filename (optional) name of MI grid file ( "${cfgr}" )" && \ + echo " -e filename (optional) name of extended weight list file ( "${weightfile}" ) " && \ + echo " (example file in GeneratorInterface/SherpaInterface/python/ExtendedSherpaWeights_cfi.py)" && \ ## echo " -P SRM path (CRAB) SE path for final results" && \ ## echo " -> ( "${MYSRMPATH}" )" && \ echo " -h display this help and exit" && echo @@ -52,6 +54,9 @@ function build_python_cff() { echo "import FWCore.ParameterSet.Config as cms" >> ${cfffilename} echo "import os" >> ${cfffilename} + if [ ! $weightlist == "" ]; then + echo "from $weightlist import *" >> ${cfffilename} + fi echo "" >> ${cfffilename} echo "source = cms.Source(\"EmptySource\")" >> ${cfffilename} echo "" >> ${cfffilename} @@ -72,6 +77,9 @@ function build_python_cff() { ## fi echo " SherpaResultDir = cms.string('Result')," >> ${cfffilename} echo " SherpaDefaultWeight = cms.double(1.0)," >> ${cfffilename} + if [ ! $weightlist == "" ]; then + echo " SherpaWeightsBlock = SherpaWeightsBlock," >> ${cfffilename} + fi echo " SherpaParameters = cms.PSet(parameterSets = cms.vstring(" >> ${cfffilename} fcnt=0 for file in `ls *.dat`; do @@ -204,6 +212,15 @@ function file_copy() { } +function parse_extended_weight() { + if [ -e $CMSSW_BASE/src/$weightlist ] + then + weightlist=$(echo "$weightlist" | sed 's:/python/:.:g' | sed 's:/:.:g' | sed 's:.py::g' ) + else + echo '$CMSSW_BASE/src/$weightlist does not exist.' + exit 1 + fi +} @@ -248,13 +265,14 @@ cfdc="" # custom data card file nam cflb="" # custom library file name cfcr="" # custom cross section file name cfgr="" # custom MI grid file name +weightfile="" # sherpa weights ordering file ##MYSRMPATH="./" # SRM path for storage of results TDIR=TMP # get & evaluate options ##while getopts :i:p:d:m:a:D:L:C:G:P:h OPT -while getopts :i:p:d:D:L:C:G:h OPT +while getopts :i:p:d:D:L:C:G:e:h OPT do case $OPT in i) datadir=$OPTARG ;; @@ -265,6 +283,7 @@ do L) cflb=$OPTARG ;; C) cfcr=$OPTARG ;; G) cfgr=$OPTARG ;; + e) weightlist=$OPTARG && parse_extended_weight;; ## P) MYSRMPATH=$OPTARG ;; h) print_help && exit 0 ;; \?) diff --git a/GeneratorInterface/SherpaInterface/python/ExtendedSherpaWeights_cfi.py b/GeneratorInterface/SherpaInterface/python/ExtendedSherpaWeights_cfi.py new file mode 100644 index 0000000000000..7bf0665db4b8e --- /dev/null +++ b/GeneratorInterface/SherpaInterface/python/ExtendedSherpaWeights_cfi.py @@ -0,0 +1,229 @@ +import FWCore.ParameterSet.Config as cms + +SherpaWeightsBlock = cms.PSet( + SherpaWeights = cms.vstring( + 'Weight', + 'MEWeight', + 'WeightNormalisation', + 'NTrials' + ), + SherpaVariationWeights = cms.vstring( + 'MUR0.5_MUF0.5_PDF260000', + 'MUR0.5_MUF1_PDF260000', + 'MUR0.5_MUF2_PDF260000', + 'MUR1_MUF0.5_PDF260000', + 'MUR1_MUF2_PDF260000', + 'MUR2_MUF0.5_PDF260000', + 'MUR2_MUF1_PDF260000', + 'MUR2_MUF2_PDF260000', + 'MUR1_MUF1_PDF260000', + 'MUR1_MUF1_PDF260001', + 'MUR1_MUF1_PDF260002', + 'MUR1_MUF1_PDF260003', + 'MUR1_MUF1_PDF260004', + 'MUR1_MUF1_PDF260005', + 'MUR1_MUF1_PDF260006', + 'MUR1_MUF1_PDF260007', + 'MUR1_MUF1_PDF260008', + 'MUR1_MUF1_PDF260009', + 'MUR1_MUF1_PDF260010', + 'MUR1_MUF1_PDF260011', + 'MUR1_MUF1_PDF260012', + 'MUR1_MUF1_PDF260013', + 'MUR1_MUF1_PDF260014', + 'MUR1_MUF1_PDF260015', + 'MUR1_MUF1_PDF260016', + 'MUR1_MUF1_PDF260017', + 'MUR1_MUF1_PDF260018', + 'MUR1_MUF1_PDF260019', + 'MUR1_MUF1_PDF260020', + 'MUR1_MUF1_PDF260021', + 'MUR1_MUF1_PDF260022', + 'MUR1_MUF1_PDF260023', + 'MUR1_MUF1_PDF260024', + 'MUR1_MUF1_PDF260025', + 'MUR1_MUF1_PDF260026', + 'MUR1_MUF1_PDF260027', + 'MUR1_MUF1_PDF260028', + 'MUR1_MUF1_PDF260029', + 'MUR1_MUF1_PDF260030', + 'MUR1_MUF1_PDF260031', + 'MUR1_MUF1_PDF260032', + 'MUR1_MUF1_PDF260033', + 'MUR1_MUF1_PDF260034', + 'MUR1_MUF1_PDF260035', + 'MUR1_MUF1_PDF260036', + 'MUR1_MUF1_PDF260037', + 'MUR1_MUF1_PDF260038', + 'MUR1_MUF1_PDF260039', + 'MUR1_MUF1_PDF260040', + 'MUR1_MUF1_PDF260041', + 'MUR1_MUF1_PDF260042', + 'MUR1_MUF1_PDF260043', + 'MUR1_MUF1_PDF260044', + 'MUR1_MUF1_PDF260045', + 'MUR1_MUF1_PDF260046', + 'MUR1_MUF1_PDF260047', + 'MUR1_MUF1_PDF260048', + 'MUR1_MUF1_PDF260049', + 'MUR1_MUF1_PDF260050', + 'MUR1_MUF1_PDF260051', + 'MUR1_MUF1_PDF260052', + 'MUR1_MUF1_PDF260053', + 'MUR1_MUF1_PDF260054', + 'MUR1_MUF1_PDF260055', + 'MUR1_MUF1_PDF260056', + 'MUR1_MUF1_PDF260057', + 'MUR1_MUF1_PDF260058', + 'MUR1_MUF1_PDF260059', + 'MUR1_MUF1_PDF260060', + 'MUR1_MUF1_PDF260061', + 'MUR1_MUF1_PDF260062', + 'MUR1_MUF1_PDF260063', + 'MUR1_MUF1_PDF260064', + 'MUR1_MUF1_PDF260065', + 'MUR1_MUF1_PDF260066', + 'MUR1_MUF1_PDF260067', + 'MUR1_MUF1_PDF260068', + 'MUR1_MUF1_PDF260069', + 'MUR1_MUF1_PDF260070', + 'MUR1_MUF1_PDF260071', + 'MUR1_MUF1_PDF260072', + 'MUR1_MUF1_PDF260073', + 'MUR1_MUF1_PDF260074', + 'MUR1_MUF1_PDF260075', + 'MUR1_MUF1_PDF260076', + 'MUR1_MUF1_PDF260077', + 'MUR1_MUF1_PDF260078', + 'MUR1_MUF1_PDF260079', + 'MUR1_MUF1_PDF260080', + 'MUR1_MUF1_PDF260081', + 'MUR1_MUF1_PDF260082', + 'MUR1_MUF1_PDF260083', + 'MUR1_MUF1_PDF260084', + 'MUR1_MUF1_PDF260085', + 'MUR1_MUF1_PDF260086', + 'MUR1_MUF1_PDF260087', + 'MUR1_MUF1_PDF260088', + 'MUR1_MUF1_PDF260089', + 'MUR1_MUF1_PDF260090', + 'MUR1_MUF1_PDF260091', + 'MUR1_MUF1_PDF260092', + 'MUR1_MUF1_PDF260093', + 'MUR1_MUF1_PDF260094', + 'MUR1_MUF1_PDF260095', + 'MUR1_MUF1_PDF260096', + 'MUR1_MUF1_PDF260097', + 'MUR1_MUF1_PDF260098', + 'MUR1_MUF1_PDF260099', + 'MUR1_MUF1_PDF260100', + 'MUR1_MUF1_PDF13100', + 'MUR1_MUF1_PDF13101', + 'MUR1_MUF1_PDF13102', + 'MUR1_MUF1_PDF13103', + 'MUR1_MUF1_PDF13104', + 'MUR1_MUF1_PDF13105', + 'MUR1_MUF1_PDF13106', + 'MUR1_MUF1_PDF13107', + 'MUR1_MUF1_PDF13108', + 'MUR1_MUF1_PDF13109', + 'MUR1_MUF1_PDF13110', + 'MUR1_MUF1_PDF13111', + 'MUR1_MUF1_PDF13112', + 'MUR1_MUF1_PDF13113', + 'MUR1_MUF1_PDF13114', + 'MUR1_MUF1_PDF13115', + 'MUR1_MUF1_PDF13116', + 'MUR1_MUF1_PDF13117', + 'MUR1_MUF1_PDF13118', + 'MUR1_MUF1_PDF13119', + 'MUR1_MUF1_PDF13120', + 'MUR1_MUF1_PDF13121', + 'MUR1_MUF1_PDF13122', + 'MUR1_MUF1_PDF13123', + 'MUR1_MUF1_PDF13124', + 'MUR1_MUF1_PDF13125', + 'MUR1_MUF1_PDF13126', + 'MUR1_MUF1_PDF13127', + 'MUR1_MUF1_PDF13128', + 'MUR1_MUF1_PDF13129', + 'MUR1_MUF1_PDF13130', + 'MUR1_MUF1_PDF13131', + 'MUR1_MUF1_PDF13132', + 'MUR1_MUF1_PDF13133', + 'MUR1_MUF1_PDF13134', + 'MUR1_MUF1_PDF13135', + 'MUR1_MUF1_PDF13136', + 'MUR1_MUF1_PDF13137', + 'MUR1_MUF1_PDF13138', + 'MUR1_MUF1_PDF13139', + 'MUR1_MUF1_PDF13140', + 'MUR1_MUF1_PDF13141', + 'MUR1_MUF1_PDF13142', + 'MUR1_MUF1_PDF13143', + 'MUR1_MUF1_PDF13144', + 'MUR1_MUF1_PDF13145', + 'MUR1_MUF1_PDF13146', + 'MUR1_MUF1_PDF13147', + 'MUR1_MUF1_PDF13148', + 'MUR1_MUF1_PDF13149', + 'MUR1_MUF1_PDF13150', + 'MUR1_MUF1_PDF13151', + 'MUR1_MUF1_PDF13152', + 'MUR1_MUF1_PDF13153', + 'MUR1_MUF1_PDF13154', + 'MUR1_MUF1_PDF13155', + 'MUR1_MUF1_PDF13156', + 'MUR1_MUF1_PDF25100', + 'MUR1_MUF1_PDF25101', + 'MUR1_MUF1_PDF25102', + 'MUR1_MUF1_PDF25103', + 'MUR1_MUF1_PDF25104', + 'MUR1_MUF1_PDF25105', + 'MUR1_MUF1_PDF25106', + 'MUR1_MUF1_PDF25107', + 'MUR1_MUF1_PDF25108', + 'MUR1_MUF1_PDF25109', + 'MUR1_MUF1_PDF25110', + 'MUR1_MUF1_PDF25111', + 'MUR1_MUF1_PDF25112', + 'MUR1_MUF1_PDF25113', + 'MUR1_MUF1_PDF25114', + 'MUR1_MUF1_PDF25115', + 'MUR1_MUF1_PDF25116', + 'MUR1_MUF1_PDF25117', + 'MUR1_MUF1_PDF25118', + 'MUR1_MUF1_PDF25119', + 'MUR1_MUF1_PDF25120', + 'MUR1_MUF1_PDF25121', + 'MUR1_MUF1_PDF25122', + 'MUR1_MUF1_PDF25123', + 'MUR1_MUF1_PDF25124', + 'MUR1_MUF1_PDF25125', + 'MUR1_MUF1_PDF25126', + 'MUR1_MUF1_PDF25127', + 'MUR1_MUF1_PDF25128', + 'MUR1_MUF1_PDF25129', + 'MUR1_MUF1_PDF25130', + 'MUR1_MUF1_PDF25131', + 'MUR1_MUF1_PDF25132', + 'MUR1_MUF1_PDF25133', + 'MUR1_MUF1_PDF25134', + 'MUR1_MUF1_PDF25135', + 'MUR1_MUF1_PDF25136', + 'MUR1_MUF1_PDF25137', + 'MUR1_MUF1_PDF25138', + 'MUR1_MUF1_PDF25139', + 'MUR1_MUF1_PDF25140', + 'MUR1_MUF1_PDF25141', + 'MUR1_MUF1_PDF25142', + 'MUR1_MUF1_PDF25143', + 'MUR1_MUF1_PDF25144', + 'MUR1_MUF1_PDF25145', + 'MUR1_MUF1_PDF25146', + 'MUR1_MUF1_PDF25147', + 'MUR1_MUF1_PDF25148', + 'MUR1_MUF1_PDF25149', + 'MUR1_MUF1_PDF25150' + ) +) diff --git a/GeneratorInterface/SherpaInterface/src/SherpaHadronizer.cc b/GeneratorInterface/SherpaInterface/src/SherpaHadronizer.cc index 595afefd7803b..38fc1995bf729 100644 --- a/GeneratorInterface/SherpaInterface/src/SherpaHadronizer.cc +++ b/GeneratorInterface/SherpaInterface/src/SherpaHadronizer.cc @@ -22,10 +22,10 @@ #include "CLHEP/Random/RandomEngine.h" -//This unnamed namespace is used (instead of static variables) to pass the +//This unnamed namespace is used (instead of static variables) to pass the //randomEngine passed to doSetRandomEngine to the External Random //Number Generator CMS_SHERPA_RNG of sherpa -//The advantage of the unnamed namespace over static variables is +//The advantage of the unnamed namespace over static variables is //that it is only accessible in this file namespace { @@ -46,15 +46,17 @@ class SherpaHadronizer : public gen::BaseHadronizer { void statistics(); bool generatePartonsAndHadronize(); bool decay(); + bool rearrangeWeights; bool residualDecay(); void finalizeEvent(); + GenLumiInfoHeader *getGenLumiInfoHeader() const override; const char *classname() const { return "SherpaHadronizer"; } - + private: virtual void doSetRandomEngine(CLHEP::HepRandomEngine* v) override; - + std::string SherpaProcess; std::string SherpaChecksum; std::string SherpaPath; @@ -67,6 +69,8 @@ class SherpaHadronizer : public gen::BaseHadronizer { SHERPA::Sherpa Generator; bool isInitialized; bool isRNGinitialized; + std::vector weightlist; + std::vector variationweightlist; }; @@ -76,12 +80,12 @@ class CMS_SHERPA_RNG: public ATOOLS::External_RNG { public: CMS_SHERPA_RNG() : randomEngine(nullptr) { - std::cout << "Use stored reference for the external RNG" << std::endl; - setRandomEngine(GetExternalEngine()); + edm::LogVerbatim("SherpaHadronizer") << "Use stored reference for the external RNG"; + setRandomEngine(GetExternalEngine()); } void setRandomEngine(CLHEP::HepRandomEngine* v) { randomEngine = v; } - -private: + +private: double Get() override; CLHEP::HepRandomEngine* randomEngine; }; @@ -94,16 +98,16 @@ void SherpaHadronizer::doSetRandomEngine(CLHEP::HepRandomEngine* v) { //First time call to this function makes the interface store the reference in the unnamed namespace if (!isRNGinitialized){ isRNGinitialized=true; - std::cout << "Store assigned reference of the randomEngine" << std::endl; + edm::LogVerbatim("SherpaHadronizer") << "Store assigned reference of the randomEngine"; SetExternalEngine(v); // Throw exception if there is no reference to an external RNG and it is not the first call! } else { - throw edm::Exception(edm::errors::LogicError) + throw edm::Exception(edm::errors::LogicError) << "The Sherpa interface got a randomEngine reference but there is no reference to the external RNG to hand it over to\n"; } } else { cmsSherpaRng->setRandomEngine(v); - } + } } SherpaHadronizer::SherpaHadronizer(const edm::ParameterSet ¶ms) : @@ -123,15 +127,35 @@ SherpaHadronizer::SherpaHadronizer(const edm::ParameterSet ¶ms) : else SherpaDefaultWeight=params.getParameter("SherpaDefaultWeight"); if (!params.exists("maxEventsToPrint")) maxEventsToPrint=0; else maxEventsToPrint=params.getParameter("maxEventsToPrint"); +// if hepmcextendedweights is used the event weights have to be reordered ( unordered list can be accessed via event->weights().write() ) +// two lists have to be provided: +// 1) SherpaWeights +// - containing nominal event weight, combined matrix element and phase space weight, event normalization, and possibly other sherpa weights +// 2) SherpaVariationsWeights +// - containing weights from scale and PDF variations ( have to be defined in the runcard ) +// - in case of unweighted events these weights are also divided by the event normalization (see list 1 ) +// Sherpa Documentation: http://sherpa.hepforge.org/doc/SHERPA-MC-2.2.0.html#Scale-and-PDF-variations + if (!params.exists("SherpaWeightsBlock")) { + rearrangeWeights=false; + } else { + rearrangeWeights=true; + edm::ParameterSet WeightsBlock = params.getParameter("SherpaWeightsBlock"); + if (WeightsBlock.exists("SherpaWeights")) + weightlist=WeightsBlock.getParameter< std::vector >("SherpaWeights"); + else + throw cms::Exception("SherpaInterface") <<"SherpaWeights does not exists in SherpaWeightsBlock" << std::endl; + if (WeightsBlock.exists("SherpaVariationWeights")) + variationweightlist=WeightsBlock.getParameter< std::vector >("SherpaVariationWeights"); + else + throw cms::Exception("SherpaInterface") <<"SherpaVariationWeights does not exists in SherpaWeightsBlock" << std::endl; + edm::LogVerbatim("SherpaHadronizer") << "SherpaHadronizer will try rearrange the event weights according to SherpaWeights and SherpaVariationWeights"; + } spf::SherpackFetcher Fetcher(params); int retval=Fetcher.Fetch(); if (retval != 0) { - std::cout << "SherpaHadronizer: Preparation of Sherpack failed ... " << std::endl; - std::cout << "SherpaHadronizer: Error code: " << retval << std::endl; - std::terminate(); - + throw cms::Exception("SherpaInterface") <<"SherpaHadronizer: Preparation of Sherpack failed ... "; } // The ids (names) of parameter sets to be read (Analysis,Run) to create Analysis.dat, Run.dat //They are given as a vstring. @@ -140,7 +164,7 @@ SherpaHadronizer::SherpaHadronizer(const edm::ParameterSet ¶ms) : for ( unsigned i=0; i pars = SherpaParameterSet.getParameter >(setNames[i]); - std::cout << "Write Sherpa parameter set " << setNames[i] <<" to "<(ATOOLS::rpa->gen.Variable("EVENT_GENERATION_MODE")) == 1){ if (evt->weights().size()>2) { - evt->weights()[0]/=evt->weights()[2]; + unweighted = true; + weight_normalization = evt->weights()[2]; + }else{ + throw cms::Exception("SherpaInterface") <<"Requested unweighted production. Missing normalization weight." << std::endl; } } + + // vector to fill new weights in correct order + std::vector newWeights; + if (rearrangeWeights){ + for ( auto &i : weightlist ) { + if (evt->weights().has_key(i)) { + newWeights.push_back(evt->weights()[i]); + } else { + throw cms::Exception("SherpaInterface") <<"Missing weights! Key " << i << " not found, please check the weight definition!" << std::endl; + } + } + for ( auto &i : variationweightlist ) { + if (evt->weights().has_key(i)) { + if(unweighted){ + newWeights.push_back(evt->weights()[i]/weight_normalization); + }else{ + newWeights.push_back(evt->weights()[i]); + } + } else { + throw cms::Exception("SherpaInterface") <<"Missing weights! Key " << i << " not found, please check the weight definition!" << std::endl; + } + + } + } + + //Change original weights for reordered ones + evt->weights().clear(); + for (auto& elem: newWeights) { + evt->weights().push_back(elem); + } + + if(unweighted){ + evt->weights()[0]/=weight_normalization; + } resetEvent(evt); return true; } @@ -314,9 +379,24 @@ double CMS_SHERPA_RNG::Get() { << "beginLuminosityBlock methods, which is not allowed.\n"; } return randomEngine->flat(); + +} +GenLumiInfoHeader *SherpaHadronizer::getGenLumiInfoHeader() const { + GenLumiInfoHeader *genLumiInfoHeader = BaseHadronizer::getGenLumiInfoHeader(); + if(rearrangeWeights){ + edm::LogPrint("SherpaHadronizer") << "The order of event weights was changed!" ; + for(auto &i: weightlist){ + genLumiInfoHeader->weightNames().push_back(i); + edm::LogVerbatim("SherpaHadronizer") << i; + } + for(auto &i: variationweightlist) { + genLumiInfoHeader->weightNames().push_back(i); + edm::LogVerbatim("SherpaHadronizer") << i; + } + } + return genLumiInfoHeader; } - #include "GeneratorInterface/ExternalDecays/interface/ExternalDecayDriver.h" typedef edm::GeneratorFilter SherpaGeneratorFilter;