File indexing completed on 2023-06-23 01:38:42
0001
0002
0003
0004
0005
0006
0007
0008
0009 #include "FWCore/Framework/interface/one/EDProducer.h"
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/Framework/interface/Run.h"
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014
0015 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0016 #include "SimDataFormats/GeneratorProducts/interface/LHERunInfoProduct.h"
0017
0018 #include "Rivet/AnalysisHandler.hh"
0019 #include "GeneratorInterface/RivetInterface/src/HiggsTemplateCrossSections.cc"
0020 #include "SimDataFormats/HTXS/interface/HiggsTemplateCrossSections.h"
0021
0022 #include <memory>
0023
0024 #include <vector>
0025 #include <cstdio>
0026 #include <cstring>
0027
0028 using namespace Rivet;
0029 using namespace edm;
0030 using namespace std;
0031
0032 class HTXSRivetProducer : public edm::one::EDProducer<edm::one::WatchRuns, edm::one::SharedResources> {
0033 public:
0034 explicit HTXSRivetProducer(const edm::ParameterSet& cfg)
0035 : _hepmcCollection(consumes<HepMCProduct>(cfg.getParameter<edm::InputTag>("HepMCCollection"))),
0036 _lheRunInfo(consumes<LHERunInfoProduct, edm::InRun>(cfg.getParameter<edm::InputTag>("LHERunInfo"))) {
0037 usesResource("Rivet");
0038 _prodMode = cfg.getParameter<string>("ProductionMode");
0039 m_HiggsProdMode = HTXS::UNKNOWN;
0040 _HTXS = nullptr;
0041 _analysisHandler = nullptr;
0042 produces<HTXS::HiggsClassification>("HiggsClassification").setBranchAlias("HiggsClassification");
0043 }
0044
0045 private:
0046 void produce(edm::Event&, const edm::EventSetup&) override;
0047
0048 void beginRun(edm::Run const& iRun, edm::EventSetup const& es) override;
0049 void endRun(edm::Run const& iRun, edm::EventSetup const& es) override;
0050
0051 edm::EDGetTokenT<edm::HepMCProduct> _hepmcCollection;
0052 edm::EDGetTokenT<LHERunInfoProduct> _lheRunInfo;
0053
0054 std::unique_ptr<Rivet::AnalysisHandler> _analysisHandler;
0055 Rivet::HiggsTemplateCrossSections* _HTXS;
0056
0057 std::string _prodMode;
0058 HTXS::HiggsProdMode m_HiggsProdMode;
0059
0060 HTXS::HiggsClassification cat_;
0061 };
0062
0063 void HTXSRivetProducer::produce(edm::Event& iEvent, const edm::EventSetup&) {
0064
0065 edm::Handle<HepMCProduct> evt;
0066
0067 bool product_exists = iEvent.getByToken(_hepmcCollection, evt);
0068 if (product_exists) {
0069
0070 const HepMC::GenEvent* myGenEvent = evt->GetEvent();
0071
0072 if (_prodMode == "AUTO") {
0073
0074 if (m_HiggsProdMode != HTXS::GGF && m_HiggsProdMode != HTXS::VBF && m_HiggsProdMode != HTXS::GG2ZH) {
0075 unsigned nWs = 0;
0076 unsigned nZs = 0;
0077 unsigned nTs = 0;
0078 unsigned nBs = 0;
0079 unsigned nHs = 0;
0080
0081 HepMC::GenVertex* HSvtx = myGenEvent->signal_process_vertex();
0082
0083 if (HSvtx) {
0084 for (auto ptcl : HepMCUtils::particles(HSvtx, HepMC::children)) {
0085 if (std::abs(ptcl->pdg_id()) == 24)
0086 ++nWs;
0087 if (ptcl->pdg_id() == 23)
0088 ++nZs;
0089 if (abs(ptcl->pdg_id()) == 6)
0090 ++nTs;
0091 if (abs(ptcl->pdg_id()) == 5)
0092 ++nBs;
0093 if (ptcl->pdg_id() == 25)
0094 ++nHs;
0095 }
0096 }
0097
0098 if (nZs == 1 && nHs == 1 && (nWs + nTs) == 0) {
0099 m_HiggsProdMode = HTXS::QQ2ZH;
0100 } else if (nWs == 1 && nHs == 1 && (nZs + nTs) == 0) {
0101 m_HiggsProdMode = HTXS::WH;
0102 } else if (nTs == 2 && nHs == 1 && nZs == 0) {
0103 m_HiggsProdMode = HTXS::TTH;
0104 } else if (nTs == 1 && nHs == 1 && nZs == 0) {
0105 m_HiggsProdMode = HTXS::TH;
0106 } else if (nBs == 2 && nHs == 1 && nZs == 0) {
0107 m_HiggsProdMode = HTXS::BBH;
0108 }
0109 }
0110 }
0111
0112 if (!_HTXS || !_HTXS->hasProjection("FS")) {
0113 _analysisHandler = std::make_unique<Rivet::AnalysisHandler>();
0114 _HTXS = new Rivet::HiggsTemplateCrossSections();
0115 _analysisHandler->addAnalysis(_HTXS);
0116
0117
0118 if (_prodMode == "GGF")
0119 m_HiggsProdMode = HTXS::GGF;
0120 else if (_prodMode == "VBF")
0121 m_HiggsProdMode = HTXS::VBF;
0122 else if (_prodMode == "WH")
0123 m_HiggsProdMode = HTXS::WH;
0124 else if (_prodMode == "ZH")
0125 m_HiggsProdMode = HTXS::QQ2ZH;
0126 else if (_prodMode == "QQ2ZH")
0127 m_HiggsProdMode = HTXS::QQ2ZH;
0128 else if (_prodMode == "GG2ZH")
0129 m_HiggsProdMode = HTXS::GG2ZH;
0130 else if (_prodMode == "TTH")
0131 m_HiggsProdMode = HTXS::TTH;
0132 else if (_prodMode == "BBH")
0133 m_HiggsProdMode = HTXS::BBH;
0134 else if (_prodMode == "TH")
0135 m_HiggsProdMode = HTXS::TH;
0136 else if (_prodMode == "AUTO") {
0137 edm::LogInfo("HTXSRivetProducer")
0138 << "Using AUTO for HiggsProdMode, found it to be: " << m_HiggsProdMode << "\n";
0139 edm::LogInfo("HTXSRivetProducer")
0140 << "(UNKNOWN=0, GGF=1, VBF=2, WH=3, QQ2ZH=4, GG2ZH=5, TTH=6, BBH=7, TH=8)" << endl;
0141 } else {
0142 throw cms::Exception("HTXSRivetProducer")
0143 << "ProductionMode must be one of: GGF,VBF,WH,ZH,QQ2ZH,GG2ZH,TTH,BBH,TH,AUTO ";
0144 }
0145 _HTXS->setHiggsProdMode(m_HiggsProdMode);
0146
0147
0148 if (m_HiggsProdMode == HTXS::UNKNOWN) {
0149 edm::LogInfo("HTXSRivetProducer") << "HTXSRivetProducer WARNING: HiggsProduction mode is UNKNOWN" << endl;
0150 }
0151
0152
0153 _analysisHandler->init(*myGenEvent);
0154 }
0155
0156
0157 Rivet::HiggsClassification rivet_cat = _HTXS->classifyEvent(*myGenEvent, m_HiggsProdMode);
0158 cat_ = HTXS::Rivet2Root(rivet_cat);
0159
0160 unique_ptr<HTXS::HiggsClassification> cat(new HTXS::HiggsClassification(cat_));
0161
0162 iEvent.put(std::move(cat), "HiggsClassification");
0163 }
0164 }
0165
0166 void HTXSRivetProducer::endRun(edm::Run const& iRun, edm::EventSetup const& es) {
0167 if (_HTXS)
0168 _HTXS->printClassificationSummary();
0169 }
0170
0171 void HTXSRivetProducer::beginRun(edm::Run const& iRun, edm::EventSetup const& es) {
0172 if (_prodMode == "AUTO") {
0173 edm::Handle<LHERunInfoProduct> run;
0174 bool product_exists = iRun.getByLabel(edm::InputTag("externalLHEProducer"), run);
0175 if (product_exists) {
0176 typedef std::vector<LHERunInfoProduct::Header>::const_iterator headers_const_iterator;
0177 LHERunInfoProduct myLHERunInfoProduct = *(run.product());
0178 for (headers_const_iterator iter = myLHERunInfoProduct.headers_begin(); iter != myLHERunInfoProduct.headers_end();
0179 iter++) {
0180 std::vector<std::string> lines = iter->lines();
0181 for (unsigned int iLine = 0; iLine < lines.size(); iLine++) {
0182 std::string line = lines.at(iLine);
0183
0184 if (line.find("gg_H_quark-mass-effects") != std::string::npos) {
0185 edm::LogInfo("HTXSRivetProducer") << iLine << " " << line << std::endl;
0186 m_HiggsProdMode = HTXS::GGF;
0187 break;
0188 }
0189 if (line.find("Process: HJ") != std::string::npos) {
0190 edm::LogInfo("HTXSRivetProducer") << iLine << " " << line << std::endl;
0191 m_HiggsProdMode = HTXS::GGF;
0192 break;
0193 }
0194 if (line.find("Process: HJJ") != std::string::npos) {
0195 edm::LogInfo("HTXSRivetProducer") << iLine << " " << line << std::endl;
0196 m_HiggsProdMode = HTXS::GGF;
0197 break;
0198 }
0199 if (line.find("VBF_H") != std::string::npos) {
0200 edm::LogInfo("HTXSRivetProducer") << iLine << " " << line << std::endl;
0201 m_HiggsProdMode = HTXS::VBF;
0202 break;
0203 }
0204 if (line.find("HZJ") != std::string::npos) {
0205 edm::LogInfo("HTXSRivetProducer") << iLine << " " << line << std::endl;
0206 m_HiggsProdMode = HTXS::QQ2ZH;
0207 break;
0208 }
0209 if (line.find("ggHZ") != std::string::npos) {
0210 edm::LogInfo("HTXSRivetProducer") << iLine << " " << line << std::endl;
0211 m_HiggsProdMode = HTXS::GG2ZH;
0212 break;
0213 }
0214
0215 if (line.find("ggh012j") != std::string::npos) {
0216 edm::LogInfo("HTXSRivetProducer") << iLine << " " << line << std::endl;
0217 m_HiggsProdMode = HTXS::GGF;
0218 break;
0219 }
0220 if (line.find("vbfh") != std::string::npos) {
0221 edm::LogInfo("HTXSRivetProducer") << iLine << " " << line << std::endl;
0222 m_HiggsProdMode = HTXS::VBF;
0223 break;
0224 }
0225 if (line.find("zh012j") != std::string::npos) {
0226 edm::LogInfo("HTXSRivetProducer") << iLine << " " << line << std::endl;
0227 m_HiggsProdMode = HTXS::QQ2ZH;
0228 break;
0229 }
0230 if (line.find("ggzh01j") != std::string::npos) {
0231 edm::LogInfo("HTXSRivetProducer") << iLine << " " << line << std::endl;
0232 m_HiggsProdMode = HTXS::GG2ZH;
0233 break;
0234 }
0235 }
0236
0237 if (m_HiggsProdMode != HTXS::UNKNOWN)
0238 break;
0239 }
0240 }
0241 }
0242 }
0243
0244 #include "FWCore/Framework/interface/MakerMacros.h"
0245
0246 DEFINE_FWK_MODULE(HTXSRivetProducer);