JetMETObjects/src/JetCorrectionUncertainty.cc

0001 #include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h"
0002 #include "CondFormats/JetMETObjects/interface/SimpleJetCorrectionUncertainty.h"
0003 #include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "Math/PtEtaPhiE4D.h"
0006 #include "Math/Vector3D.h"
0007 #include "Math/LorentzVector.h"
0008 #include <vector>
0009 #include <string>
0010
0011 /////////////////////////////////////////////////////////////////////////
0012 JetCorrectionUncertainty::JetCorrectionUncertainty() {
0013   mJetEta = -9999;
0014   mJetPt = -9999;
0015   mJetPhi = -9999;
0016   mJetE = -9999;
0017   mJetEMF = -9999;
0018   mLepPx = -9999;
0019   mLepPy = -9999;
0020   mLepPz = -9999;
0021   mIsJetEset = false;
0022   mIsJetPtset = false;
0023   mIsJetPhiset = false;
0024   mIsJetEtaset = false;
0025   mIsJetEMFset = false;
0026   mIsLepPxset = false;
0027   mIsLepPyset = false;
0028   mIsLepPzset = false;
0029   mAddLepToJet = false;
0030 }
0031 /////////////////////////////////////////////////////////////////////////
0032 JetCorrectionUncertainty::JetCorrectionUncertainty(const std::string& fDataFile) : mUncertainty{fDataFile} {
0033   mJetEta = -9999;
0034   mJetPt = -9999;
0035   mJetPhi = -9999;
0036   mJetE = -9999;
0037   mJetEMF = -9999;
0038   mLepPx = -9999;
0039   mLepPy = -9999;
0040   mLepPz = -9999;
0041   mIsJetEset = false;
0042   mIsJetPtset = false;
0043   mIsJetPhiset = false;
0044   mIsJetEtaset = false;
0045   mIsJetEMFset = false;
0046   mIsLepPxset = false;
0047   mIsLepPyset = false;
0048   mIsLepPzset = false;
0049   mAddLepToJet = false;
0050 }
0051 /////////////////////////////////////////////////////////////////////////
0052 JetCorrectionUncertainty::JetCorrectionUncertainty(const JetCorrectorParameters& fParameters)
0053     : mUncertainty{fParameters} {
0054   mJetEta = -9999;
0055   mJetPt = -9999;
0056   mJetPhi = -9999;
0057   mJetE = -9999;
0058   mJetEMF = -9999;
0059   mLepPx = -9999;
0060   mLepPy = -9999;
0061   mLepPz = -9999;
0062   mIsJetEset = false;
0063   mIsJetPtset = false;
0064   mIsJetPhiset = false;
0065   mIsJetEtaset = false;
0066   mIsJetEMFset = false;
0067   mIsLepPxset = false;
0068   mIsLepPyset = false;
0069   mIsLepPzset = false;
0070   mAddLepToJet = false;
0071 }
0072 /////////////////////////////////////////////////////////////////////////
0073 void JetCorrectionUncertainty::setParameters(const std::string& fDataFile) {
0074   mUncertainty = SimpleJetCorrectionUncertainty(fDataFile);
0075 }
0076 /////////////////////////////////////////////////////////////////////////
0077 float JetCorrectionUncertainty::getUncertainty(bool fDirection) {
0078   float result;
0079   std::vector<float> vx, vy;
0080   vx = fillVector(mUncertainty.parameters().definitions().binVar());
0081   vy = fillVector(mUncertainty.parameters().definitions().parVar());
0082   result = mUncertainty.uncertainty(vx, vy[0], fDirection);
0083   mIsJetEset = false;
0084   mIsJetPtset = false;
0085   mIsJetPhiset = false;
0086   mIsJetEtaset = false;
0087   mIsJetEMFset = false;
0088   mIsLepPxset = false;
0089   mIsLepPyset = false;
0090   mIsLepPzset = false;
0091   return result;
0092 }
0093 //------------------------------------------------------------------------
0094 //--- Reads the parameter names and fills a vector of floats -------------
0095 //------------------------------------------------------------------------
0096 std::vector<float> JetCorrectionUncertainty::fillVector(const std::vector<std::string>& fNames) {
0097   std::vector<float> result;
0098   for (unsigned i = 0; i < fNames.size(); i++) {
0099     if (fNames[i] == "JetEta") {
0100       if (!mIsJetEtaset) {
0101         edm::LogError("JetCorrectionUncertainty::") << " jet eta is not set";
0102         result.push_back(-999.0);
0103       } else {
0104         result.push_back(mJetEta);
0105       }
0106     } else if (fNames[i] == "JetPt") {
0107       if (!mIsJetPtset) {
0108         edm::LogError("JetCorrectionUncertainty::") << " jet pt is not set";
0109         result.push_back(-999.0);
0110       } else {
0111         result.push_back(mJetPt);
0112       }
0113     } else if (fNames[i] == "JetPhi") {
0114       if (!mIsJetPhiset) {
0115         edm::LogError("JetCorrectionUncertainty::") << " jet phi is not set";
0116         result.push_back(-999.0);
0117       } else {
0118         result.push_back(mJetPhi);
0119       }
0120     } else if (fNames[i] == "JetE") {
0121       if (!mIsJetEset) {
0122         edm::LogError("JetCorrectionUncertainty::") << " jet energy is not set";
0123         result.push_back(-999.0);
0124       } else {
0125         result.push_back(mJetE);
0126       }
0127     } else if (fNames[i] == "JetEMF") {
0128       if (!mIsJetEMFset) {
0129         edm::LogError("JetCorrectionUncertainty::") << " jet emf is not set";
0130         result.push_back(-999.0);
0131       } else {
0132         result.push_back(mJetEMF);
0133       }
0134     } else if (fNames[i] == "LepPx") {
0135       if (!mIsLepPxset) {
0136         edm::LogError("JetCorrectionUncertainty::") << " lepton px is not set";
0137         result.push_back(-999.0);
0138       } else {
0139         result.push_back(mLepPx);
0140       }
0141     } else if (fNames[i] == "LepPy") {
0142       if (!mIsLepPyset) {
0143         edm::LogError("JetCorrectionUncertainty::") << " lepton py is not set";
0144         result.push_back(-999.0);
0145       } else {
0146         result.push_back(mLepPy);
0147       }
0148     } else if (fNames[i] == "LepPz") {
0149       if (!mIsLepPzset) {
0150         edm::LogError("JetCorrectionUncertainty::") << " lepton pz is not set";
0151         result.push_back(-999.0);
0152       } else {
0153         result.push_back(mLepPz);
0154       }
0155     }
0156
0157     else {
0158       edm::LogError("JetCorrectionUncertainty::") << " unknown parameter " << fNames[i];
0159       result.push_back(-999.0);
0160     }
0161   }
0162   return result;
0163 }
0164 //------------------------------------------------------------------------
0165 //--- Calculate the PtRel (needed for the SLB) ---------------------------
0166 //------------------------------------------------------------------------
0167 float JetCorrectionUncertainty::getPtRel() {
0168   typedef ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiE4D<float> > PtEtaPhiELorentzVector;
0169   typedef ROOT::Math::DisplacementVector3D<ROOT::Math::Cartesian3D<float> > XYZVector;
0170   PtEtaPhiELorentzVector jet;
0171   XYZVector lep;
0172   jet.SetPt(mJetPt);
0173   jet.SetEta(mJetEta);
0174   jet.SetPhi(mJetPhi);
0175   jet.SetE(mJetE);
0176   lep.SetXYZ(mLepPx, mLepPy, mLepPz);
0177   float lj_x = (mAddLepToJet) ? lep.X() + jet.Px() : jet.Px();
0178   float lj_y = (mAddLepToJet) ? lep.Y() + jet.Py() : jet.Py();
0179   float lj_z = (mAddLepToJet) ? lep.Z() + jet.Pz() : jet.Pz();
0180   // absolute values squared
0181   float lj2 = lj_x * lj_x + lj_y * lj_y + lj_z * lj_z;
0182   float pTrel2 = -999.0;
0183   if (lj2 > 0) {
0184     float lep2 = lep.X() * lep.X() + lep.Y() * lep.Y() + lep.Z() * lep.Z();
0185     // projection vec(mu) to lepjet axis
0186     float lepXlj = lep.X() * lj_x + lep.Y() * lj_y + lep.Z() * lj_z;
0187     // absolute value squared and normalized
0188     float pLrel2 = lepXlj * lepXlj / lj2;
0189     // lep2 = pTrel2 + pLrel2
0190     pTrel2 = lep2 - pLrel2;
0191   } else
0192     edm::LogError("JetCorrectionUncertainty") << " not positive lepton-jet momentum: " << lj2;
0193   return (pTrel2 > 0) ? std::sqrt(pTrel2) : 0.0;
0194 }
0195 //------------------------------------------------------------------------
0196 //--- Setters ------------------------------------------------------------
0197 //------------------------------------------------------------------------
0198 void JetCorrectionUncertainty::setJetEta(float fEta) {
0199   mJetEta = fEta;
0200   mIsJetEtaset = true;
0201 }
0202 //------------------------------------------------------------------------
0203 void JetCorrectionUncertainty::setJetPt(float fPt) {
0204   mJetPt = fPt;
0205   mIsJetPtset = true;
0206 }
0207 //------------------------------------------------------------------------
0208 void JetCorrectionUncertainty::setJetPhi(float fPhi) {
0209   mJetPhi = fPhi;
0210   mIsJetPhiset = true;
0211 }
0212 //------------------------------------------------------------------------
0213 void JetCorrectionUncertainty::setJetE(float fE) {
0214   mJetE = fE;
0215   mIsJetEset = true;
0216 }
0217 //------------------------------------------------------------------------
0218 void JetCorrectionUncertainty::setJetEMF(float fEMF) {
0219   mJetEMF = fEMF;
0220   mIsJetEMFset = true;
0221 }
0222 //------------------------------------------------------------------------
0223 void JetCorrectionUncertainty::setLepPx(float fPx) {
0224   mLepPx = fPx;
0225   mIsLepPxset = true;
0226 }
0227 //------------------------------------------------------------------------
0228 void JetCorrectionUncertainty::setLepPy(float fPy) {
0229   mLepPy = fPy;
0230   mIsLepPyset = true;
0231 }
0232 //------------------------------------------------------------------------
0233 void JetCorrectionUncertainty::setLepPz(float fPz) {
0234   mLepPz = fPz;
0235   mIsLepPzset = true;
0236 }
0237 //------------------------------------------------------------------------