METPUSubtraction/interface/PFMEtSignInterfaceBase.h

0001 #ifndef RecoMET_METPUSubtraction_PFMEtSignInterfaceBase_h
0002 #define RecoMET_METPUSubtraction_PFMEtSignInterfaceBase_h
0003
0004 /** \class PFMEtSignInterfaceBase
0005  *
0006  * Auxiliary class interfacing the TauAnalysis software to
0007  *  RecoMET/METAlgorithms/interface/significanceAlgo.h
0008  * for computing (PF)MEt significance
0009  * (see CMS AN-10/400 for description of the (PF)MEt significance computation)
0010  *
0011  * \author Christian Veelken, LLR
0012  *
0013  */
0014
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "FWCore/Utilities/interface/InputTag.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018
0019 #include "RecoMET/METAlgorithms/interface/SignAlgoResolutions.h"
0020
0021 #include "DataFormats/Candidate/interface/Candidate.h"
0022 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0023 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0024 #include "DataFormats/JetReco/interface/PFJet.h"
0025 #include "DataFormats/METReco/interface/MET.h"
0026 #include "DataFormats/METReco/interface/SigInputObj.h"
0027 #include "DataFormats/MuonReco/interface/Muon.h"
0028 #include "DataFormats/PatCandidates/interface/Electron.h"
0029 #include "DataFormats/PatCandidates/interface/Photon.h"
0030 #include "DataFormats/PatCandidates/interface/Muon.h"
0031 #include "DataFormats/PatCandidates/interface/Tau.h"
0032 #include "DataFormats/PatCandidates/interface/Jet.h"
0033 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0034 #include "DataFormats/TauReco/interface/PFTau.h"
0035 #include "DataFormats/TrackReco/interface/Track.h"
0036
0037 #include <TFile.h>
0038 #include <TH2.h>
0039
0040 class PFMEtSignInterfaceBase {
0041 public:
0042   PFMEtSignInterfaceBase(const edm::ParameterSet&);
0043   ~PFMEtSignInterfaceBase();
0044
0045   template <typename T>
0046   metsig::SigInputObj compResolution(const T* particle) const {
0047     double pt = particle->pt();
0048     double eta = particle->eta();
0049     double phi = particle->phi();
0050
0051     if (dynamic_cast<const reco::GsfElectron*>(particle) != nullptr ||
0052         dynamic_cast<const pat::Electron*>(particle) != nullptr) {
0053       std::string particleType = "electron";
0054       // WARNING: SignAlgoResolutions::PFtype2 needs to be kept in sync with reco::PFCandidate::e !!
0055       double dpt = pfMEtResolution_->eval(metsig::PFtype2, metsig::ET, pt, phi, eta);
0056       double dphi = pfMEtResolution_->eval(metsig::PFtype2, metsig::PHI, pt, phi, eta);
0057       //std::cout << "electron: pt = " << pt << ", eta = " << eta << ", phi = " << phi
0058       //          << " --> dpt = " << dpt << ", dphi = " << dphi << std::endl;
0059       return metsig::SigInputObj(particleType, pt, phi, dpt, dphi);
0060     } else if (dynamic_cast<const reco::Photon*>(particle) != nullptr ||
0061                dynamic_cast<const pat::Photon*>(particle) != nullptr) {
0062       // CV: assume resolutions for photons to be identical to electron resolutions
0063       std::string particleType = "electron";
0064       // WARNING: SignAlgoResolutions::PFtype2 needs to be kept in sync with reco::PFCandidate::e !!
0065       double dpt = pfMEtResolution_->eval(metsig::PFtype2, metsig::ET, pt, phi, eta);
0066       double dphi = pfMEtResolution_->eval(metsig::PFtype2, metsig::PHI, pt, phi, eta);
0067       //std::cout << "photon: pt = " << pt << ", eta = " << eta << ", phi = " << phi
0068       //          << " --> dpt = " << dpt << ", dphi = " << dphi << std::endl;
0069       return metsig::SigInputObj(particleType, pt, phi, dpt, dphi);
0070     } else if (dynamic_cast<const reco::Muon*>(particle) != nullptr ||
0071                dynamic_cast<const pat::Muon*>(particle) != nullptr) {
0072       std::string particleType = "muon";
0073       double dpt, dphi;
0074       const reco::Track* muonTrack = nullptr;
0075       if (dynamic_cast<const pat::Muon*>(particle) != nullptr) {
0076         const pat::Muon* muon = dynamic_cast<const pat::Muon*>(particle);
0077         if (muon->track().isNonnull() && muon->track().isAvailable())
0078           muonTrack = muon->track().get();
0079       } else if (dynamic_cast<const reco::Muon*>(particle) != nullptr) {
0080         const reco::Muon* muon = dynamic_cast<const reco::Muon*>(particle);
0081         if (muon->track().isNonnull() && muon->track().isAvailable())
0082           muonTrack = muon->track().get();
0083       } else
0084         assert(0);
0085       if (muonTrack) {
0086         dpt = muonTrack->ptError();
0087         dphi = pt * muonTrack->phiError();  // CV: pt*dphi is indeed correct
0088       } else {
0089         // WARNING: SignAlgoResolutions::PFtype3 needs to be kept in sync with reco::PFCandidate::mu !!
0090         dpt = pfMEtResolution_->eval(metsig::PFtype3, metsig::ET, pt, phi, eta);
0091         dphi = pfMEtResolution_->eval(metsig::PFtype3, metsig::PHI, pt, phi, eta);
0092       }
0093       //std::cout << "muon: pt = " << pt << ", eta = " << eta << ", phi = " << phi
0094       //      << " --> dpt = " << dpt << ", dphi = " << dphi << std::endl;
0095       return metsig::SigInputObj(particleType, pt, phi, dpt, dphi);
0096     } else if (dynamic_cast<const reco::PFTau*>(particle) != nullptr ||
0097                dynamic_cast<const pat::Tau*>(particle) != nullptr) {
0098       // CV: use PFJet resolutions for PFTaus for now...
0099       //    (until PFTau specific resolutions are available)
0100       if (dynamic_cast<const pat::Tau*>(particle) != nullptr) {
0101         const pat::Tau* pfTau = dynamic_cast<const pat::Tau*>(particle);
0102         return pfMEtResolution_->evalPFJet(pfTau->pfJetRef().get());
0103       } else if (dynamic_cast<const reco::PFTau*>(particle) != nullptr) {
0104         const reco::PFTau* pfTau = dynamic_cast<const reco::PFTau*>(particle);
0105         return pfMEtResolution_->evalPFJet(pfTau->jetRef().get());
0106       } else
0107         assert(0);
0108     } else if (dynamic_cast<const reco::PFJet*>(particle) != nullptr ||
0109                dynamic_cast<const pat::Jet*>(particle) != nullptr) {
0110       metsig::SigInputObj pfJetResolution;
0111       if (dynamic_cast<const reco::PFJet*>(particle) != nullptr) {
0112         const reco::PFJet* pfJet = dynamic_cast<const reco::PFJet*>(particle);
0113         pfJetResolution = pfMEtResolution_->evalPFJet(pfJet);
0114       } else if (dynamic_cast<const pat::Jet*>(particle) != nullptr) {
0115         const pat::Jet* jet = dynamic_cast<const pat::Jet*>(particle);
0116         if (jet->isPFJet()) {
0117           reco::PFJet pfJet(jet->p4(), jet->vertex(), jet->pfSpecific(), jet->getJetConstituents());
0118           pfJetResolution = pfMEtResolution_->evalPFJet(&pfJet);
0119         } else
0120           throw cms::Exception("addPFMEtSignObjects") << "PAT jet not of PF-type !!\n";
0121       } else
0122         assert(0);
0123       //std::cout << "pfJet: pt = " << pt << ", eta = " << eta << ", phi = " << phi << std::endl;
0124       // CV: apply additional jet energy resolution corrections
0125       //     not included in (PF)MEt significance algorithm yet
0126       //    (cf. CMS AN-11/400 vs. CMS AN-11/330)
0127       if (lut_ && pt > 10.) {
0128         double x = std::abs(eta);
0129         double y = pt;
0130         if (x > lut_->GetXaxis()->GetXmin() && x < lut_->GetXaxis()->GetXmax() && y > lut_->GetYaxis()->GetXmin() &&
0131             y < lut_->GetYaxis()->GetXmax()) {
0132           int binIndex = lut_->FindBin(x, y);
0133           double addJERcorrFactor = lut_->GetBinContent(binIndex);
0134           //std::cout << " addJERcorrFactor = " << addJERcorrFactor << std::endl;
0135           pfJetResolution.set(pfJetResolution.get_type(),
0136                               pfJetResolution.get_energy(),
0137                               pfJetResolution.get_phi(),
0138                               addJERcorrFactor * pfJetResolution.get_sigma_e(),
0139                               pfJetResolution.get_sigma_tan());
0140         }
0141       }
0142       return pfJetResolution;
0143     } else if (dynamic_cast<const reco::PFCandidate*>(particle) != nullptr) {
0144       const reco::PFCandidate* pfCandidate = dynamic_cast<const reco::PFCandidate*>(particle);
0145       //std::cout << "pfCandidate: pt = " << pt << ", eta = " << eta << ", phi = " << phi << std::endl;
0146       return pfMEtResolution_->evalPF(pfCandidate);
0147     } else
0148       throw cms::Exception("addPFMEtSignObjects")
0149           << "Invalid type of particle:"
0150           << " valid types = { reco::GsfElectron/pat::Electron, reco::Photon/pat::Photon, reco::Muon/pat::Muon, "
0151              "reco::PFTau/pat::Tau,"
0152           << " reco::PFJet/pat::Jet, reco::PFCandidate } !!\n";
0153   }
0154
0155   template <typename T>
0156   std::vector<metsig::SigInputObj> compResolution(const std::list<T*>& particles) const {
0157     LogDebug("compResolution") << " particles: entries = " << particles.size() << std::endl;
0158
0159     std::vector<metsig::SigInputObj> pfMEtSignObjects;
0160     addPFMEtSignObjects(pfMEtSignObjects, particles);
0161
0162     return pfMEtSignObjects;
0163   }
0164
0165   reco::METCovMatrix operator()(const std::vector<metsig::SigInputObj>&) const;
0166
0167   template <typename T>
0168   reco::METCovMatrix operator()(const std::list<T*>& particles) const {
0169     std::vector<metsig::SigInputObj> pfMEtSignObjects = compResolution(particles);
0170     return this->operator()(pfMEtSignObjects);
0171   }
0172
0173 protected:
0174   template <typename T>
0175   void addPFMEtSignObjects(std::vector<metsig::SigInputObj>& metSignObjects, const std::list<T*>& particles) const {
0176     for (typename std::list<T*>::const_iterator particle = particles.begin(); particle != particles.end(); ++particle) {
0177       metSignObjects.push_back(this->compResolution(*particle));
0178     }
0179   }
0180
0181 private:
0182   metsig::SignAlgoResolutions* pfMEtResolution_;
0183
0184   // CV: look-up table for additional jet energy resolution corrections
0185   //     not included in (PF)MEt significance algorithm yet
0186   //    (cf. CMS AN-11/400 vs. CMS AN-11/330)
0187   TFile* inputFile_;
0188   TH2* lut_;
0189
0190   int verbosity_;
0191 };
0192
0193 #endif