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   // These should never fail
0058   assert(cmsswbase);
0059   assert(cmsswrelease);
0060   if (!std::getenv("RIVET_REF_PATH")) {
0061     const std::string rivetref = string(cmsswbase) +
0062                                  "/src/GeneratorInterface/RivetInterface/data:" + string(cmsswrelease) +
0063                                  "/src/GeneratorInterface/RivetInterface/data:.";
0064     char* rivetrefCstr = strdup(rivetref.c_str());
0065     setenv("RIVET_REF_PATH", rivetrefCstr, 1);
0066     free(rivetrefCstr);
0067   }
0068   if (!std::getenv("RIVET_INFO_PATH")) {
0069     const std::string rivetinfo = string(cmsswbase) +
0070                                   "/src/GeneratorInterface/RivetInterface/data:" + string(cmsswrelease) +
0071                                   "/src/GeneratorInterface/RivetInterface/data:.";
0072     char* rivetinfoCstr = strdup(rivetinfo.c_str());
0073     setenv("RIVET_INFO_PATH", rivetinfoCstr, 1);
0074     free(rivetinfoCstr);
0075   }
0076 }
0077
0078 void RivetAnalyzer::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0079   if (_useLHEweights) {
0080     edm::Handle<LHERunInfoProduct> lheRunInfoHandle;
0081     iRun.getByLabel(_lheLabel, lheRunInfoHandle);
0082     typedef std::vector<LHERunInfoProduct::Header>::const_iterator headers_const_iterator;
0083
0084     std::regex reg("<weight.*> ?(.*?) ?<\\/weight>");
0085     for (headers_const_iterator iter = lheRunInfoHandle->headers_begin(); iter != lheRunInfoHandle->headers_end();
0086          iter++) {
0087       std::vector<std::string> lines = iter->lines();
0088       for (unsigned int iLine = 0; iLine < lines.size(); iLine++) {
0089         std::smatch match;
0090         std::regex_search(lines.at(iLine), match, reg);
0091         if (!match.empty()) {
0092           _lheWeightNames.push_back(match[1]);
0093         }
0094       }
0095     }
0096   }
0097 }
0098
0099 void RivetAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0100   //finalize weight names on the first event
0101   if (_isFirstEvent) {
0102     auto genLumiInfoHandle = iEvent.getLuminosityBlock().getHandle(_genLumiInfoToken);
0103     if (genLumiInfoHandle.isValid()) {
0104       _weightNames = genLumiInfoHandle->weightNames();
0105     }
0106
0107     // need to reset the default weight name (or plotting will fail)
0108     if (!_weightNames.empty()) {
0109       _weightNames[0] = "";
0110     } else {  // Summer16 samples have 1 weight stored in HepMC but no weightNames
0111       _weightNames.push_back("");
0112     }
0113     if (_useLHEweights) {
0114       // Some samples have weights but no weight names -> assign generic names lheN
0115       if (_lheWeightNames.empty()) {
0116         edm::Handle<LHEEventProduct> lheEventHandle;
0117         iEvent.getByToken(_LHECollection, lheEventHandle);
0118         for (unsigned int i = 0; i < lheEventHandle->weights().size(); i++) {
0119           _lheWeightNames.push_back("lhe" + std::to_string(i + 1));  // start with 1 to match LHE weight IDs
0120         }
0121       }
0122       _weightNames.insert(_weightNames.end(), _lheWeightNames.begin(), _lheWeightNames.end());
0123     }
0124     // clean weight names to be accepted by Rivet plotting
0125     for (const std::string& wn : _weightNames) {
0126       _cleanedWeightNames.push_back(std::regex_replace(wn, std::regex("[^A-Za-z\\d\\._=]"), "_"));
0127     }
0128     runinfo = make_shared<HepMC3::GenRunInfo>();
0129     runinfo->set_weight_names(_cleanedWeightNames);
0130   }
0131
0132   //get the hepmc product from the event
0133   edm::Handle<HepMC3Product> evt;
0134   iEvent.getByToken(_hepmcCollection, evt);
0135
0136   // get HepMC GenEvent
0137   const HepMC3::GenEventData* genEventData = evt->GetEvent();
0138   std::unique_ptr<HepMC3::GenEvent> genEvent = std::make_unique<HepMC3::GenEvent>();
0139   genEvent->read_data(*genEventData);
0140
0141   std::vector<double> mergedWeights;
0142   for (unsigned int i = 0; i < genEvent->weights().size(); i++) {
0143     mergedWeights.push_back(genEvent->weights()[i]);
0144   }
0145
0146   if (_useLHEweights) {
0147     edm::Handle<LHEEventProduct> lheEventHandle;
0148     iEvent.getByToken(_LHECollection, lheEventHandle);
0149     for (unsigned int i = 0; i < _lheWeightNames.size(); i++) {
0150       mergedWeights.push_back(genEvent->weights()[0] * lheEventHandle->weights().at(i).wgt /
0151                               lheEventHandle->originalXWGTUP());
0152     }
0153   }
0154
0155   double xsection = _xsection > 0 ? _xsection : genEvent->cross_section()->xsecs()[0];
0156   HepMC3::GenCrossSectionPtr xsec = make_shared<HepMC3::GenCrossSection>();
0157   xsec->set_cross_section(std::vector<double>(mergedWeights.size(), xsection),
0158                           std::vector<double>(mergedWeights.size(), 0.));
0159   genEvent->set_cross_section(xsec);
0160   genEvent->set_run_info(runinfo);
0161   genEvent->weights() = mergedWeights;
0162
0163   //apply the beams initialization on the first event
0164   if (_isFirstEvent) {
0165     _analysisHandler = std::make_unique<Rivet::AnalysisHandler>();
0166     _analysisHandler->addAnalyses(_analysisNames);
0167
0168     /// Set analysis handler weight options
0169     _analysisHandler->setCheckBeams(!_setIgnoreBeams);
0170     _analysisHandler->skipMultiWeights(_skipMultiWeights);
0171     _analysisHandler->matchWeightNames(_selectMultiWeights);
0172     _analysisHandler->unmatchWeightNames(_deselectMultiWeights);
0173     _analysisHandler->setNominalWeightName(_setNominalWeightName);
0174     _analysisHandler->setWeightCap(_weightCap);
0175     _analysisHandler->setNLOSmearing(_NLOSmearing);
0176
0177     _analysisHandler->init(*genEvent);
0178
0179     _isFirstEvent = false;
0180   }
0181
0182   //run the analysis
0183   _analysisHandler->analyze(const_cast<GenEvent&>(*genEvent));
0184 }
0185
0186 void RivetAnalyzer::endRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0187   if (_doFinalize)
0188     _analysisHandler->finalize();
0189   _analysisHandler->writeData(_outFileName);
0190
0191   return;
0192 }
0193
0194 void RivetAnalyzer::endJob() {}
0195
0196 DEFINE_FWK_MODULE(RivetAnalyzer);