Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-06-23 01:38:42 

0001 /* \class HTXSRivetProducer
0002  *
0003  * \author David Sperka, University of Florida
0004  *
0005  * $Id: HTXSRivetProducer.cc,v 1.1 2016/09/27 13:07:29 dsperka Exp $
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   //get the hepmc product from the event
0065   edm::Handle<HepMCProduct> evt;
0066 
0067   bool product_exists = iEvent.getByToken(_hepmcCollection, evt);
0068   if (product_exists) {
0069     // get HepMC GenEvent
0070     const HepMC::GenEvent* myGenEvent = evt->GetEvent();
0071 
0072     if (_prodMode == "AUTO") {
0073       // for these prod modes, don't change what is set in BeginRun
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       // set the production mode if not done already
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       // at this point the production mode must be known
0148       if (m_HiggsProdMode == HTXS::UNKNOWN) {
0149         edm::LogInfo("HTXSRivetProducer") << "HTXSRivetProducer WARNING: HiggsProduction mode is UNKNOWN" << endl;
0150       }
0151 
0152       // initialize rivet analysis
0153       _analysisHandler->init(*myGenEvent);
0154     }
0155 
0156     // classify the event
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           // POWHEG
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           // MC@NLO
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);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!