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       _produceDQM(pset.getParameter<bool>("ProduceDQMOutput")),
0027       _lheLabel(pset.getParameter<edm::InputTag>("LHECollection")),
0028       _xsection(-1.) {
0029   usesResource("Rivet");
0030
0031   _hepmcCollection = consumes<HepMCProduct>(pset.getParameter<edm::InputTag>("HepMCCollection"));
0032   _genLumiInfoToken = consumes<GenLumiInfoHeader, edm::InLumi>(pset.getParameter<edm::InputTag>("genLumiInfo"));
0033
0034   _useLHEweights = pset.getParameter<bool>("useLHEweights");
0035   if (_useLHEweights) {
0036     _lheRunInfoToken = consumes<LHERunInfoProduct, edm::InRun>(_lheLabel);
0037     _LHECollection = consumes<LHEEventProduct>(_lheLabel);
0038   }
0039
0040   _weightCap = pset.getParameter<double>("weightCap");
0041   _NLOSmearing = pset.getParameter<double>("NLOSmearing");
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   if (_produceDQM) {
0051     // book stuff needed for DQM
0052     dbe = nullptr;
0053     dbe = edm::Service<DQMStore>().operator->();
0054   }
0055 }
0056
0057 RivetAnalyzer::~RivetAnalyzer() {}
0058
0059 void RivetAnalyzer::beginJob() {
0060   //set the environment, very ugly but rivet is monolithic when it comes to paths
0061   char* cmsswbase = std::getenv("CMSSW_BASE");
0062   char* cmsswrelease = std::getenv("CMSSW_RELEASE_BASE");
0063   if (!std::getenv("RIVET_REF_PATH")) {
0064     const std::string rivetref = string(cmsswbase) +
0065                                  "/src/GeneratorInterface/RivetInterface/data:" + string(cmsswrelease) +
0066                                  "/src/GeneratorInterface/RivetInterface/data:.";
0067     char* rivetrefCstr = strdup(rivetref.c_str());
0068     setenv("RIVET_REF_PATH", rivetrefCstr, 1);
0069     free(rivetrefCstr);
0070   }
0071   if (!std::getenv("RIVET_INFO_PATH")) {
0072     const std::string rivetinfo = string(cmsswbase) +
0073                                   "/src/GeneratorInterface/RivetInterface/data:" + string(cmsswrelease) +
0074                                   "/src/GeneratorInterface/RivetInterface/data:.";
0075     char* rivetinfoCstr = strdup(rivetinfo.c_str());
0076     setenv("RIVET_INFO_PATH", rivetinfoCstr, 1);
0077     free(rivetinfoCstr);
0078   }
0079 }
0080
0081 void RivetAnalyzer::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0082   if (_useLHEweights) {
0083     edm::Handle<LHERunInfoProduct> lheRunInfoHandle;
0084     iRun.getByLabel(_lheLabel, lheRunInfoHandle);
0085     typedef std::vector<LHERunInfoProduct::Header>::const_iterator headers_const_iterator;
0086
0087     std::regex reg("<weight.*> ?(.*?) ?<\\/weight>");
0088     for (headers_const_iterator iter = lheRunInfoHandle->headers_begin(); iter != lheRunInfoHandle->headers_end();
0089          iter++) {
0090       std::vector<std::string> lines = iter->lines();
0091       for (unsigned int iLine = 0; iLine < lines.size(); iLine++) {
0092         std::smatch match;
0093         std::regex_search(lines.at(iLine), match, reg);
0094         if (!match.empty()) {
0095           _lheWeightNames.push_back(match[1]);
0096         }
0097       }
0098     }
0099   }
0100 }
0101
0102 void RivetAnalyzer::beginLuminosityBlock(const edm::LuminosityBlock& iLumi, const edm::EventSetup& iSetup) {
0103   edm::Handle<GenLumiInfoHeader> genLumiInfoHandle;
0104   if (iLumi.getByToken(_genLumiInfoToken, genLumiInfoHandle)) {
0105     _weightNames = genLumiInfoHandle->weightNames();
0106   }
0107
0108   // need to reset the default weight name (or plotting will fail)
0109   if (!_weightNames.empty()) {
0110     _weightNames[0] = "";
0111   } else {  // Summer16 samples have 1 weight stored in HepMC but no weightNames
0112     _weightNames.push_back("");
0113   }
0114 }
0115
0116 void RivetAnalyzer::endLuminosityBlock(const edm::LuminosityBlock& iLumi, const edm::EventSetup& iSetup) { return; }
0117
0118 void RivetAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0119   //finalize weight names on the first event
0120   if (_isFirstEvent) {
0121     if (_useLHEweights) {
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   }
0129
0130   //get the hepmc product from the event
0131   edm::Handle<HepMCProduct> evt;
0132   iEvent.getByToken(_hepmcCollection, evt);
0133
0134   // get HepMC GenEvent
0135   const HepMC::GenEvent* myGenEvent = evt->GetEvent();
0136   std::unique_ptr<HepMC::GenEvent> tmpGenEvtPtr;
0137   //if you want to use an external weight or set the cross section we have to clone the GenEvent and change the weight
0138   tmpGenEvtPtr = std::make_unique<HepMC::GenEvent>(*(evt->GetEvent()));
0139
0140   if (_xsection > 0) {
0141     HepMC::GenCrossSection xsec;
0142     xsec.set_cross_section(_xsection);
0143     tmpGenEvtPtr->set_cross_section(xsec);
0144   }
0145
0146   std::vector<double> mergedWeights;
0147   for (unsigned int i = 0; i < tmpGenEvtPtr->weights().size(); i++) {
0148     mergedWeights.push_back(tmpGenEvtPtr->weights()[i]);
0149   }
0150
0151   if (_useLHEweights) {
0152     edm::Handle<LHEEventProduct> lheEventHandle;
0153     iEvent.getByToken(_LHECollection, lheEventHandle);
0154     for (unsigned int i = 0; i < _lheWeightNames.size(); i++) {
0155       mergedWeights.push_back(tmpGenEvtPtr->weights()[0] * lheEventHandle->weights().at(i).wgt /
0156                               lheEventHandle->originalXWGTUP());
0157     }
0158   }
0159
0160   tmpGenEvtPtr->weights().clear();
0161   for (unsigned int i = 0; i < _cleanedWeightNames.size(); i++) {
0162     tmpGenEvtPtr->weights()[_cleanedWeightNames[i]] = mergedWeights[i];
0163   }
0164   myGenEvent = tmpGenEvtPtr.get();
0165
0166   //apply the beams initialization on the first event
0167   if (_isFirstEvent) {
0168     _analysisHandler = std::make_unique<Rivet::AnalysisHandler>();
0169     _analysisHandler->addAnalyses(_analysisNames);
0170
0171     /// Set analysis handler weight options
0172     _analysisHandler->skipMultiWeights(_skipMultiWeights);
0173     _analysisHandler->selectMultiWeights(_selectMultiWeights);
0174     _analysisHandler->deselectMultiWeights(_deselectMultiWeights);
0175     _analysisHandler->setNominalWeightName(_setNominalWeightName);
0176     _analysisHandler->setWeightCap(_weightCap);
0177     _analysisHandler->setNLOSmearing(_NLOSmearing);
0178
0179     _analysisHandler->init(*myGenEvent);
0180
0181     _isFirstEvent = false;
0182   }
0183
0184   //run the analysis
0185   _analysisHandler->analyze(*myGenEvent);
0186 }
0187
0188 void RivetAnalyzer::endRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0189   if (_doFinalize)
0190     _analysisHandler->finalize();
0191   else {
0192     //if we don't finalize we just want to do the transformation from histograms to DPS
0193     ////normalizeTree(_analysisHandler->tree());
0194     //normalizeTree();
0195   }
0196   _analysisHandler->writeData(_outFileName);
0197
0198   return;
0199 }
0200
0201 //from Rivet 2.X: Analysis.hh (cls 18Feb2014)
0202 /// List of registered analysis data objects
0203 //const vector<AnalysisObjectPtr>& analysisObjects() const {
0204 //return _analysisobjects;
0205 //}
0206
0207 void RivetAnalyzer::endJob() {}
0208
0209 void RivetAnalyzer::normalizeTree() {
0210   using namespace YODA;
0211   std::vector<string> analyses = _analysisHandler->analysisNames();
0212
0213   //tree.ls(".", true);
0214   const string tmpdir = "/RivetNormalizeTmp";
0215   //tree.mkdir(tmpdir);
0216   for (const string& analysis : analyses) {
0217     if (_produceDQM) {
0218       dbe->setCurrentFolder("Rivet/" + analysis);
0219       //global variables that are always present
0220       //sumOfWeights
0221       TH1F nevent("nEvt", "n analyzed Events", 1, 0., 1.);
0222       nevent.SetBinContent(1, _analysisHandler->sumW());
0223       _mes.push_back(dbe->book1D("nEvt", &nevent));
0224     }
0225     //cross section
0226     //TH1F xsection("xSection", "Cross Section", 1, 0., 1.);
0227     //xsection.SetBinContent(1,_analysisHandler->crossSection());
0228     //_mes.push_back(dbe->book1D("xSection",&xsection));
0229     //now loop over the histograms
0230
0231     /*
0232     const vector<string> paths = tree.listObjectNames("/"+analysis, true); // args set recursive listing
0233     std::cout << "Number of objects in YODA tree for analysis " << analysis << " = " << paths.size() << std::endl;
0234     foreach (const string& path, paths) {
0235       IManagedObject* hobj = tree.find(path);
0236       if (hobj) {
0237         // Weird seg fault on SLC4 when trying to dyn cast an IProfile ptr to a IHistogram
0238         // Fix by attempting to cast to IProfile first, only try IHistogram if it fails.
0239         IHistogram1D* histo = 0;
0240         IProfile1D* prof = dynamic_cast<IProfile1D*>(hobj);
0241         if (!prof) histo = dynamic_cast<IHistogram1D*>(hobj);
0242
0243         std::cout << "Converting histo " << path << " to DPS" << std::endl;
0244         tree.mv(path, tmpdir);
0245         const size_t lastslash = path.find_last_of("/");
0246         const string basename = path.substr(lastslash+1, path.length() - (lastslash+1));
0247         const string tmppath = tmpdir + "/" + basename;
0248
0249         // If it's a normal histo:
0250         if (histo) {
0251           IHistogram1D* tmphisto = dynamic_cast<IHistogram1D*>(tree.find(tmppath));
0252           if (tmphisto) {
0253             _analysisHandler->datapointsetFactory().create(path, *tmphisto);
0254           }
0255           //now convert to root and then ME
0256     //need aida2flat (from Rivet 1.X) & flat2root here
0257           TH1F* h = aida2root<IHistogram1D, TH1F>(histo, basename);
0258           if (_produceDQM)
0259             _mes.push_back(dbe->book1D(h->GetName(), h));
0260           delete h;
0261           tree.rm(tmppath);
0262         }
0263         // If it's a profile histo:
0264         else if (prof) {
0265           IProfile1D* tmpprof = dynamic_cast<IProfile1D*>(tree.find(tmppath));
0266           if (tmpprof) {
0267             _analysisHandler->datapointsetFactory().create(path, *tmpprof);
0268           }
0269           //now convert to root and then ME
0270     //need aida2flat (from Rivet 1.X) & flat2root here
0271           TProfile* p = aida2root<IProfile1D, TProfile>(prof, basename);
0272           if (_produceDQM)
0273             _mes.push_back(dbe->bookProfile(p->GetName(), p));
0274           delete p;
0275           tree.rm(tmppath);
0276         }
0277       }
0278     }
0279     */
0280   }
0281   //tree.rmdir(tmpdir);
0282 }
0283
0284 DEFINE_FWK_MODULE(RivetAnalyzer);