#include <PFMEtSignInterfaceBase.h>

Public Member Functions
template<typename T >
metsig::SigInputObj	compResolution (const T *particle) const

template<typename T >
std::vector< metsig::SigInputObj >	compResolution (const std::list< T * > &particles) const

reco::METCovMatrix	operator() (const std::vector< metsig::SigInputObj > &) const

template<typename T >
reco::METCovMatrix	operator() (const std::list< T * > &particles) const

	PFMEtSignInterfaceBase (const edm::ParameterSet &)

	~PFMEtSignInterfaceBase ()

Protected Member Functions
template<typename T >
void	addPFMEtSignObjects (std::vector< metsig::SigInputObj > &metSignObjects, const std::list< T * > &particles) const

Private Attributes
TFile *	inputFile_

TH2 *	lut_

metsig::SignAlgoResolutions *	pfMEtResolution_

int	verbosity_

Detailed Description

Auxiliary class interfacing the TauAnalysis software to RecoMET/METAlgorithms/interface/significanceAlgo.h for computing (PF)MEt significance (see CMS AN-10/400 for description of the (PF)MEt significance computation)

Author: Christian Veelken, LLR

Definition at line 39 of file PFMEtSignInterfaceBase.h.

Constructor & Destructor Documentation

PFMEtSignInterfaceBase::PFMEtSignInterfaceBase ( const edm::ParameterSet & cfg )

Definition at line 15 of file PFMEtSignInterfaceBase.cc.

References Exception, edm::ParameterSet::exists(), edm::ParameterSet::getParameter(), inputFile_, InefficientDoubleROC::inputFileName, edm::FileInPath::Local, lut_, pfMEtResolution_, AlCaHLTBitMon_QueryRunRegistry::string, and verbosity_.

     : pfMEtResolution_(nullptr), inputFile_(nullptr), lut_(nullptr) {
   pfMEtResolution_ = new metsig::SignAlgoResolutions(cfg);
 
   if (cfg.exists("addJERcorr")) {
     edm::ParameterSet cfgJERcorr = cfg.getParameter<edm::ParameterSet>("addJERcorr");
     edm::FileInPath inputFileName = cfgJERcorr.getParameter<edm::FileInPath>("inputFileName");
     std::string lutName = cfgJERcorr.getParameter<std::string>("lutName");
     if (inputFileName.location() != edm::FileInPath::Local)
       throw cms::Exception("PFMEtSignInterfaceBase") << " Failed to find File = " << inputFileName << " !!\n";
 
     inputFile_ = new TFile(inputFileName.fullPath().data());
     lut_ = dynamic_cast<TH2*>(inputFile_->Get(lutName.data()));
     if (!lut_)
       throw cms::Exception("PFMEtSignInterfaceBase") << " Failed to load LUT = " << lutName.data()
                                                      << " from file = " << inputFileName.fullPath().data() << " !!\n";
   }
 
   verbosity_ = cfg.exists("verbosity") ? cfg.getParameter<int>("verbosity") : 0;
 }

PFMEtSignInterfaceBase::~PFMEtSignInterfaceBase ( )

Definition at line 36 of file PFMEtSignInterfaceBase.cc.

References inputFile_, lut_, and pfMEtResolution_.

                                                 {
   delete pfMEtResolution_;
   delete inputFile_;
   delete lut_;
 }

Member Function Documentation

template<typename T >

void PFMEtSignInterfaceBase::addPFMEtSignObjects	(	std::vector< metsig::SigInputObj > &	metSignObjects,
		const std::list< T * > &	particles
	)		const

inlineprotected

Definition at line 174 of file PFMEtSignInterfaceBase.h.

References compResolution().

Referenced by compResolution().

                                                                                                              {
     for (typename std::list<T*>::const_iterator particle = particles.begin(); particle != particles.end(); ++particle) {
       metSignObjects.push_back(this->compResolution(*particle));
     }
   }

template<typename T >

metsig::SigInputObj PFMEtSignInterfaceBase::compResolution ( const T * particle ) const

inline

Definition at line 45 of file PFMEtSignInterfaceBase.h.

