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