Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-01-09 23:33:35 

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/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   //get the hepmc product from the event
0066   edm::Handle<HepMC3Product> evt;
0067 
0068   bool product_exists = iEvent.getByToken(_hepmcCollection, evt);
0069   if (product_exists) {
0070     // get HepMC GenEvent
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       // for these prod modes, don't change what is set in BeginRun
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       // set the production mode if not done already
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       // at this point the production mode must be known
0151       if (m_HiggsProdMode == HTXS::UNKNOWN) {
0152         edm::LogInfo("HTXSRivetProducer") << "HTXSRivetProducer WARNING: HiggsProduction mode is UNKNOWN" << endl;
0153       }
0154 
0155       // initialize rivet analysis
0156       _analysisHandler->init(*myGenEvent);
0157     }
0158 
0159     // Create the Rivet event wrapper
0160     const Rivet::Event event(const_cast<GenEvent&>(*myGenEvent));
0161     // classify the event
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           // POWHEG
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           // MC@NLO
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);
[ 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!