References funct::abs(), metsig::ET, PVValHelper::eta, metsig::SignAlgoResolutions::eval(), metsig::SignAlgoResolutions::evalPF(), metsig::SignAlgoResolutions::evalPFJet(), Exception, edm::RefToBase< T >::get(), edm::Ref< C, T, F >::get(), metsig::SigInputObj::get_energy(), metsig::SigInputObj::get_phi(), metsig::SigInputObj::get_sigma_e(), metsig::SigInputObj::get_sigma_tan(), metsig::SigInputObj::get_type(), reco::Jet::getJetConstituents(), edm::Ref< C, T, F >::isAvailable(), edm::Ref< C, T, F >::isNonnull(), pat::Jet::isPFJet(), metsig::jet, reco::PFTau::jetRef(), lut_, reco::LeafCandidate::p4(), PbPb_ZMuSkimMuonDPG_cff::particleType, pat::Tau::pfJetRef(), pfMEtResolution_, pat::Jet::pfSpecific(), metsig::PFtype2, metsig::PFtype3, phi, metsig::PHI, reco::TrackBase::phiError(), DiDispStaMuonMonitor_cfi::pt, reco::TrackBase::ptError(), metsig::SigInputObj::set(), AlCaHLTBitMon_QueryRunRegistry::string, reco::Muon::track(), pat::Muon::track(), reco::LeafCandidate::vertex(), x, and y.

