METPUSubtraction/src/PFMEtSignInterfaceBase.cc

0001 #include "RecoMET/METPUSubtraction/interface/PFMEtSignInterfaceBase.h"
0002
0003 #include "FWCore/Utilities/interface/Exception.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005
0006 #include "RecoMET/METAlgorithms/interface/significanceAlgo.h"
0007
0008 #include <TVectorD.h>
0009
0010 const double defaultPFMEtResolutionX = 10.;
0011 const double defaultPFMEtResolutionY = 10.;
0012
0013 const double epsilon = 1.e-9;
0014
0015 PFMEtSignInterfaceBase::PFMEtSignInterfaceBase(const edm::ParameterSet& cfg)
0016     : pfMEtResolution_(nullptr), inputFile_(nullptr), lut_(nullptr) {
0017   pfMEtResolution_ = new metsig::SignAlgoResolutions(cfg);
0018
0019   if (cfg.exists("addJERcorr")) {
0020     edm::ParameterSet cfgJERcorr = cfg.getParameter<edm::ParameterSet>("addJERcorr");
0021     edm::FileInPath inputFileName = cfgJERcorr.getParameter<edm::FileInPath>("inputFileName");
0022     std::string lutName = cfgJERcorr.getParameter<std::string>("lutName");
0023     if (inputFileName.location() != edm::FileInPath::Local)
0024       throw cms::Exception("PFMEtSignInterfaceBase") << " Failed to find File = " << inputFileName << " !!\n";
0025
0026     inputFile_ = new TFile(inputFileName.fullPath().data());
0027     lut_ = dynamic_cast<TH2*>(inputFile_->Get(lutName.data()));
0028     if (!lut_)
0029       throw cms::Exception("PFMEtSignInterfaceBase") << " Failed to load LUT = " << lutName.data()
0030                                                      << " from file = " << inputFileName.fullPath().data() << " !!\n";
0031   }
0032
0033   verbosity_ = cfg.exists("verbosity") ? cfg.getParameter<int>("verbosity") : 0;
0034 }
0035
0036 PFMEtSignInterfaceBase::~PFMEtSignInterfaceBase() {
0037   delete pfMEtResolution_;
0038   delete inputFile_;
0039   delete lut_;
0040 }
0041
0042 reco::METCovMatrix PFMEtSignInterfaceBase::operator()(const std::vector<metsig::SigInputObj>& pfMEtSignObjects) const {
0043   // if ( this->verbosity_ ) {
0044   //   std::cout << "<PFMEtSignInterfaceBase::operator()>:" << std::endl;
0045   //   std::cout << " pfMEtSignObjects: entries = " << pfMEtSignObjects.size() << std::endl;
0046   //   double dpt2Sum = 0.;
0047   //   for ( std::vector<metsig::SigInputObj>::const_iterator pfMEtSignObject = pfMEtSignObjects.begin();
0048   //      pfMEtSignObject != pfMEtSignObjects.end(); ++pfMEtSignObject ) {
0049   //     std::cout << pfMEtSignObject->get_type() << ": pt = " << pfMEtSignObject->get_energy() << ","
0050   //        << " phi = " << pfMEtSignObject->get_phi() << " --> dpt = " << pfMEtSignObject->get_sigma_e() << std::endl;
0051   //     dpt2Sum += pfMEtSignObject->get_sigma_e();
0052   //   }
0053   //   std::cout << "--> sqrt(sum(dpt^2)) = " << sqrt(dpt2Sum) << std::endl;
0054   // }
0055
0056   reco::METCovMatrix pfMEtCov;
0057   if (pfMEtSignObjects.size() >= 2) {
0058     metsig::significanceAlgo pfMEtSignAlgorithm;
0059     pfMEtSignAlgorithm.addObjects(pfMEtSignObjects);
0060     pfMEtCov = pfMEtSignAlgorithm.getSignifMatrix();
0061
0062     double det = 0;
0063     pfMEtCov.Det(det);
0064
0065     // if ( this->verbosity_ && std::abs(det) > epsilon ) {
0066     //   //keep TMatrixD as it is much easier to find
0067     //   //eigenvectors and values than with SMatrix;
0068     //   //not used anyway, except for debugging
0069     //   TMatrixD tmpMatrix(2,2);
0070     //   tmpMatrix(0,0) = pfMEtCov(0,0);
0071     //   tmpMatrix(0,1) = pfMEtCov(0,1);
0072     //   tmpMatrix(1,0) = pfMEtCov(1,0);
0073     //   tmpMatrix(1,1) = pfMEtCov(1,1);
0074
0075     //   TVectorD eigenValues(2);
0076     //   TMatrixD eigenVectors = tmpMatrix.EigenVectors(eigenValues);
0077     //   // CV: eigenvectors are stored in columns
0078     //   //     and are sorted such that the one corresponding to the highest eigenvalue is in the **first** column
0079     //   for ( unsigned iEigenVector = 0; iEigenVector < 2; ++iEigenVector ) {
0080     //  std::cout << "eigenVector #" << iEigenVector << " (eigenValue = " << eigenValues(iEigenVector) << "):"
0081     //        << " x = " << eigenVectors(0, iEigenVector) << ", y = " << eigenVectors(1, iEigenVector) << std::endl;
0082     //   }
0083     // }
0084
0085     //--- substitute (PF)MEt resolution matrix by default values
0086     //    in case resolution matrix cannot be inverted
0087     if (std::abs(det) < epsilon) {
0088       edm::LogWarning("PFMEtSignInterfaceBase::operator()")
0089           << "Inversion of PFMEt covariance matrix failed, det = " << det
0090           << " --> replacing covariance matrix by resolution defaults !!";
0091       pfMEtCov(0, 0) = defaultPFMEtResolutionX * defaultPFMEtResolutionX;
0092       pfMEtCov(0, 1) = 0.;
0093       pfMEtCov(1, 0) = 0.;
0094       pfMEtCov(1, 1) = defaultPFMEtResolutionY * defaultPFMEtResolutionY;
0095     }
0096   } else {
0097     pfMEtCov(0, 0) = 0.;
0098     pfMEtCov(0, 1) = 0.;
0099     pfMEtCov(1, 0) = 0.;
0100     pfMEtCov(1, 1) = 0.;
0101   }
0102
0103   return pfMEtCov;
0104 }