RivetInterface/plugins/RivetAnalyzer.cc

0001 #include "GeneratorInterface/RivetInterface/interface/RivetAnalyzer.h"
0002
0003 #include "FWCore/Framework/interface/Event.h"
0004 #include "FWCore/Framework/interface/MakerMacros.h"
0005 #include "FWCore/Framework/interface/Run.h"
0006
0007 #include "DataFormats/Common/interface/Handle.h"
0008 #include "FWCore/ServiceRegistry/interface/Service.h"
0009
0010 #include "Rivet/Run.hh"
0011 #include "Rivet/AnalysisHandler.hh"
0012 #include "Rivet/Analysis.hh"
0013
0014 #include <regex>
0015
0016 using namespace Rivet;
0017 using namespace edm;
0018
0019 RivetAnalyzer::RivetAnalyzer(const edm::ParameterSet& pset)
0020     : _isFirstEvent(true),
0021       _outFileName(pset.getParameter<std::string>("OutputFile")),
0022       _analysisNames(pset.getParameter<std::vector<std::string> >("AnalysisNames")),
0023       //decide whether to finalize the plots or not.
0024       //deciding not to finalize them can be useful for further harvesting of many jobs
0025       _doFinalize(pset.getParameter<bool>("DoFinalize")),
0026       _lheLabel(pset.getParameter<edm::InputTag>("LHECollection")),
0027       _xsection(-1.) {
0028   usesResource("Rivet");
0029
0030   _hepmcCollection = consumes<HepMC3Product>(pset.getParameter<edm::InputTag>("HepMCCollection"));
0031   _genLumiInfoToken = consumes<GenLumiInfoHeader, edm::InLumi>(pset.getParameter<edm::InputTag>("genLumiInfo"));
0032
0033   _useLHEweights = pset.getParameter<bool>("useLHEweights");
0034   if (_useLHEweights) {
0035     _lheRunInfoToken = consumes<LHERunInfoProduct, edm::InRun>(_lheLabel);
0036     _LHECollection = consumes<LHEEventProduct>(_lheLabel);
0037   }
0038
0039   _weightCap = pset.getParameter<double>("weightCap");
0040   _NLOSmearing = pset.getParameter<double>("NLOSmearing");
0041   _setIgnoreBeams = pset.getParameter<bool>("setIgnoreBeams");
0042   _skipMultiWeights = pset.getParameter<bool>("skipMultiWeights");
0043   _selectMultiWeights = pset.getParameter<std::string>("selectMultiWeights");
0044   _deselectMultiWeights = pset.getParameter<std::string>("deselectMultiWeights");
0045   _setNominalWeightName = pset.getParameter<std::string>("setNominalWeightName");
0046
0047   //set user cross section if needed
0048   _xsection = pset.getParameter<double>("CrossSection");
0049 }
0050
0051 RivetAnalyzer::~RivetAnalyzer() {}
0052
0053 void RivetAnalyzer::beginJob() {
0054   //set the environment, very ugly but rivet is monolithic when it comes to paths
0055   char* cmsswbase = std::getenv("CMSSW_BASE");
0056   char* cmsswrelease = std::getenv("CMSSW_RELEASE_BASE");
0057   if (!std::getenv("RIVET_REF_PATH")) {
0058     const std::string rivetref = string(cmsswbase) +
0059                                  "/src/GeneratorInterface/RivetInterface/data:" + string(cmsswrelease) +
0060                                  "/src/GeneratorInterface/RivetInterface/data:.";
0061     char* rivetrefCstr = strdup(rivetref.c_str());
0062     setenv("RIVET_REF_PATH", rivetrefCstr, 1);
0063     free(rivetrefCstr);
0064   }
0065   if (!std::getenv("RIVET_INFO_PATH")) {
0066     const std::string rivetinfo = string(cmsswbase) +
0067                                   "/src/GeneratorInterface/RivetInterface/data:" + string(cmsswrelease) +
0068                                   "/src/GeneratorInterface/RivetInterface/data:.";
0069     char* rivetinfoCstr = strdup(rivetinfo.c_str());
0070     setenv("RIVET_INFO_PATH", rivetinfoCstr, 1);
0071     free(rivetinfoCstr);
0072   }
0073 }
0074
0075 void RivetAnalyzer::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0076   if (_useLHEweights) {
0077     edm::Handle<LHERunInfoProduct> lheRunInfoHandle;
0078     iRun.getByLabel(_lheLabel, lheRunInfoHandle);
0079     typedef std::vector<LHERunInfoProduct::Header>::const_iterator headers_const_iterator;
0080
0081     std::regex reg("<weight.*> ?(.*?) ?<\\/weight>");
0082     for (headers_const_iterator iter = lheRunInfoHandle->headers_begin(); iter != lheRunInfoHandle->headers_end();
0083          iter++) {
0084       std::vector<std::string> lines = iter->lines();
0085       for (unsigned int iLine = 0; iLine < lines.size(); iLine++) {
0086         std::smatch match;
0087         std::regex_search(lines.at(iLine), match, reg);
0088         if (!match.empty()) {
0089           _lheWeightNames.push_back(match[1]);
0090         }
0091       }
0092     }
0093   }
0094 }
0095
0096 void RivetAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0097   //finalize weight names on the first event
0098   if (_isFirstEvent) {
0099     auto genLumiInfoHandle = iEvent.getLuminosityBlock().getHandle(_genLumiInfoToken);
0100     if (genLumiInfoHandle.isValid()) {
0101       _weightNames = genLumiInfoHandle->weightNames();
0102     }
0103
0104     // need to reset the default weight name (or plotting will fail)
0105     if (!_weightNames.empty()) {
0106       _weightNames[0] = "";
0107     } else {  // Summer16 samples have 1 weight stored in HepMC but no weightNames
0108       _weightNames.push_back("");
0109     }
0110     if (_useLHEweights) {
0111       // Some samples have weights but no weight names -> assign generic names lheN
0112       if (_lheWeightNames.empty()) {
0113         edm::Handle<LHEEventProduct> lheEventHandle;
0114         iEvent.getByToken(_LHECollection, lheEventHandle);
0115         for (unsigned int i = 0; i < lheEventHandle->weights().size(); i++) {
0116           _lheWeightNames.push_back("lhe" + std::to_string(i + 1));  // start with 1 to match LHE weight IDs
0117         }
0118       }
0119       _weightNames.insert(_weightNames.end(), _lheWeightNames.begin(), _lheWeightNames.end());
0120     }
0121     // clean weight names to be accepted by Rivet plotting
0122     for (const std::string& wn : _weightNames) {
0123       _cleanedWeightNames.push_back(std::regex_replace(wn, std::regex("[^A-Za-z\\d\\._=]"), "_"));
0124     }
0125     runinfo = make_shared<HepMC3::GenRunInfo>();
0126     runinfo->set_weight_names(_cleanedWeightNames);
0127   }
0128
0129   //get the hepmc product from the event
0130   edm::Handle<HepMC3Product> evt;
0131   iEvent.getByToken(_hepmcCollection, evt);
0132
0133   // get HepMC GenEvent
0134   const HepMC3::GenEventData* genEventData = evt->GetEvent();
0135   std::unique_ptr<HepMC3::GenEvent> genEvent = std::make_unique<HepMC3::GenEvent>();
0136   genEvent->read_data(*genEventData);
0137
0138   std::vector<double> mergedWeights;
0139   for (unsigned int i = 0; i < genEvent->weights().size(); i++) {
0140     mergedWeights.push_back(genEvent->weights()[i]);
0141   }
0142
0143   if (_useLHEweights) {
0144     edm::Handle<LHEEventProduct> lheEventHandle;
0145     iEvent.getByToken(_LHECollection, lheEventHandle);
0146     for (unsigned int i = 0; i < _lheWeightNames.size(); i++) {
0147       mergedWeights.push_back(genEvent->weights()[0] * lheEventHandle->weights().at(i).wgt /
0148                               lheEventHandle->originalXWGTUP());
0149     }
0150   }
0151
0152   double xsection = _xsection > 0 ? _xsection : genEvent->cross_section()->xsecs()[0];
0153   HepMC3::GenCrossSectionPtr xsec = make_shared<HepMC3::GenCrossSection>();
0154   xsec->set_cross_section(std::vector<double>(mergedWeights.size(), xsection),
0155                           std::vector<double>(mergedWeights.size(), 0.));
0156   genEvent->set_cross_section(xsec);
0157   genEvent->set_run_info(runinfo);
0158   genEvent->weights() = mergedWeights;
0159
0160   //apply the beams initialization on the first event
0161   if (_isFirstEvent) {
0162     _analysisHandler = std::make_unique<Rivet::AnalysisHandler>();
0163     _analysisHandler->addAnalyses(_analysisNames);
0164
0165     /// Set analysis handler weight options
0166     _analysisHandler->setCheckBeams(!_setIgnoreBeams);
0167     _analysisHandler->skipMultiWeights(_skipMultiWeights);
0168     _analysisHandler->matchWeightNames(_selectMultiWeights);
0169     _analysisHandler->unmatchWeightNames(_deselectMultiWeights);
0170     _analysisHandler->setNominalWeightName(_setNominalWeightName);
0171     _analysisHandler->setWeightCap(_weightCap);
0172     _analysisHandler->setNLOSmearing(_NLOSmearing);
0173
0174     _analysisHandler->init(*genEvent);
0175
0176     _isFirstEvent = false;
0177   }
0178
0179   //run the analysis
0180   _analysisHandler->analyze(const_cast<GenEvent&>(*genEvent));
0181 }
0182
0183 void RivetAnalyzer::endRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0184   if (_doFinalize)
0185     _analysisHandler->finalize();
0186   _analysisHandler->writeData(_outFileName);
0187
0188   return;
0189 }
0190
0191 void RivetAnalyzer::endJob() {}
0192
0193 DEFINE_FWK_MODULE(RivetAnalyzer);