Referenced by addPFMEtSignObjects(), operator()(), NoPileUpPFMEtProducer::produce(), and NoPileUpPFMEtDataProducer::produce().

                                                             {
     double pt = particle->pt();
     double eta = particle->eta();
     double phi = particle->phi();
 
     if (dynamic_cast<const reco::GsfElectron*>(particle) != nullptr ||
         dynamic_cast<const pat::Electron*>(particle) != nullptr) {
       std::string particleType = "electron";
       // WARNING: SignAlgoResolutions::PFtype2 needs to be kept in sync with reco::PFCandidate::e !!
       double dpt = pfMEtResolution_->eval(metsig::PFtype2, metsig::ET, pt, phi, eta);
       double dphi = pfMEtResolution_->eval(metsig::PFtype2, metsig::PHI, pt, phi, eta);
       //std::cout << "electron: pt = " << pt << ", eta = " << eta << ", phi = " << phi
       //          << " --> dpt = " << dpt << ", dphi = " << dphi << std::endl;
       return metsig::SigInputObj(particleType, pt, phi, dpt, dphi);
     } else if (dynamic_cast<const reco::Photon*>(particle) != nullptr ||
                dynamic_cast<const pat::Photon*>(particle) != nullptr) {
       // CV: assume resolutions for photons to be identical to electron resolutions
       std::string particleType = "electron";
       // WARNING: SignAlgoResolutions::PFtype2 needs to be kept in sync with reco::PFCandidate::e !!
       double dpt = pfMEtResolution_->eval(metsig::PFtype2, metsig::ET, pt, phi, eta);
       double dphi = pfMEtResolution_->eval(metsig::PFtype2, metsig::PHI, pt, phi, eta);
       //std::cout << "photon: pt = " << pt << ", eta = " << eta << ", phi = " << phi
       //          << " --> dpt = " << dpt << ", dphi = " << dphi << std::endl;
       return metsig::SigInputObj(particleType, pt, phi, dpt, dphi);
     } else if (dynamic_cast<const reco::Muon*>(particle) != nullptr ||
                dynamic_cast<const pat::Muon*>(particle) != nullptr) {
       std::string particleType = "muon";
       double dpt, dphi;
       const reco::Track* muonTrack = nullptr;
       if (dynamic_cast<const pat::Muon*>(particle) != nullptr) {
         const pat::Muon* muon = dynamic_cast<const pat::Muon*>(particle);
         if (muon->track().isNonnull() && muon->track().isAvailable())
           muonTrack = muon->track().get();
       } else if (dynamic_cast<const reco::Muon*>(particle) != nullptr) {
         const reco::Muon* muon = dynamic_cast<const reco::Muon*>(particle);
         if (muon->track().isNonnull() && muon->track().isAvailable())
           muonTrack = muon->track().get();
       } else
         assert(0);
       if (muonTrack) {
         dpt = muonTrack->ptError();
         dphi = pt * muonTrack->phiError();  // CV: pt*dphi is indeed correct
       } else {
         // WARNING: SignAlgoResolutions::PFtype3 needs to be kept in sync with reco::PFCandidate::mu !!
         dpt = pfMEtResolution_->eval(metsig::PFtype3, metsig::ET, pt, phi, eta);
         dphi = pfMEtResolution_->eval(metsig::PFtype3, metsig::PHI, pt, phi, eta);
       }
       //std::cout << "muon: pt = " << pt << ", eta = " << eta << ", phi = " << phi
       //      << " --> dpt = " << dpt << ", dphi = " << dphi << std::endl;
       return metsig::SigInputObj(particleType, pt, phi, dpt, dphi);
     } else if (dynamic_cast<const reco::PFTau*>(particle) != nullptr ||
                dynamic_cast<const pat::Tau*>(particle) != nullptr) {
       // CV: use PFJet resolutions for PFTaus for now...
       //    (until PFTau specific resolutions are available)
       if (dynamic_cast<const pat::Tau*>(particle) != nullptr) {
         const pat::Tau* pfTau = dynamic_cast<const pat::Tau*>(particle);
         return pfMEtResolution_->evalPFJet(pfTau->pfJetRef().get());
       } else if (dynamic_cast<const reco::PFTau*>(particle) != nullptr) {
         const reco::PFTau* pfTau = dynamic_cast<const reco::PFTau*>(particle);
         return pfMEtResolution_->evalPFJet(pfTau->jetRef().get());
       } else
         assert(0);
     } else if (dynamic_cast<const reco::PFJet*>(particle) != nullptr ||
                dynamic_cast<const pat::Jet*>(particle) != nullptr) {
       metsig::SigInputObj pfJetResolution;
       if (dynamic_cast<const reco::PFJet*>(particle) != nullptr) {
         const reco::PFJet* pfJet = dynamic_cast<const reco::PFJet*>(particle);
         pfJetResolution = pfMEtResolution_->evalPFJet(pfJet);
       } else if (dynamic_cast<const pat::Jet*>(particle) != nullptr) {
         const pat::Jet* jet = dynamic_cast<const pat::Jet*>(particle);
         if (jet->isPFJet()) {
           reco::PFJet pfJet(jet->p4(), jet->vertex(), jet->pfSpecific(), jet->getJetConstituents());
           pfJetResolution = pfMEtResolution_->evalPFJet(&pfJet);
         } else
           throw cms::Exception("addPFMEtSignObjects") << "PAT jet not of PF-type !!\n";
       } else
         assert(0);
       //std::cout << "pfJet: pt = " << pt << ", eta = " << eta << ", phi = " << phi << std::endl;
       // CV: apply additional jet energy resolution corrections
       //     not included in (PF)MEt significance algorithm yet
       //    (cf. CMS AN-11/400 vs. CMS AN-11/330)
       if (lut_ && pt > 10.) {
         double x = std::abs(eta);
         double y = pt;
         if (x > lut_->GetXaxis()->GetXmin() && x < lut_->GetXaxis()->GetXmax() && y > lut_->GetYaxis()->GetXmin() &&
             y < lut_->GetYaxis()->GetXmax()) {
           int binIndex = lut_->FindBin(x, y);
           double addJERcorrFactor = lut_->GetBinContent(binIndex);
           //std::cout << " addJERcorrFactor = " << addJERcorrFactor << std::endl;
           pfJetResolution.set(pfJetResolution.get_type(),
                               pfJetResolution.get_energy(),
                               pfJetResolution.get_phi(),
                               addJERcorrFactor * pfJetResolution.get_sigma_e(),
                               pfJetResolution.get_sigma_tan());
         }
       }
       return pfJetResolution;
     } else if (dynamic_cast<const reco::PFCandidate*>(particle) != nullptr) {
       const reco::PFCandidate* pfCandidate = dynamic_cast<const reco::PFCandidate*>(particle);
       //std::cout << "pfCandidate: pt = " << pt << ", eta = " << eta << ", phi = " << phi << std::endl;
       return pfMEtResolution_->evalPF(pfCandidate);
     } else
       throw cms::Exception("addPFMEtSignObjects")
           << "Invalid type of particle:"
           << " valid types = { reco::GsfElectron/pat::Electron, reco::Photon/pat::Photon, reco::Muon/pat::Muon, "
              "reco::PFTau/pat::Tau,"
           << " reco::PFJet/pat::Jet, reco::PFCandidate } !!\n";
   }

template<typename T >

std::vector<metsig::SigInputObj> PFMEtSignInterfaceBase::compResolution ( const std::list< T * > & particles ) const

inline

Definition at line 155 of file PFMEtSignInterfaceBase.h.

References addPFMEtSignObjects(), LogDebug, and operator()().

                                                                                   {
     LogDebug("compResolution") << " particles: entries = " << particles.size() << std::endl;
 
     std::vector<metsig::SigInputObj> pfMEtSignObjects;
     addPFMEtSignObjects(pfMEtSignObjects, particles);
 
     return pfMEtSignObjects;
   }

