TopKinFitter/interface/CovarianceMatrix.h

0001 #ifndef CovarianceMatrix_h
0002 #define CovarianceMatrix_h
0003
0004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006
0007 #include "DataFormats/PatCandidates/interface/PATObject.h"
0008 #include "DataFormats/PatCandidates/interface/Electron.h"
0009 #include "DataFormats/PatCandidates/interface/Muon.h"
0010 #include "DataFormats/PatCandidates/interface/Jet.h"
0011 #include "DataFormats/PatCandidates/interface/MET.h"
0012
0013 #include "TopQuarkAnalysis/TopKinFitter/interface/TopKinFitter.h"
0014 #include "TopQuarkAnalysis/TopObjectResolutions/interface/MET.h"
0015 #include "TopQuarkAnalysis/TopObjectResolutions/interface/Jet.h"
0016 #include "TopQuarkAnalysis/TopObjectResolutions/interface/Muon.h"
0017 #include "TopQuarkAnalysis/TopObjectResolutions/interface/Electron.h"
0018
0019 /*
0020   \class   CovarianceMatrix CovarianceMatrix.h "TopQuarkAnalysis/TopKinFitter/interface/CovarianceMatrix.h"
0021
0022   \brief   Helper class used to setup covariance matrices for given objects and known resolutions
0023
0024   More details to be added here...
0025
0026 **/
0027
0028 class CovarianceMatrix {
0029 public:
0030   enum ObjectType { kUdscJet, kBJet, kMuon, kElectron, kMet };
0031
0032   /// default constructor
0033   CovarianceMatrix() {}
0034   /// constructor for the fully-hadronic channel
0035   CovarianceMatrix(const std::vector<edm::ParameterSet>& udscResolutions,
0036                    const std::vector<edm::ParameterSet>& bResolutions,
0037                    const std::vector<double>& jetEnergyResolutionScaleFactors,
0038                    const std::vector<double>& jetEnergyResolutionEtaBinning);
0039   /// constructor for the lepton+jets channel
0040   CovarianceMatrix(const std::vector<edm::ParameterSet>& udscResolutions,
0041                    const std::vector<edm::ParameterSet>& bResolutions,
0042                    const std::vector<edm::ParameterSet>& lepResolutions,
0043                    const std::vector<edm::ParameterSet>& metResolutions,
0044                    const std::vector<double>& jetEnergyResolutionScaleFactors,
0045                    const std::vector<double>& jetEnergyResolutionEtaBinning);
0046   // destructor
0047   ~CovarianceMatrix() {}
0048
0049   /// return covariance matrix for a PAT object
0050   template <class T>
0051   TMatrixD setupMatrix(const pat::PATObject<T>& object,
0052                        const TopKinFitter::Param param,
0053                        const std::string& resolutionProvider = "");
0054   /// return covariance matrix for a plain 4-vector
0055   TMatrixD setupMatrix(const TLorentzVector& object, const ObjectType objType, const TopKinFitter::Param param);
0056   /// get resolution for a given component of an object
0057   double getResolution(const TLorentzVector& object, const ObjectType objType, const std::string& whichResolution = "");
0058   /// get resolution for a given PAT object
0059   template <class T>
0060   double getResolution(const pat::PATObject<T>& object, const std::string& whichResolution, const bool isBJet = false) {
0061     return getResolution(TLorentzVector(object.px(), object.py(), object.pz(), object.energy()),
0062                          getObjectType(object, isBJet),
0063                          whichResolution);
0064   }
0065
0066 private:
0067   /// vector of strings for the binning of the resolutions
0068   std::vector<std::string> binsUdsc_, binsB_, binsLep_, binsMet_;
0069   /// vectors for the resolution functions
0070   std::vector<std::string> funcEtUdsc_, funcEtB_, funcEtLep_, funcEtMet_;
0071   std::vector<std::string> funcEtaUdsc_, funcEtaB_, funcEtaLep_, funcEtaMet_;
0072   std::vector<std::string> funcPhiUdsc_, funcPhiB_, funcPhiLep_, funcPhiMet_;
0073   /// scale factors for the jet energy resolution
0074   const std::vector<double> jetEnergyResolutionScaleFactors_;
0075   const std::vector<double> jetEnergyResolutionEtaBinning_;
0076
0077   /// determine type for a given PAT object
0078   template <class T>
0079   ObjectType getObjectType(const pat::PATObject<T>& object, const bool isBJet = false);
0080   /// get eta dependent smear factor for a PAT object
0081   template <class T>
0082   double getEtaDependentScaleFactor(const pat::PATObject<T>& object);
0083   /// get eta-dependent scale factor for a plain 4-vector
0084   double getEtaDependentScaleFactor(const TLorentzVector& object);
0085 };
0086
0087 template <class T>
0088 TMatrixD CovarianceMatrix::setupMatrix(const pat::PATObject<T>& object,
0089                                        const TopKinFitter::Param param,
0090                                        const std::string& resolutionProvider) {
0091   // This part is for pat objects with resolutions embedded
0092   if (object.hasKinResolution()) {
0093     TMatrixD CovM3(3, 3);
0094     CovM3.Zero();
0095     TMatrixD CovM4(4, 4);
0096     CovM4.Zero();
0097     TMatrixD* CovM = &CovM3;
0098     switch (param) {
0099       case TopKinFitter::kEtEtaPhi:
0100         CovM3(0, 0) = pow(object.resolEt(resolutionProvider), 2);
0101         if constexpr (std::is_same_v<T, reco::Jet>)
0102           CovM3(0, 0) *= getEtaDependentScaleFactor(object);
0103         else if constexpr (std::is_same_v<T, reco::MET>)
0104           CovM3(1, 1) = pow(9999., 2);
0105         else
0106           CovM3(1, 1) = pow(object.resolEta(resolutionProvider), 2);
0107         CovM3(2, 2) = pow(object.resolPhi(resolutionProvider), 2);
0108         CovM = &CovM3;
0109         break;
0110       case TopKinFitter::kEtThetaPhi:
0111         CovM3(0, 0) = pow(object.resolEt(resolutionProvider), 2);
0112         if constexpr (std::is_same_v<T, reco::Jet>)
0113           CovM3(0, 0) *= getEtaDependentScaleFactor(object);
0114         CovM3(1, 1) = pow(object.resolTheta(resolutionProvider), 2);
0115         CovM3(2, 2) = pow(object.resolPhi(resolutionProvider), 2);
0116         CovM = &CovM3;
0117         break;
0118       case TopKinFitter::kEMom:
0119         CovM4(0, 0) = pow(1, 2);
0120         CovM4(1, 1) = pow(1, 2);
0121         CovM4(2, 2) = pow(1, 2);
0122         CovM4(3, 3) = pow(1, 2);
0123         CovM = &CovM4;
0124         break;
0125     }
0126     return *CovM;
0127   }
0128   // This part is for objects without resolutions embedded
0129   else {
0130     const ObjectType objType = getObjectType(object, (resolutionProvider == "bjets"));
0131     const TLorentzVector p4(object.px(), object.py(), object.pz(), object.energy());
0132     return setupMatrix(p4, objType, param);
0133   }
0134 }
0135
0136 template <class T>
0137 CovarianceMatrix::ObjectType CovarianceMatrix::getObjectType(const pat::PATObject<T>& object, const bool isBJet) {
0138   ObjectType objType;
0139   // jets
0140   if constexpr (std::is_same_v<T, reco::Jet>) {
0141     if (isBJet)
0142       objType = kBJet;
0143     else
0144       objType = kUdscJet;
0145   }
0146   // muons
0147   else if constexpr (std::is_same_v<T, reco::Muon>)
0148     objType = kMuon;
0149   // electrons
0150   else if constexpr (std::is_same_v<T, reco::GsfElectron>)
0151     objType = kElectron;
0152   // MET
0153   else if constexpr (std::is_same_v<T, reco::MET>)
0154     objType = kMet;
0155   // catch anything else
0156   else
0157     throw cms::Exception("UnsupportedObject") << "The object given is not supported!\n";
0158   return objType;
0159 }
0160
0161 template <class T>
0162 double CovarianceMatrix::getEtaDependentScaleFactor(const pat::PATObject<T>& object) {
0163   double etaDependentScaleFactor = 1.;
0164   for (unsigned int i = 0; i < jetEnergyResolutionEtaBinning_.size(); i++) {
0165     if (std::abs(object.eta()) >= jetEnergyResolutionEtaBinning_[i] && jetEnergyResolutionEtaBinning_[i] >= 0.) {
0166       if (i == jetEnergyResolutionEtaBinning_.size() - 1) {
0167         edm::LogWarning("CovarianceMatrix") << "object eta (" << std::abs(object.eta()) << ") beyond last eta bin ("
0168                                             << jetEnergyResolutionEtaBinning_[i] << ") using scale factor 1.0!";
0169         etaDependentScaleFactor = 1.;
0170         break;
0171       }
0172       etaDependentScaleFactor = jetEnergyResolutionScaleFactors_[i];
0173     } else
0174       break;
0175   }
0176   return etaDependentScaleFactor;
0177 }
0178
0179 #endif