File indexing completed on 2022-05-04 02:52: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
0024
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
0048 _xsection = pset.getParameter<double>("CrossSection");
0049
0050 if (_produceDQM) {
0051
0052 dbe = nullptr;
0053 dbe = edm::Service<DQMStore>().operator->();
0054 }
0055 }
0056
0057 RivetAnalyzer::~RivetAnalyzer() {}
0058
0059 void RivetAnalyzer::beginJob() {
0060
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
0109 if (!_weightNames.empty()) {
0110 _weightNames[0] = "";
0111 } else {
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
0120 if (_isFirstEvent) {
0121 if (_useLHEweights) {
0122 _weightNames.insert(_weightNames.end(), _lheWeightNames.begin(), _lheWeightNames.end());
0123 }
0124
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
0131 edm::Handle<HepMCProduct> evt;
0132 iEvent.getByToken(_hepmcCollection, evt);
0133
0134
0135 const HepMC::GenEvent* myGenEvent = evt->GetEvent();
0136 std::unique_ptr<HepMC::GenEvent> tmpGenEvtPtr;
0137
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
0167 if (_isFirstEvent) {
0168 _analysisHandler = std::make_unique<Rivet::AnalysisHandler>();
0169 _analysisHandler->addAnalyses(_analysisNames);
0170
0171
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
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
0193
0194
0195 }
0196 _analysisHandler->writeData(_outFileName);
0197
0198 return;
0199 }
0200
0201
0202
0203
0204
0205
0206
0207 void RivetAnalyzer::endJob() {}
0208
0209 void RivetAnalyzer::normalizeTree() {
0210 using namespace YODA;
0211 std::vector<string> analyses = _analysisHandler->analysisNames();
0212
0213
0214 const string tmpdir = "/RivetNormalizeTmp";
0215
0216 for (const string& analysis : analyses) {
0217 if (_produceDQM) {
0218 dbe->setCurrentFolder("Rivet/" + analysis);
0219
0220
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
0226
0227
0228
0229
0230
0231
0232
0233
0234
0235
0236
0237
0238
0239
0240
0241
0242
0243
0244
0245
0246
0247
0248
0249
0250
0251
0252
0253
0254
0255
0256
0257
0258
0259
0260
0261
0262
0263
0264
0265
0266
0267
0268
0269
0270
0271
0272
0273
0274
0275
0276
0277
0278
0279
0280 }
0281
0282 }
0283
0284 DEFINE_FWK_MODULE(RivetAnalyzer);