PFMEtSignInterfaceBase.cc

Go to the documentation of this file.

#include "RecoMET/METPUSubtraction/interface/PFMEtSignInterfaceBase.h"

#include "FWCore/Utilities/interface/Exception.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "RecoMET/METAlgorithms/interface/significanceAlgo.h"

#include <TVectorD.h>

const double defaultPFMEtResolutionX = 10.;
const double defaultPFMEtResolutionY = 10.;

const double epsilon = 1.e-9;

PFMEtSignInterfaceBase::PFMEtSignInterfaceBase(const edm::ParameterSet& cfg)
  : 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()
{
  delete pfMEtResolution_;
  delete inputFile_;
  delete lut_;
}

reco::METCovMatrix PFMEtSignInterfaceBase::operator()(const std::vector<metsig::SigInputObj>& pfMEtSignObjects) const
{
  // 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;
}