PFMEtSignInterfaceBase Class Reference

#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 41 of file PFMEtSignInterfaceBase.h.

Constructor & Destructor Documentation

PFMEtSignInterfaceBase::PFMEtSignInterfaceBase

(

const edm::ParameterSet &

cfg

)

Definition at line 15 of file PFMEtSignInterfaceBase.cc.

References edm::hlt::Exception, edm::ParameterSet::exists(), edm::ParameterSet::getParameter(), inputFile_, PFRecoTauDiscriminationAgainstElectronMVA5_cfi::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 41 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 178 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 49 of file PFMEtSignInterfaceBase.h.

References funct::abs(), metsig::ET, eta(), metsig::SignAlgoResolutions::eval(), metsig::SignAlgoResolutions::evalPF(), metsig::SignAlgoResolutions::evalPFJet(), edm::hlt::Exception, 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_, trackingTruthProducerFastSim_cfi::muon, reco::LeafCandidate::p4(), objects.autophobj::particleType, pfJet, pat::Tau::pfJetRef(), pfMEtResolution_, pat::Jet::pfSpecific(), metsig::PFtype2, metsig::PFtype3, phi, metsig::PHI, reco::TrackBase::phiError(), EnergyCorrector::pt, reco::TrackBase::ptError(), metsig::SigInputObj::set(), AlCaHLTBitMon_QueryRunRegistry::string, reco::Muon::track(), pat::Muon::track(), reco::LeafCandidate::vertex(), x, and detailsBasic3DVector::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) != 0 ||
     dynamic_cast<const pat::Electron*>(particle) != 0 ) {
      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) != 0 ||
        dynamic_cast<const pat::Photon*>(particle) != 0 ) {
      // 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) != 0 ||
        dynamic_cast<const pat::Muon*>(particle) != 0 ) {
      std::string particleType = "muon";
      double dpt, dphi;
      const reco::Track* muonTrack = nullptr;
    if ( dynamic_cast<const pat::Muon*>(particle) != 0 ) {
      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) != 0 ) {
      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) != 0 ||
        dynamic_cast<const pat::Tau*>(particle) != 0 ) {
      // CV: use PFJet resolutions for PFTaus for now...
      //    (until PFTau specific resolutions are available)
      if ( dynamic_cast<const pat::Tau*>(particle) != 0 ) {
    const pat::Tau* pfTau = dynamic_cast<const pat::Tau*>(particle);
    //std::cout << "tau: pt = " << pt << ", eta = " << eta << ", phi = " << phi << std::endl;
    return pfMEtResolution_->evalPFJet(pfTau->pfJetRef().get());
      } else if ( dynamic_cast<const reco::PFTau*>(particle) != 0  ) {
    const reco::PFTau* pfTau = dynamic_cast<const reco::PFTau*>(particle);
    //std::cout << "tau: pt = " << pt << ", eta = " << eta << ", phi = " << phi << std::endl;
    return pfMEtResolution_->evalPFJet(pfTau->jetRef().get());
      } else assert(0);
    } else if ( dynamic_cast<const reco::PFJet*>(particle) != 0 ||
        dynamic_cast<const pat::Jet*>(particle) != 0 ) {
      metsig::SigInputObj pfJetResolution;
      if ( dynamic_cast<const reco::PFJet*>(particle) != 0 ) {
    const reco::PFJet* pfJet = dynamic_cast<const reco::PFJet*>(particle);
    pfJetResolution = pfMEtResolution_->evalPFJet(pfJet);
      } else if ( dynamic_cast<const pat::Jet*>(particle) != 0 ) {
    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) != 0 ) {
      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(), and LogDebug.

  {
    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 48 of file PFMEtSignInterfaceBase.cc.

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

Referenced by 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 169 of file PFMEtSignInterfaceBase.h.

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

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

Member Data Documentation

TFile* PFMEtSignInterfaceBase::inputFile_

private

Definition at line 194 of file PFMEtSignInterfaceBase.h.

Referenced by PFMEtSignInterfaceBase(), and ~PFMEtSignInterfaceBase().

TH2* PFMEtSignInterfaceBase::lut_

private

Definition at line 195 of file PFMEtSignInterfaceBase.h.

Referenced by compResolution(), PFMEtSignInterfaceBase(), and ~PFMEtSignInterfaceBase().

metsig::SignAlgoResolutions* PFMEtSignInterfaceBase::pfMEtResolution_

private

Definition at line 189 of file PFMEtSignInterfaceBase.h.

Referenced by compResolution(), PFMEtSignInterfaceBase(), and ~PFMEtSignInterfaceBase().

int PFMEtSignInterfaceBase::verbosity_

private

Definition at line 197 of file PFMEtSignInterfaceBase.h.

Referenced by PFMEtSignInterfaceBase().