CovarianceMatrix.cc

Go to the documentation of this file.

#include "TopQuarkAnalysis/TopKinFitter/interface/CovarianceMatrix.h"

#include "CommonTools/Utils/interface/StringObjectFunction.h"
#include "CommonTools/Utils/interface/StringCutObjectSelector.h"

CovarianceMatrix::CovarianceMatrix(const std::vector<edm::ParameterSet>& udscResolutions, const std::vector<edm::ParameterSet>& bResolutions,
                   const std::vector<double>& jetEnergyResolutionScaleFactors, const std::vector<double>& jetEnergyResolutionEtaBinning):
  jetEnergyResolutionScaleFactors_(jetEnergyResolutionScaleFactors), jetEnergyResolutionEtaBinning_(jetEnergyResolutionEtaBinning)
{
  for(std::vector<edm::ParameterSet>::const_iterator iSet = udscResolutions.begin(); iSet != udscResolutions.end(); ++iSet){
    if(iSet->exists("bin")) binsUdsc_.push_back(iSet->getParameter<std::string>("bin"));
    else if(udscResolutions.size()==1) binsUdsc_.push_back("");
    else throw cms::Exception("Configuration") << "Parameter 'bin' is needed if more than one PSet is specified!\n";
    
    funcEtUdsc_.push_back(iSet->getParameter<std::string>("et"));
    funcEtaUdsc_.push_back(iSet->getParameter<std::string>("eta"));
    funcPhiUdsc_.push_back(iSet->getParameter<std::string>("phi"));
  }
  for(std::vector<edm::ParameterSet>::const_iterator iSet = bResolutions.begin(); iSet != bResolutions.end(); ++iSet){
    if(iSet->exists("bin")) binsB_.push_back(iSet->getParameter<std::string>("bin"));
    else if(bResolutions.size()==1) binsB_.push_back("");
    else throw cms::Exception("Configuration") << "Parameter 'bin' is needed if more than one PSet is specified!\n";
    
    funcEtB_.push_back(iSet->getParameter<std::string>("et"));
    funcEtaB_.push_back(iSet->getParameter<std::string>("eta"));
    funcPhiB_.push_back(iSet->getParameter<std::string>("phi"));
  }
}

CovarianceMatrix::CovarianceMatrix(const std::vector<edm::ParameterSet>& udscResolutions, const std::vector<edm::ParameterSet>& bResolutions,
                   const std::vector<edm::ParameterSet>& lepResolutions, const std::vector<edm::ParameterSet>& metResolutions,
                   const std::vector<double>& jetEnergyResolutionScaleFactors, const std::vector<double>& jetEnergyResolutionEtaBinning):
  jetEnergyResolutionScaleFactors_(jetEnergyResolutionScaleFactors), jetEnergyResolutionEtaBinning_(jetEnergyResolutionEtaBinning)
{
  if(jetEnergyResolutionScaleFactors_.size()+1!=jetEnergyResolutionEtaBinning_.size())
    throw cms::Exception("Configuration") << "The number of scale factors does not fit to the number of eta bins!\n";
  for(unsigned int i=0; i<jetEnergyResolutionEtaBinning_.size(); i++)
    if(jetEnergyResolutionEtaBinning_[i]<0. && i<jetEnergyResolutionEtaBinning_.size()-1)
      throw cms::Exception("Configuration") << "eta binning in absolut values required!\n";

  for(std::vector<edm::ParameterSet>::const_iterator iSet = udscResolutions.begin(); iSet != udscResolutions.end(); ++iSet){
    if(iSet->exists("bin")) binsUdsc_.push_back(iSet->getParameter<std::string>("bin"));
    else if(udscResolutions.size()==1) binsUdsc_.push_back("");
    else throw cms::Exception("Configuration") << "Parameter 'bin' is needed if more than one PSet is specified!\n";
    
    funcEtUdsc_.push_back(iSet->getParameter<std::string>("et"));
    funcEtaUdsc_.push_back(iSet->getParameter<std::string>("eta"));
    funcPhiUdsc_.push_back(iSet->getParameter<std::string>("phi"));
  }
  for(std::vector<edm::ParameterSet>::const_iterator iSet = bResolutions.begin(); iSet != bResolutions.end(); ++iSet){
    if(iSet->exists("bin")) binsB_.push_back(iSet->getParameter<std::string>("bin"));
    else if(bResolutions.size()==1) binsB_.push_back("");
    else throw cms::Exception("Configuration") << "Parameter 'bin' is needed if more than one PSet is specified!\n";
    
    funcEtB_.push_back(iSet->getParameter<std::string>("et"));
    funcEtaB_.push_back(iSet->getParameter<std::string>("eta"));
    funcPhiB_.push_back(iSet->getParameter<std::string>("phi"));
  }
  for(std::vector<edm::ParameterSet>::const_iterator iSet = lepResolutions.begin(); iSet != lepResolutions.end(); ++iSet){
    if(iSet->exists("bin")) binsLep_.push_back(iSet->getParameter<std::string>("bin"));
    else if(lepResolutions.size()==1) binsLep_.push_back("");
    else throw cms::Exception("Configuration") << "Parameter 'bin' is needed if more than one PSet is specified!\n";
    
    funcEtLep_.push_back(iSet->getParameter<std::string>("et"));
    funcEtaLep_.push_back(iSet->getParameter<std::string>("eta"));
    funcPhiLep_.push_back(iSet->getParameter<std::string>("phi"));
  }
  for(std::vector<edm::ParameterSet>::const_iterator iSet = metResolutions.begin(); iSet != metResolutions.end(); ++iSet){
    if(iSet->exists("bin")) binsMet_.push_back(iSet->getParameter<std::string>("bin"));
    else if(metResolutions.size()==1) binsMet_.push_back("");
    else throw cms::Exception("Configuration") << "Parameter 'bin' is needed if more than one PSet is specified!\n";
    
    funcEtMet_.push_back(iSet->getParameter<std::string>("et"));
    funcEtaMet_.push_back(iSet->getParameter<std::string>("eta"));
    funcPhiMet_.push_back(iSet->getParameter<std::string>("phi"));
  }
}

