Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 09:48:58

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<HepMCProduct>(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   }
0126 
0127   //get the hepmc product from the event
0128   edm::Handle<HepMCProduct> evt;
0129   iEvent.getByToken(_hepmcCollection, evt);
0130 
0131   // get HepMC GenEvent
0132   const HepMC::GenEvent* myGenEvent = evt->GetEvent();
0133   std::unique_ptr<HepMC::GenEvent> tmpGenEvtPtr;
0134   //if you want to use an external weight or set the cross section we have to clone the GenEvent and change the weight
0135   tmpGenEvtPtr = std::make_unique<HepMC::GenEvent>(*(evt->GetEvent()));
0136 
0137   if (_xsection > 0) {
0138     HepMC::GenCrossSection xsec;
0139     xsec.set_cross_section(_xsection);
0140     tmpGenEvtPtr->set_cross_section(xsec);
0141   }
0142 
0143   std::vector<double> mergedWeights;
0144   for (unsigned int i = 0; i < tmpGenEvtPtr->weights().size(); i++) {
0145     mergedWeights.push_back(tmpGenEvtPtr->weights()[i]);
0146   }
0147 
0148   if (_useLHEweights) {
0149     edm::Handle<LHEEventProduct> lheEventHandle;
0150     iEvent.getByToken(_LHECollection, lheEventHandle);
0151     for (unsigned int i = 0; i < _lheWeightNames.size(); i++) {
0152       mergedWeights.push_back(tmpGenEvtPtr->weights()[0] * lheEventHandle->weights().at(i).wgt /
0153                               lheEventHandle->originalXWGTUP());
0154     }
0155   }
0156 
0157   tmpGenEvtPtr->weights().clear();
0158   for (unsigned int i = 0; i < _cleanedWeightNames.size(); i++) {
0159     tmpGenEvtPtr->weights()[_cleanedWeightNames[i]] = mergedWeights[i];
0160   }
0161   myGenEvent = tmpGenEvtPtr.get();
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->setIgnoreBeams(_setIgnoreBeams);
0170     _analysisHandler->skipMultiWeights(_skipMultiWeights);
0171     _analysisHandler->selectMultiWeights(_selectMultiWeights);
0172     _analysisHandler->deselectMultiWeights(_deselectMultiWeights);
0173     _analysisHandler->setNominalWeightName(_setNominalWeightName);
0174     _analysisHandler->setWeightCap(_weightCap);
0175     _analysisHandler->setNLOSmearing(_NLOSmearing);
0176 
0177     _analysisHandler->init(*myGenEvent);
0178 
0179     _isFirstEvent = false;
0180   }
0181 
0182   //run the analysis
0183   _analysisHandler->analyze(*myGenEvent);
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);