METProducers/src/METSignificanceProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    METProducers
0004 // Class:      METSignificanceProducer
0005 //
0006 //
0007
0008 //____________________________________________________________________________||
0009 #include "RecoMET/METProducers/interface/METSignificanceProducer.h"
0010
0011 //____________________________________________________________________________||
0012 namespace cms {
0013
0014   //____________________________________________________________________________||
0015   METSignificanceProducer::METSignificanceProducer(const edm::ParameterSet& iConfig) {
0016     std::vector<edm::InputTag> srcLeptonsTags = iConfig.getParameter<std::vector<edm::InputTag>>("srcLeptons");
0017     for (std::vector<edm::InputTag>::const_iterator it = srcLeptonsTags.begin(); it != srcLeptonsTags.end(); it++) {
0018       lepTokens_.push_back(consumes<edm::View<reco::Candidate>>(*it));
0019     }
0020
0021     pfjetsToken_ = consumes<edm::View<reco::Jet>>(iConfig.getParameter<edm::InputTag>("srcPfJets"));
0022
0023     metToken_ = consumes<edm::View<reco::MET>>(iConfig.getParameter<edm::InputTag>("srcMet"));
0024     pfCandidatesToken_ = consumes<edm::View<reco::Candidate>>(iConfig.getParameter<edm::InputTag>("srcPFCandidates"));
0025
0026     jetSFType_ = iConfig.getParameter<std::string>("srcJetSF");
0027     jetResPtType_ = iConfig.getParameter<std::string>("srcJetResPt");
0028     jetResPhiType_ = iConfig.getParameter<std::string>("srcJetResPhi");
0029     jetSFTypeToken_ = esConsumes(edm::ESInputTag("", jetSFType_));
0030     jetResPtTypeToken_ = esConsumes(edm::ESInputTag("", jetResPtType_));
0031     jetResPhiTypeToken_ = esConsumes(edm::ESInputTag("", jetResPhiType_));
0032     rhoToken_ = consumes<double>(iConfig.getParameter<edm::InputTag>("srcRho"));
0033
0034     edm::InputTag srcWeights = iConfig.getParameter<edm::InputTag>("srcWeights");
0035     if (!srcWeights.label().empty())
0036       weightsToken_ = consumes<edm::ValueMap<float>>(srcWeights);
0037
0038     metSigAlgo_ = new metsig::METSignificance(iConfig);
0039
0040     produces<double>("METSignificance");
0041     produces<math::Error<2>::type>("METCovariance");
0042   }
0043
0044   METSignificanceProducer::~METSignificanceProducer() { delete metSigAlgo_; }
0045
0046   //____________________________________________________________________________||
0047   void METSignificanceProducer::produce(edm::Event& event, const edm::EventSetup& setup) {
0048     //
0049     // met
0050     //
0051     edm::Handle<edm::View<reco::MET>> metHandle;
0052     event.getByToken(metToken_, metHandle);
0053     const reco::MET& met = (*metHandle)[0];
0054
0055     //
0056     // candidates
0057     //
0058     edm::Handle<reco::CandidateView> pfCandidates;
0059     event.getByToken(pfCandidatesToken_, pfCandidates);
0060
0061     //
0062     // leptons
0063     //
0064     std::vector<edm::Handle<reco::CandidateView>> leptons;
0065     for (std::vector<edm::EDGetTokenT<edm::View<reco::Candidate>>>::const_iterator srcLeptons_i = lepTokens_.begin();
0066          srcLeptons_i != lepTokens_.end();
0067          ++srcLeptons_i) {
0068       edm::Handle<reco::CandidateView> leptons_i;
0069       event.getByToken(*srcLeptons_i, leptons_i);
0070       leptons.push_back(leptons_i);
0071     }
0072
0073     //
0074     // jets
0075     //
0076     edm::Handle<edm::View<reco::Jet>> jets;
0077     event.getByToken(pfjetsToken_, jets);
0078
0079     edm::Handle<double> rho;
0080     event.getByToken(rhoToken_, rho);
0081
0082     edm::Handle<edm::ValueMap<float>> weights;
0083     if (!weightsToken_.isUninitialized())
0084       event.getByToken(weightsToken_, weights);
0085
0086     JME::JetResolution resPtObj = JME::JetResolution::get(setup, jetResPtTypeToken_);
0087     JME::JetResolution resPhiObj = JME::JetResolution::get(setup, jetResPhiTypeToken_);
0088     JME::JetResolutionScaleFactor resSFObj = JME::JetResolutionScaleFactor::get(setup, jetSFTypeToken_);
0089
0090     //
0091     // compute the significance
0092     //
0093     double sumPtUnclustered = 0;
0094     const edm::ValueMap<float>* weightsPtr = nullptr;
0095     if (met.isWeighted()) {
0096       if (weightsToken_.isUninitialized())
0097         throw cms::Exception("InvalidInput")
0098             << "MET is weighted (e.g. PUPPI), but no weights given in METSignificanceProducer\n";
0099       weightsPtr = &*weights;
0100     }
0101     reco::METCovMatrix cov = metSigAlgo_->getCovariance(*jets,
0102                                                         leptons,
0103                                                         pfCandidates,
0104                                                         *rho,
0105                                                         resPtObj,
0106                                                         resPhiObj,
0107                                                         resSFObj,
0108                                                         event.isRealData(),
0109                                                         sumPtUnclustered,
0110                                                         weightsPtr);
0111     double sig = metSigAlgo_->getSignificance(cov, met);
0112
0113     auto significance = std::make_unique<double>();
0114     (*significance) = sig;
0115
0116     auto covPtr = std::make_unique<math::Error<2>::type>();
0117     (*covPtr)(0, 0) = cov(0, 0);
0118     (*covPtr)(1, 0) = cov(1, 0);
0119     (*covPtr)(1, 1) = cov(1, 1);
0120
0121     event.put(std::move(covPtr), "METCovariance");
0122     event.put(std::move(significance), "METSignificance");
0123   }
0124
0125   //____________________________________________________________________________||
0126   DEFINE_FWK_MODULE(METSignificanceProducer);
0127 }  // namespace cms
0128
0129 //____________________________________________________________________________||