double CovarianceMatrix::getResolution(const TLorentzVector& object, const ObjectType objType, const std::string& whichResolution)
{
  std::vector<std::string> * bins_, * funcEt_, * funcEta_, * funcPhi_;

  switch(objType) {
  case kUdscJet:
    bins_    = &binsUdsc_;
    funcEt_  = &funcEtUdsc_;
    funcEta_ = &funcEtaUdsc_;
    funcPhi_ = &funcPhiUdsc_;
    break;
  case kBJet:
    bins_    = &binsB_;
    funcEt_  = &funcEtB_;
    funcEta_ = &funcEtaB_;
    funcPhi_ = &funcPhiB_;
    break;
  case kMuon:
    bins_    = &binsLep_;
    funcEt_  = &funcEtLep_;
    funcEta_ = &funcEtaLep_;
    funcPhi_ = &funcPhiLep_;
    break;
  case kElectron:
    bins_    = &binsLep_;
    funcEt_  = &funcEtLep_;
    funcEta_ = &funcEtaLep_;
    funcPhi_ = &funcPhiLep_;
    break;
  case kMet:
    bins_    = &binsMet_;
    funcEt_  = &funcEtMet_;
    funcEta_ = &funcEtaMet_;
    funcPhi_ = &funcPhiMet_;
    break;
  default:
    throw cms::Exception("UnsupportedObject") << "The object given is not supported!\n";
  }

  int selectedBin=-1;
  reco::LeafCandidate candidate;
  for(unsigned int i=0; i<bins_->size(); ++i){
    StringCutObjectSelector<reco::LeafCandidate> select_(bins_->at(i));
    candidate = reco::LeafCandidate( 0, reco::LeafCandidate::LorentzVector(object.Px(), object.Py(), object.Pz(), object.Energy()) );
    if(select_(candidate)){
      selectedBin = i;
      break;
    }
  }
  double res = 0;
  if(selectedBin>=0){
    if(whichResolution == "et")       res = StringObjectFunction<reco::LeafCandidate>(funcEt_ ->at(selectedBin)).operator()(candidate);
    else if(whichResolution == "eta") res = StringObjectFunction<reco::LeafCandidate>(funcEta_->at(selectedBin)).operator()(candidate);
    else if(whichResolution == "phi") res = StringObjectFunction<reco::LeafCandidate>(funcPhi_->at(selectedBin)).operator()(candidate);
    else throw cms::Exception("ProgrammingError") << "Only 'et', 'eta' and 'phi' resolutions supported!\n";
  }
  return res;
}