reco::METCovMatrix PFMEtSignInterfaceBase::operator() ( const std::vector< metsig::SigInputObj > & pfMEtSignObjects ) const

Definition at line 42 of file PFMEtSignInterfaceBase.cc.

References funct::abs(), metsig::significanceAlgo::addObjects(), defaultPFMEtResolutionX, defaultPFMEtResolutionY, epsilon, and metsig::significanceAlgo::getSignifMatrix().

Referenced by compResolution(), and operator()().

                                                                                                               {
   // if ( this->verbosity_ ) {
   //   std::cout << "<PFMEtSignInterfaceBase::operator()>:" << std::endl;
   //   std::cout << " pfMEtSignObjects: entries = " << pfMEtSignObjects.size() << std::endl;
   //   double dpt2Sum = 0.;
   //   for ( std::vector<metsig::SigInputObj>::const_iterator pfMEtSignObject = pfMEtSignObjects.begin();
   //      pfMEtSignObject != pfMEtSignObjects.end(); ++pfMEtSignObject ) {
   //     std::cout << pfMEtSignObject->get_type() << ": pt = " << pfMEtSignObject->get_energy() << ","
   //        << " phi = " << pfMEtSignObject->get_phi() << " --> dpt = " << pfMEtSignObject->get_sigma_e() << std::endl;
   //     dpt2Sum += pfMEtSignObject->get_sigma_e();
   //   }
   //   std::cout << "--> sqrt(sum(dpt^2)) = " << sqrt(dpt2Sum) << std::endl;
   // }
 
   reco::METCovMatrix pfMEtCov;
   if (pfMEtSignObjects.size() >= 2) {
     metsig::significanceAlgo pfMEtSignAlgorithm;
     pfMEtSignAlgorithm.addObjects(pfMEtSignObjects);
     pfMEtCov = pfMEtSignAlgorithm.getSignifMatrix();
 
     double det = 0;
     pfMEtCov.Det(det);
 
     // if ( this->verbosity_ && std::abs(det) > epsilon ) {
     //   //keep TMatrixD as it is much easier to find
     //   //eigenvectors and values than with SMatrix;
     //   //not used anyway, except for debugging
     //   TMatrixD tmpMatrix(2,2);
     //   tmpMatrix(0,0) = pfMEtCov(0,0);
     //   tmpMatrix(0,1) = pfMEtCov(0,1);
     //   tmpMatrix(1,0) = pfMEtCov(1,0);
     //   tmpMatrix(1,1) = pfMEtCov(1,1);
 
     //   TVectorD eigenValues(2);
     //   TMatrixD eigenVectors = tmpMatrix.EigenVectors(eigenValues);
     //   // CV: eigenvectors are stored in columns
     //   //     and are sorted such that the one corresponding to the highest eigenvalue is in the **first** column
     //   for ( unsigned iEigenVector = 0; iEigenVector < 2; ++iEigenVector ) {
     //  std::cout << "eigenVector #" << iEigenVector << " (eigenValue = " << eigenValues(iEigenVector) << "):"
     //        << " x = " << eigenVectors(0, iEigenVector) << ", y = " << eigenVectors(1, iEigenVector) << std::endl;
     //   }
     // }
 
     //--- substitute (PF)MEt resolution matrix by default values
     //    in case resolution matrix cannot be inverted
     if (std::abs(det) < epsilon) {
       edm::LogWarning("PFMEtSignInterfaceBase::operator()")
           << "Inversion of PFMEt covariance matrix failed, det = " << det
           << " --> replacing covariance matrix by resolution defaults !!";
       pfMEtCov(0, 0) = defaultPFMEtResolutionX * defaultPFMEtResolutionX;
       pfMEtCov(0, 1) = 0.;
       pfMEtCov(1, 0) = 0.;
       pfMEtCov(1, 1) = defaultPFMEtResolutionY * defaultPFMEtResolutionY;
     }
   } else {
     pfMEtCov(0, 0) = 0.;
     pfMEtCov(0, 1) = 0.;
     pfMEtCov(1, 0) = 0.;
     pfMEtCov(1, 1) = 0.;
   }
 
   return pfMEtCov;
 }

template<typename T >

reco::METCovMatrix PFMEtSignInterfaceBase::operator() ( const std::list< T * > & particles ) const

inline

Definition at line 167 of file PFMEtSignInterfaceBase.h.

References compResolution(), and operator()().

                                                                   {
     std::vector<metsig::SigInputObj> pfMEtSignObjects = compResolution(particles);
     return this->operator()(pfMEtSignObjects);
   }