METAlgorithms/src/SignAlgoResolutions.cc

0001 #include "RecoMET/METAlgorithms/interface/SignAlgoResolutions.h"
0002 // -*- C++ -*-
0003 //
0004 // Package:    METAlgorithms
0005 // Class:      SignAlgoResolutions
0006 //
0007 /**\class METSignificance SignAlgoResolutions.cc RecoMET/METAlgorithms/src/SignAlgoResolutions.cc
0008
0009  Description: <one line class summary>
0010
0011  Implementation:
0012      <Notes on implementation>
0013 */
0014 //
0015 // Original Author:  Kyle Story, Freya Blekman (Cornell University)
0016 //         Created:  Fri Apr 18 11:58:33 CEST 2008
0017 //
0018 //
0019 #include "FWCore/Framework/interface/EventSetup.h"
0020 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0021 #include "DataFormats/TrackReco/interface/Track.h"
0022 #include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
0023 #include "DataFormats/ParticleFlowReco/interface/PFBlockElement.h"
0024
0025 #include <cmath>
0026
0027 #include <cstdlib>
0028 #include <iostream>
0029 #include <sstream>
0030 #include <string>
0031 #include <sys/stat.h>
0032
0033 metsig::SignAlgoResolutions::SignAlgoResolutions(const edm::ParameterSet &iConfig)
0034     : functionmap_(), ptResol_(nullptr), phiResol_(nullptr) {
0035   addResolutions(iConfig);
0036 }
0037
0038 double metsig::SignAlgoResolutions::eval(const resolutionType &type,
0039                                          const resolutionFunc &func,
0040                                          const double &et,
0041                                          const double &phi,
0042                                          const double &eta) const {
0043   // derive p from et and eta;
0044   double theta = 2 * atan(exp(-eta));
0045   double p = et / sin(theta);  // rough assumption: take e and p equivalent...
0046   return eval(type, func, et, phi, eta, p);
0047 }
0048
0049 double metsig::SignAlgoResolutions::eval(const resolutionType &type,
0050                                          const resolutionFunc &func,
0051                                          const double &et,
0052                                          const double &phi,
0053                                          const double &eta,
0054                                          const double &p) const {
0055   functionPars x(4);
0056   x[0] = et;
0057   x[1] = phi;
0058   x[2] = eta;
0059   x[3] = p;
0060   //  std::cout << "getting function of type " << type << " " << func << " " << x[0] << " " << x[1] << " " << x[2] << " " << x[3] << std::endl;
0061   return getfunc(type, func, x);
0062 }
0063
0064 metsig::SigInputObj metsig::SignAlgoResolutions::evalPF(const reco::PFCandidate *candidate) const {
0065   double eta = candidate->eta();
0066   double phi = candidate->phi();
0067   double et = candidate->energy() * sin(candidate->theta());
0068   resolutionType thetype;
0069   std::string name;
0070   int type = candidate->particleId();
0071   switch (type) {
0072     case 1:
0073       thetype = PFtype1;
0074       name = "PFChargedHadron";
0075       break;
0076     case 2:
0077       thetype = PFtype2;
0078       name = "PFChargedEM";
0079       break;
0080     case 3:
0081       thetype = PFtype3;
0082       name = "PFMuon";
0083       break;
0084     case 4:
0085       thetype = PFtype4;
0086       name = "PFNeutralEM";
0087       break;
0088     case 5:
0089       thetype = PFtype5;
0090       name = "PFNeutralHadron";
0091       break;
0092     case 6:
0093       thetype = PFtype6;
0094       name = "PFtype6";
0095       break;
0096     case 7:
0097       thetype = PFtype7;
0098       name = "PFtype7";
0099       break;
0100     default:
0101       thetype = PFtype7;
0102       name = "PFunknown";
0103       break;
0104   }
0105
0106   double d_et = 0, d_phi = 0;  //d_phi here is the error on phi component of the et
0107   reco::TrackRef trackRef = candidate->trackRef();
0108   if (!trackRef.isNull()) {
0109     d_phi = et * trackRef->phiError();
0110     d_et = (type == 2) ? ElectronPtResolution(candidate) : trackRef->ptError();
0111     //if(type==2) std::cout << eval(thetype,ET,et,phi,eta) << "    " << trackRef->ptError() << "    "<< ElectronPtResolution(candidate) << std::endl;
0112   } else {
0113     d_et = eval(thetype, ET, et, phi, eta);
0114     d_phi = eval(thetype, PHI, et, phi, eta);
0115   }
0116
0117   metsig::SigInputObj resultingobj(name, et, phi, d_et, d_phi);
0118   return resultingobj;
0119 }
0120
0121 metsig::SigInputObj metsig::SignAlgoResolutions::evalPFJet(const reco::Jet *jet) const {
0122   double jpt = jet->pt();
0123   double jphi = jet->phi();
0124   double jeta = jet->eta();
0125   double jdeltapt = 999.;
0126   double jdeltapphi = 999.;
0127
0128   if (jpt < ptResolThreshold_ && jpt < 20.) {  //use temporary fix for low pT jets
0129     double feta = TMath::Abs(jeta);
0130     int ieta = feta < 5. ? int(feta / 0.5) : 9;  //bin size = 0.5
0131     int ipt = jpt > 3. ? int(jpt - 3. / 2) : 0;  //bin size =2, starting from ptmin=3GeV
0132     jdeltapt = jdpt[ieta][ipt];
0133     jdeltapphi = jpt * jdphi[ieta][ipt];
0134   } else {
0135     //use the resolution functions at |eta|=5 to avoid crash for jets with large eta.
0136     if (jeta > 5)
0137       jeta = 5;
0138     if (jeta < -5)
0139       jeta = -5;
0140     double jptForEval = jpt > ptResolThreshold_ ? jpt : ptResolThreshold_;
0141     jdeltapt = jpt * ptResol_->parameterEtaEval("sigma", jeta, jptForEval);
0142     jdeltapphi = jpt * phiResol_->parameterEtaEval("sigma", jeta, jptForEval);
0143   }
0144
0145   std::string inputtype = "jet";
0146   metsig::SigInputObj obj_jet(inputtype, jpt, jphi, jdeltapt, jdeltapphi);
0147   //std::cout << "RESOLUTIONS JET: " << jpt << "   " << jphi<< "   " <<jdeltapt << "   " << jdeltapphi << std::endl;
0148   return obj_jet;
0149 }
0150
0151 void metsig::SignAlgoResolutions::addResolutions(const edm::ParameterSet &iConfig) {
0152   using namespace std;
0153
0154   // Jet Resolutions - for now load from the files. Migrate to EventSetup asap.
0155   metsig::SignAlgoResolutions::initializeJetResolutions(iConfig);
0156
0157   ptResolThreshold_ = iConfig.getParameter<double>("ptresolthreshold");
0158
0159   //get temporary low pT pfjet resolutions
0160   for (int ieta = 0; ieta < 10; ieta++) {
0161     jdpt[ieta] = iConfig.getParameter<std::vector<double> >(Form("jdpt%d", ieta));
0162     jdphi[ieta] = iConfig.getParameter<std::vector<double> >(Form("jdphi%d", ieta));
0163   }
0164
0165   // for now: do this by hand - this can obviously also be done via ESSource etc.
0166   functionPars etparameters(3, 0);
0167   functionPars phiparameters(1, 0);
0168   // set the parameters per function:
0169   // ECAL, BARREL:
0170   std::vector<double> ebet = iConfig.getParameter<std::vector<double> >("EB_EtResPar");
0171   std::vector<double> ebphi = iConfig.getParameter<std::vector<double> >("EB_PhiResPar");
0172
0173   etparameters[0] = ebet[0];
0174   etparameters[1] = ebet[1];
0175   etparameters[2] = ebet[2];
0176   phiparameters[0] = ebphi[0];
0177   addfunction(caloEB, ET, etparameters);
0178   addfunction(caloEB, PHI, phiparameters);
0179   // ECAL, ENDCAP:
0180   std::vector<double> eeet = iConfig.getParameter<std::vector<double> >("EE_EtResPar");
0181   std::vector<double> eephi = iConfig.getParameter<std::vector<double> >("EE_PhiResPar");
0182
0183   etparameters[0] = eeet[0];
0184   etparameters[1] = eeet[1];
0185   etparameters[2] = eeet[2];
0186   phiparameters[0] = eephi[0];
0187   addfunction(caloEE, ET, etparameters);
0188   addfunction(caloEE, PHI, phiparameters);
0189   // HCAL, BARREL:
0190   std::vector<double> hbet = iConfig.getParameter<std::vector<double> >("HB_EtResPar");
0191   std::vector<double> hbphi = iConfig.getParameter<std::vector<double> >("HB_PhiResPar");
0192
0193   etparameters[0] = hbet[0];
0194   etparameters[1] = hbet[1];
0195   etparameters[2] = hbet[2];
0196   phiparameters[0] = hbphi[0];
0197   addfunction(caloHB, ET, etparameters);
0198   addfunction(caloHB, PHI, phiparameters);
0199   // HCAL, ENDCAP:
0200   std::vector<double> heet = iConfig.getParameter<std::vector<double> >("HE_EtResPar");
0201   std::vector<double> hephi = iConfig.getParameter<std::vector<double> >("HE_PhiResPar");
0202
0203   etparameters[0] = heet[0];
0204   etparameters[1] = heet[1];
0205   etparameters[2] = heet[2];
0206   phiparameters[0] = hephi[0];
0207   addfunction(caloHE, ET, etparameters);
0208   addfunction(caloHE, PHI, phiparameters);
0209   // HCAL, Outer
0210   std::vector<double> hoet = iConfig.getParameter<std::vector<double> >("HO_EtResPar");
0211   std::vector<double> hophi = iConfig.getParameter<std::vector<double> >("HO_PhiResPar");
0212
0213   etparameters[0] = hoet[0];
0214   etparameters[1] = hoet[1];
0215   etparameters[2] = hoet[2];
0216   phiparameters[0] = hophi[0];
0217   addfunction(caloHO, ET, etparameters);
0218   addfunction(caloHO, PHI, phiparameters);
0219   // HCAL, Forward
0220   std::vector<double> hfet = iConfig.getParameter<std::vector<double> >("HF_EtResPar");
0221   std::vector<double> hfphi = iConfig.getParameter<std::vector<double> >("HF_PhiResPar");
0222
0223   etparameters[0] = hfet[0];
0224   etparameters[1] = hfet[1];
0225   etparameters[2] = hfet[2];
0226   phiparameters[0] = hfphi[0];
0227   addfunction(caloHF, ET, etparameters);
0228   addfunction(caloHF, PHI, phiparameters);
0229
0230   // PF objects:
0231   // type 1:
0232   std::vector<double> pf1et = iConfig.getParameter<std::vector<double> >("PF_EtResType1");
0233   std::vector<double> pf1phi = iConfig.getParameter<std::vector<double> >("PF_PhiResType1");
0234   etparameters[0] = pf1et[0];
0235   etparameters[1] = pf1et[1];
0236   etparameters[2] = pf1et[2];
0237   phiparameters[0] = pf1phi[0];
0238   addfunction(PFtype1, ET, etparameters);
0239   addfunction(PFtype1, PHI, phiparameters);
0240
0241   // PF objects:
0242   // type 2:
0243   std::vector<double> pf2et = iConfig.getParameter<std::vector<double> >("PF_EtResType2");
0244   std::vector<double> pf2phi = iConfig.getParameter<std::vector<double> >("PF_PhiResType2");
0245   etparameters[0] = pf2et[0];
0246   etparameters[1] = pf2et[1];
0247   etparameters[2] = pf2et[2];
0248   phiparameters[0] = pf2phi[0];
0249   addfunction(PFtype2, ET, etparameters);
0250   addfunction(PFtype2, PHI, phiparameters);
0251
0252   // PF objects:
0253   // type 3:
0254   std::vector<double> pf3et = iConfig.getParameter<std::vector<double> >("PF_EtResType3");
0255   std::vector<double> pf3phi = iConfig.getParameter<std::vector<double> >("PF_PhiResType3");
0256   etparameters[0] = pf3et[0];
0257   etparameters[1] = pf3et[1];
0258   etparameters[2] = pf3et[2];
0259   phiparameters[0] = pf3phi[0];
0260   addfunction(PFtype3, ET, etparameters);
0261   addfunction(PFtype3, PHI, phiparameters);
0262
0263   // PF objects:
0264   // type 4:
0265   std::vector<double> pf4et = iConfig.getParameter<std::vector<double> >("PF_EtResType4");
0266   std::vector<double> pf4phi = iConfig.getParameter<std::vector<double> >("PF_PhiResType4");
0267   etparameters[0] = pf4et[0];
0268   etparameters[1] = pf4et[1];
0269   etparameters[2] = pf4et[2];
0270   //phiparameters[0]=pf4phi[0];
0271   addfunction(PFtype4, ET, etparameters);
0272   addfunction(PFtype4, PHI, pf4phi);  //use the same functional form for photon phi error as for pT, pass whole vector
0273
0274   // PF objects:
0275   // type 5:
0276   std::vector<double> pf5et = iConfig.getParameter<std::vector<double> >("PF_EtResType5");
0277   std::vector<double> pf5phi = iConfig.getParameter<std::vector<double> >("PF_PhiResType5");
0278   etparameters[0] = pf5et[0];
0279   etparameters[1] = pf5et[1];
0280   etparameters[2] = pf5et[2];
0281   phiparameters[0] = pf5phi[0];
0282   addfunction(PFtype5, ET, etparameters);
0283   addfunction(PFtype5, PHI, pf5phi);
0284
0285   // PF objects:
0286   // type 6:
0287   std::vector<double> pf6et = iConfig.getParameter<std::vector<double> >("PF_EtResType6");
0288   std::vector<double> pf6phi = iConfig.getParameter<std::vector<double> >("PF_PhiResType6");
0289   etparameters[0] = pf6et[0];
0290   etparameters[1] = pf6et[1];
0291   etparameters[2] = pf6et[2];
0292   phiparameters[0] = pf6phi[0];
0293   addfunction(PFtype6, ET, etparameters);
0294   addfunction(PFtype6, PHI, phiparameters);
0295
0296   // PF objects:
0297   // type 7:
0298   std::vector<double> pf7et = iConfig.getParameter<std::vector<double> >("PF_EtResType7");
0299   std::vector<double> pf7phi = iConfig.getParameter<std::vector<double> >("PF_PhiResType7");
0300   etparameters[0] = pf7et[0];
0301   etparameters[1] = pf7et[1];
0302   etparameters[2] = pf7et[2];
0303   phiparameters[0] = pf7phi[0];
0304   addfunction(PFtype7, ET, etparameters);
0305   addfunction(PFtype7, PHI, phiparameters);
0306
0307   return;
0308 }
0309
0310 void metsig::SignAlgoResolutions::addfunction(resolutionType type,
0311                                               resolutionFunc func,
0312                                               const functionPars &parameters) {
0313   functionCombo mypair(type, func);
0314   functionmap_[mypair] = parameters;
0315 }
0316
0317 double metsig::SignAlgoResolutions::getfunc(const metsig::resolutionType &type,
0318                                             const metsig::resolutionFunc &func,
0319                                             functionPars &x) const {
0320   double result = 0;
0321   functionCombo mypair(type, func);
0322
0323   if (functionmap_.count(mypair) == 0) {
0324     return result;
0325   }
0326
0327   functionPars values = (functionmap_.find(mypair))->second;
0328   switch (func) {
0329     case metsig::ET:
0330       return EtFunction(x, values);
0331     case metsig::PHI:
0332       return PhiFunction(x, values);
0333     case metsig::TRACKP:
0334       return PFunction(x, values);
0335     case metsig::CONSTPHI:
0336       return PhiConstFunction(x, values);
0337   }
0338
0339   //  std::cout << "returning function " << type << " " << func << " " << result << " " << x[0] << std::endl;
0340
0341   return result;
0342 }
0343
0344 double metsig::SignAlgoResolutions::EtFunction(const functionPars &x, const functionPars &par) const {
0345   if (par.size() < 3)
0346     return 0.;
0347   if (x.empty())
0348     return 0.;
0349   double et = x[0];
0350   if (et <= 0.)
0351     return 0.;
0352   double result = et * sqrt((par[2] * par[2]) + (par[1] * par[1] / et) + (par[0] * par[0] / (et * et)));
0353   return result;
0354 }
0355
0356 double metsig::SignAlgoResolutions::PhiFunction(const functionPars &x, const functionPars &par) const {
0357   double et = x[0];
0358   if (et <= 0.) {
0359     return 0.;
0360   }
0361
0362   //if 1 parameter is C provided, returns C*pT, if three parameters N, S, C are provided, it returns the usual resolution value, as for sigmaPt
0363   if (par.size() != 1 && par.size() != 3) {  //only 1 or 3 parameters supported for phi function
0364     return 0.;
0365   } else if (par.size() == 1) {
0366     return par[0] * et;
0367   } else {
0368     return et * sqrt((par[2] * par[2]) + (par[1] * par[1] / et) + (par[0] * par[0] / (et * et)));
0369   }
0370 }
0371 double metsig::SignAlgoResolutions::PFunction(const functionPars &x, const functionPars &par) const {
0372   // not currently implemented
0373   return 0;
0374 }
0375
0376 double metsig::SignAlgoResolutions::PhiConstFunction(const functionPars &x, const functionPars &par) const {
0377   return par[0];
0378 }
0379
0380 void metsig::SignAlgoResolutions::initializeJetResolutions(const edm::ParameterSet &iConfig) {
0381   using namespace std;
0382
0383   // only reinitialize the resolutsion if the pointers are zero
0384   if (ptResol_ == nullptr) {
0385     string resolutionsAlgo = iConfig.getParameter<std::string>("resolutionsAlgo");
0386     string resolutionsEra = iConfig.getParameter<std::string>("resolutionsEra");
0387
0388     string cmssw_base(std::getenv("CMSSW_BASE"));
0389     string cmssw_release_base(std::getenv("CMSSW_RELEASE_BASE"));
0390     string path = cmssw_base + "/src/CondFormats/JetMETObjects/data";
0391     struct stat st;
0392     if (stat(path.c_str(), &st) != 0) {
0393       path = cmssw_release_base + "/src/CondFormats/JetMETObjects/data";
0394     }
0395     if (stat(path.c_str(), &st) != 0) {
0396       cerr << "ERROR: tried to set path but failed, abort." << endl;
0397     }
0398     const string &era(resolutionsEra);
0399     const string &alg(resolutionsAlgo);
0400     string ptFileName = path + "/" + era + "_PtResolution_" + alg + ".txt";
0401     string phiFileName = path + "/" + era + "_PhiResolution_" + alg + ".txt";
0402
0403     ptResol_ = new JetResolution(ptFileName, false);
0404     phiResol_ = new JetResolution(phiFileName, false);
0405   }
0406 }
0407
0408 double metsig::SignAlgoResolutions::ElectronPtResolution(const reco::PFCandidate *c) const {
0409   double eta = c->eta();
0410   double energy = c->energy();
0411   double dEnergy = pfresol_->getEnergyResolutionEm(energy, eta);
0412
0413   return dEnergy / cosh(eta);
0414 }