TMatrixD CovarianceMatrix::setupMatrix(const TLorentzVector& object, const ObjectType objType, const TopKinFitter::Param param)
{
  TMatrixD CovM3 (3,3); CovM3.Zero();
  TMatrixD CovM4 (4,4); CovM4.Zero();
  const double pt  = object.Pt();
  const double eta = object.Eta();
  switch(objType) {
  case kUdscJet:
    // if object is a non-b jet
    {
      res::HelperJet jetRes;
      switch(param) {
      case TopKinFitter::kEMom :
    CovM4(0,0) = pow(jetRes.a (pt, eta, res::HelperJet::kUds), 2);
    CovM4(1,1) = pow(jetRes.b (pt, eta, res::HelperJet::kUds), 2);
    CovM4(2,2) = pow(jetRes.c (pt, eta, res::HelperJet::kUds), 2);
    CovM4(3,3) = pow(jetRes.d (pt, eta, res::HelperJet::kUds), 2);
    return CovM4;
      case TopKinFitter::kEtEtaPhi : 
    if(binsUdsc_.empty()){
      CovM3(0,0) = pow(jetRes.et (pt, eta, res::HelperJet::kUds), 2);
      CovM3(0,0)*= pow(getEtaDependentScaleFactor(object)       , 2);
      CovM3(1,1) = pow(jetRes.eta(pt, eta, res::HelperJet::kUds), 2);
      CovM3(2,2) = pow(jetRes.phi(pt, eta, res::HelperJet::kUds), 2);
    }
    else{
      CovM3(0,0) = pow(getResolution(object, objType, "et") , 2);
      CovM3(0,0)*= pow(getEtaDependentScaleFactor(object)   , 2);
      CovM3(1,1) = pow(getResolution(object, objType, "eta"), 2);
      CovM3(2,2) = pow(getResolution(object, objType, "phi"), 2);
    }   
    return CovM3;
      case TopKinFitter::kEtThetaPhi :
    CovM3(0,0) = pow(jetRes.et   (pt, eta, res::HelperJet::kUds), 2);
    CovM3(0,0)*= pow(getEtaDependentScaleFactor(object)         , 2);
    CovM3(1,1) = pow(jetRes.theta(pt, eta, res::HelperJet::kUds), 2);
    CovM3(2,2) = pow(jetRes.phi  (pt, eta, res::HelperJet::kUds), 2);
    return CovM3;
      }
    }
    break;
  case kBJet:
    // if object is a b jet
    {
      res::HelperJet jetRes;
      switch(param) {
      case TopKinFitter::kEMom :
    CovM4(0,0) = pow(jetRes.a (pt, eta, res::HelperJet::kB), 2);
    CovM4(1,1) = pow(jetRes.b (pt, eta, res::HelperJet::kB), 2);
    CovM4(2,2) = pow(jetRes.c (pt, eta, res::HelperJet::kB), 2);
    CovM4(3,3) = pow(jetRes.d (pt, eta, res::HelperJet::kB), 2);
    return CovM4;
      case TopKinFitter::kEtEtaPhi : 
    if(binsUdsc_.empty()){
      CovM3(0,0) = pow(jetRes.et (pt, eta, res::HelperJet::kB), 2);
      CovM3(0,0)*= pow(getEtaDependentScaleFactor(object)     , 2);
      CovM3(1,1) = pow(jetRes.eta(pt, eta, res::HelperJet::kB), 2);
      CovM3(2,2) = pow(jetRes.phi(pt, eta, res::HelperJet::kB), 2);
    }
    else{
      CovM3(0,0) = pow(getResolution(object, objType, "et") , 2); 
      CovM3(0,0)*= pow(getEtaDependentScaleFactor(object)   , 2);
      CovM3(1,1) = pow(getResolution(object, objType, "eta"), 2); 
      CovM3(2,2) = pow(getResolution(object, objType, "phi"), 2);
    }
    return CovM3;
      case TopKinFitter::kEtThetaPhi :
    CovM3(0,0) = pow(jetRes.et   (pt, eta, res::HelperJet::kB), 2);
    CovM3(0,0)*= pow(getEtaDependentScaleFactor(object)       , 2);
    CovM3(1,1) = pow(jetRes.theta(pt, eta, res::HelperJet::kB), 2);
    CovM3(2,2) = pow(jetRes.phi  (pt, eta, res::HelperJet::kB), 2);
    return CovM3;
      }
    }
    break;
  case kMuon:
    // if object is a muon
    {
      res::HelperMuon muonRes;
      switch(param){
      case TopKinFitter::kEMom :
    CovM3(0,0) = pow(muonRes.a (pt, eta), 2);
    CovM3(1,1) = pow(muonRes.b (pt, eta), 2); 
    CovM3(2,2) = pow(muonRes.c (pt, eta), 2);
    return CovM3;
      case TopKinFitter::kEtEtaPhi :
    if(binsLep_.empty()){
      CovM3(0,0) = pow(muonRes.et (pt, eta), 2);
      CovM3(1,1) = pow(muonRes.eta(pt, eta), 2); 
      CovM3(2,2) = pow(muonRes.phi(pt, eta), 2);
    }
    else{
      CovM3(0,0) = pow(getResolution(object, objType, "et") , 2);
      CovM3(1,1) = pow(getResolution(object, objType, "eta"), 2);
      CovM3(2,2) = pow(getResolution(object, objType, "phi"), 2);
    }
    return CovM3;
      case TopKinFitter::kEtThetaPhi :
    CovM3(0,0) = pow(muonRes.et   (pt, eta), 2);
    CovM3(1,1) = pow(muonRes.theta(pt, eta), 2); 
    CovM3(2,2) = pow(muonRes.phi  (pt, eta), 2);
    return CovM3;
      }
    }
    break;
  case kElectron:
    {
      // if object is an electron
      res::HelperElectron elecRes;
      switch(param){
      case TopKinFitter::kEMom :
    CovM3(0,0) = pow(elecRes.a (pt, eta), 2);
    CovM3(1,1) = pow(elecRes.b (pt, eta), 2); 
    CovM3(2,2) = pow(elecRes.c (pt, eta), 2);
    return CovM3;
      case TopKinFitter::kEtEtaPhi :
    if(binsLep_.empty()){
      CovM3(0,0) = pow(elecRes.et (pt, eta), 2);
      CovM3(1,1) = pow(elecRes.eta(pt, eta), 2); 
      CovM3(2,2) = pow(elecRes.phi(pt, eta), 2);
    }
    else{
      CovM3(0,0) = pow(getResolution(object, objType, "et") , 2);
      CovM3(1,1) = pow(getResolution(object, objType, "eta"), 2);
      CovM3(2,2) = pow(getResolution(object, objType, "phi"), 2);
    }
    return CovM3;
      case TopKinFitter::kEtThetaPhi :
    CovM3(0,0) = pow(elecRes.et   (pt, eta), 2);
    CovM3(1,1) = pow(elecRes.theta(pt, eta), 2); 
    CovM3(2,2) = pow(elecRes.phi  (pt, eta), 2);
    return CovM3;
      }
    }
    break;
  case kMet:
    // if object is met
    {
      res::HelperMET metRes;
      switch(param){
      case TopKinFitter::kEMom :
    CovM3(0,0) = pow(metRes.a(pt), 2);
    CovM3(1,1) = pow(metRes.b(pt), 2);
    CovM3(2,2) = pow(metRes.c(pt), 2);
    return CovM3;
      case TopKinFitter::kEtEtaPhi :
    if(binsMet_.empty()){
      CovM3(0,0) = pow(metRes.et(pt) , 2);
      CovM3(1,1) = pow(        9999. , 2);
      CovM3(2,2) = pow(metRes.phi(pt), 2);
    }
    else{
      CovM3(0,0) = pow(getResolution(object, objType, "et") , 2);
      CovM3(1,1) = pow(getResolution(object, objType, "eta"), 2);
      CovM3(2,2) = pow(getResolution(object, objType, "phi"), 2);
    }
    return CovM3;
      case TopKinFitter::kEtThetaPhi :
    CovM3(0,0) = pow(metRes.et(pt) , 2);
    CovM3(1,1) = pow(        9999. , 2);
    CovM3(2,2) = pow(metRes.phi(pt), 2);
    return CovM3;
      }
    }
    break;
  }
  cms::Exception("Logic") << "Something went wrong while trying to setup a covariance matrix!\n";
  return CovM4; //should never get here
}

double CovarianceMatrix::getEtaDependentScaleFactor(const TLorentzVector& object)
{
  double etaDependentScaleFactor = 1.;
  for(unsigned int i=0; i<jetEnergyResolutionEtaBinning_.size(); i++){
    if(std::abs(object.Eta())>=jetEnergyResolutionEtaBinning_[i] && jetEnergyResolutionEtaBinning_[i]>=0.){
      if(i==jetEnergyResolutionEtaBinning_.size()-1) {
    edm::LogWarning("CovarianceMatrix") << "object eta ("<<std::abs(object.Eta())<<") beyond last eta bin ("<<jetEnergyResolutionEtaBinning_[i]<<") using scale factor 1.0!";
    etaDependentScaleFactor=1.;
    break;
      }
      etaDependentScaleFactor=jetEnergyResolutionScaleFactors_[i];
    }
    else
      break;
  }
  return etaDependentScaleFactor;
}