Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-10-08 23:09:58 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!