Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:26:14 

0001 #include "TopQuarkAnalysis/TopKinFitter/interface/CovarianceMatrix.h"
0002 
0003 #include "CommonTools/Utils/interface/StringObjectFunction.h"
0004 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 
0007 CovarianceMatrix::CovarianceMatrix(const std::vector<edm::ParameterSet>& udscResolutions,
0008                                    const std::vector<edm::ParameterSet>& bResolutions,
0009                                    const std::vector<double>& jetEnergyResolutionScaleFactors,
0010                                    const std::vector<double>& jetEnergyResolutionEtaBinning)
0011     : jetEnergyResolutionScaleFactors_(jetEnergyResolutionScaleFactors),
0012       jetEnergyResolutionEtaBinning_(jetEnergyResolutionEtaBinning) {
0013   for (std::vector<edm::ParameterSet>::const_iterator iSet = udscResolutions.begin(); iSet != udscResolutions.end();
0014        ++iSet) {
0015     if (iSet->exists("bin"))
0016       binsUdsc_.push_back(iSet->getParameter<std::string>("bin"));
0017     else if (udscResolutions.size() == 1)
0018       binsUdsc_.push_back("");
0019     else
0020       throw cms::Exception("Configuration") << "Parameter 'bin' is needed if more than one PSet is specified!\n";
0021 
0022     funcEtUdsc_.push_back(iSet->getParameter<std::string>("et"));
0023     funcEtaUdsc_.push_back(iSet->getParameter<std::string>("eta"));
0024     funcPhiUdsc_.push_back(iSet->getParameter<std::string>("phi"));
0025   }
0026   for (std::vector<edm::ParameterSet>::const_iterator iSet = bResolutions.begin(); iSet != bResolutions.end(); ++iSet) {
0027     if (iSet->exists("bin"))
0028       binsB_.push_back(iSet->getParameter<std::string>("bin"));
0029     else if (bResolutions.size() == 1)
0030       binsB_.push_back("");
0031     else
0032       throw cms::Exception("Configuration") << "Parameter 'bin' is needed if more than one PSet is specified!\n";
0033 
0034     funcEtB_.push_back(iSet->getParameter<std::string>("et"));
0035     funcEtaB_.push_back(iSet->getParameter<std::string>("eta"));
0036     funcPhiB_.push_back(iSet->getParameter<std::string>("phi"));
0037   }
0038 }
0039 
0040 CovarianceMatrix::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     : jetEnergyResolutionScaleFactors_(jetEnergyResolutionScaleFactors),
0047       jetEnergyResolutionEtaBinning_(jetEnergyResolutionEtaBinning) {
0048   if (jetEnergyResolutionScaleFactors_.size() + 1 != jetEnergyResolutionEtaBinning_.size())
0049     throw cms::Exception("Configuration") << "The number of scale factors does not fit to the number of eta bins!\n";
0050   for (unsigned int i = 0; i < jetEnergyResolutionEtaBinning_.size(); i++)
0051     if (jetEnergyResolutionEtaBinning_[i] < 0. && i < jetEnergyResolutionEtaBinning_.size() - 1)
0052       throw cms::Exception("Configuration") << "eta binning in absolut values required!\n";
0053 
0054   for (std::vector<edm::ParameterSet>::const_iterator iSet = udscResolutions.begin(); iSet != udscResolutions.end();
0055        ++iSet) {
0056     if (iSet->exists("bin"))
0057       binsUdsc_.push_back(iSet->getParameter<std::string>("bin"));
0058     else if (udscResolutions.size() == 1)
0059       binsUdsc_.push_back("");
0060     else
0061       throw cms::Exception("Configuration") << "Parameter 'bin' is needed if more than one PSet is specified!\n";
0062 
0063     funcEtUdsc_.push_back(iSet->getParameter<std::string>("et"));
0064     funcEtaUdsc_.push_back(iSet->getParameter<std::string>("eta"));
0065     funcPhiUdsc_.push_back(iSet->getParameter<std::string>("phi"));
0066   }
0067   for (std::vector<edm::ParameterSet>::const_iterator iSet = bResolutions.begin(); iSet != bResolutions.end(); ++iSet) {
0068     if (iSet->exists("bin"))
0069       binsB_.push_back(iSet->getParameter<std::string>("bin"));
0070     else if (bResolutions.size() == 1)
0071       binsB_.push_back("");
0072     else
0073       throw cms::Exception("Configuration") << "Parameter 'bin' is needed if more than one PSet is specified!\n";
0074 
0075     funcEtB_.push_back(iSet->getParameter<std::string>("et"));
0076     funcEtaB_.push_back(iSet->getParameter<std::string>("eta"));
0077     funcPhiB_.push_back(iSet->getParameter<std::string>("phi"));
0078   }
0079   for (std::vector<edm::ParameterSet>::const_iterator iSet = lepResolutions.begin(); iSet != lepResolutions.end();
0080        ++iSet) {
0081     if (iSet->exists("bin"))
0082       binsLep_.push_back(iSet->getParameter<std::string>("bin"));
0083     else if (lepResolutions.size() == 1)
0084       binsLep_.push_back("");
0085     else
0086       throw cms::Exception("Configuration") << "Parameter 'bin' is needed if more than one PSet is specified!\n";
0087 
0088     funcEtLep_.push_back(iSet->getParameter<std::string>("et"));
0089     funcEtaLep_.push_back(iSet->getParameter<std::string>("eta"));
0090     funcPhiLep_.push_back(iSet->getParameter<std::string>("phi"));
0091   }
0092   for (std::vector<edm::ParameterSet>::const_iterator iSet = metResolutions.begin(); iSet != metResolutions.end();
0093        ++iSet) {
0094     if (iSet->exists("bin"))
0095       binsMet_.push_back(iSet->getParameter<std::string>("bin"));
0096     else if (metResolutions.size() == 1)
0097       binsMet_.push_back("");
0098     else
0099       throw cms::Exception("Configuration") << "Parameter 'bin' is needed if more than one PSet is specified!\n";
0100 
0101     funcEtMet_.push_back(iSet->getParameter<std::string>("et"));
0102     funcEtaMet_.push_back(iSet->getParameter<std::string>("eta"));
0103     funcPhiMet_.push_back(iSet->getParameter<std::string>("phi"));
0104   }
0105 }
0106 
0107 double CovarianceMatrix::getResolution(const TLorentzVector& object,
0108                                        const ObjectType objType,
0109                                        const std::string& whichResolution) {
0110   std::vector<std::string>*bins_, *funcEt_, *funcEta_, *funcPhi_;
0111 
0112   switch (objType) {
0113     case kUdscJet:
0114       bins_ = &binsUdsc_;
0115       funcEt_ = &funcEtUdsc_;
0116       funcEta_ = &funcEtaUdsc_;
0117       funcPhi_ = &funcPhiUdsc_;
0118       break;
0119     case kBJet:
0120       bins_ = &binsB_;
0121       funcEt_ = &funcEtB_;
0122       funcEta_ = &funcEtaB_;
0123       funcPhi_ = &funcPhiB_;
0124       break;
0125     case kMuon:
0126       bins_ = &binsLep_;
0127       funcEt_ = &funcEtLep_;
0128       funcEta_ = &funcEtaLep_;
0129       funcPhi_ = &funcPhiLep_;
0130       break;
0131     case kElectron:
0132       bins_ = &binsLep_;
0133       funcEt_ = &funcEtLep_;
0134       funcEta_ = &funcEtaLep_;
0135       funcPhi_ = &funcPhiLep_;
0136       break;
0137     case kMet:
0138       bins_ = &binsMet_;
0139       funcEt_ = &funcEtMet_;
0140       funcEta_ = &funcEtaMet_;
0141       funcPhi_ = &funcPhiMet_;
0142       break;
0143     default:
0144       throw cms::Exception("UnsupportedObject") << "The object given is not supported!\n";
0145   }
0146 
0147   int selectedBin = -1;
0148   reco::LeafCandidate candidate;
0149   for (unsigned int i = 0; i < bins_->size(); ++i) {
0150     StringCutObjectSelector<reco::LeafCandidate> select_(bins_->at(i));
0151     candidate = reco::LeafCandidate(
0152         0, reco::LeafCandidate::LorentzVector(object.Px(), object.Py(), object.Pz(), object.Energy()));
0153     if (select_(candidate)) {
0154       selectedBin = i;
0155       break;
0156     }
0157   }
0158   double res = 0;
0159   if (selectedBin >= 0) {
0160     if (whichResolution == "et")
0161       res = StringObjectFunction<reco::LeafCandidate>(funcEt_->at(selectedBin)).operator()(candidate);
0162     else if (whichResolution == "eta")
0163       res = StringObjectFunction<reco::LeafCandidate>(funcEta_->at(selectedBin)).operator()(candidate);
0164     else if (whichResolution == "phi")
0165       res = StringObjectFunction<reco::LeafCandidate>(funcPhi_->at(selectedBin)).operator()(candidate);
0166     else
0167       throw cms::Exception("ProgrammingError") << "Only 'et', 'eta' and 'phi' resolutions supported!\n";
0168   }
0169   return res;
0170 }
0171 
0172 TMatrixD CovarianceMatrix::setupMatrix(const TLorentzVector& object,
0173                                        const ObjectType objType,
0174                                        const TopKinFitter::Param param) {
0175   TMatrixD CovM3(3, 3);
0176   CovM3.Zero();
0177   TMatrixD CovM4(4, 4);
0178   CovM4.Zero();
0179   const double pt = object.Pt();
0180   const double eta = object.Eta();
0181   switch (objType) {
0182     case kUdscJet:
0183       // if object is a non-b jet
0184       {
0185         res::HelperJet jetRes;
0186         switch (param) {
0187           case TopKinFitter::kEMom:
0188             CovM4(0, 0) = pow(jetRes.a(pt, eta, res::HelperJet::kUds), 2);
0189             CovM4(1, 1) = pow(jetRes.b(pt, eta, res::HelperJet::kUds), 2);
0190             CovM4(2, 2) = pow(jetRes.c(pt, eta, res::HelperJet::kUds), 2);
0191             CovM4(3, 3) = pow(jetRes.d(pt, eta, res::HelperJet::kUds), 2);
0192             return CovM4;
0193           case TopKinFitter::kEtEtaPhi:
0194             if (binsUdsc_.empty()) {
0195               CovM3(0, 0) = pow(jetRes.et(pt, eta, res::HelperJet::kUds), 2);
0196               CovM3(0, 0) *= pow(getEtaDependentScaleFactor(object), 2);
0197               CovM3(1, 1) = pow(jetRes.eta(pt, eta, res::HelperJet::kUds), 2);
0198               CovM3(2, 2) = pow(jetRes.phi(pt, eta, res::HelperJet::kUds), 2);
0199             } else {
0200               CovM3(0, 0) = pow(getResolution(object, objType, "et"), 2);
0201               CovM3(0, 0) *= pow(getEtaDependentScaleFactor(object), 2);
0202               CovM3(1, 1) = pow(getResolution(object, objType, "eta"), 2);
0203               CovM3(2, 2) = pow(getResolution(object, objType, "phi"), 2);
0204             }
0205             return CovM3;
0206           case TopKinFitter::kEtThetaPhi:
0207             CovM3(0, 0) = pow(jetRes.et(pt, eta, res::HelperJet::kUds), 2);
0208             CovM3(0, 0) *= pow(getEtaDependentScaleFactor(object), 2);
0209             CovM3(1, 1) = pow(jetRes.theta(pt, eta, res::HelperJet::kUds), 2);
0210             CovM3(2, 2) = pow(jetRes.phi(pt, eta, res::HelperJet::kUds), 2);
0211             return CovM3;
0212         }
0213       }
0214       break;
0215     case kBJet:
0216       // if object is a b jet
0217       {
0218         res::HelperJet jetRes;
0219         switch (param) {
0220           case TopKinFitter::kEMom:
0221             CovM4(0, 0) = pow(jetRes.a(pt, eta, res::HelperJet::kB), 2);
0222             CovM4(1, 1) = pow(jetRes.b(pt, eta, res::HelperJet::kB), 2);
0223             CovM4(2, 2) = pow(jetRes.c(pt, eta, res::HelperJet::kB), 2);
0224             CovM4(3, 3) = pow(jetRes.d(pt, eta, res::HelperJet::kB), 2);
0225             return CovM4;
0226           case TopKinFitter::kEtEtaPhi:
0227             if (binsUdsc_.empty()) {
0228               CovM3(0, 0) = pow(jetRes.et(pt, eta, res::HelperJet::kB), 2);
0229               CovM3(0, 0) *= pow(getEtaDependentScaleFactor(object), 2);
0230               CovM3(1, 1) = pow(jetRes.eta(pt, eta, res::HelperJet::kB), 2);
0231               CovM3(2, 2) = pow(jetRes.phi(pt, eta, res::HelperJet::kB), 2);
0232             } else {
0233               CovM3(0, 0) = pow(getResolution(object, objType, "et"), 2);
0234               CovM3(0, 0) *= pow(getEtaDependentScaleFactor(object), 2);
0235               CovM3(1, 1) = pow(getResolution(object, objType, "eta"), 2);
0236               CovM3(2, 2) = pow(getResolution(object, objType, "phi"), 2);
0237             }
0238             return CovM3;
0239           case TopKinFitter::kEtThetaPhi:
0240             CovM3(0, 0) = pow(jetRes.et(pt, eta, res::HelperJet::kB), 2);
0241             CovM3(0, 0) *= pow(getEtaDependentScaleFactor(object), 2);
0242             CovM3(1, 1) = pow(jetRes.theta(pt, eta, res::HelperJet::kB), 2);
0243             CovM3(2, 2) = pow(jetRes.phi(pt, eta, res::HelperJet::kB), 2);
0244             return CovM3;
0245         }
0246       }
0247       break;
0248     case kMuon:
0249       // if object is a muon
0250       {
0251         res::HelperMuon muonRes;
0252         switch (param) {
0253           case TopKinFitter::kEMom:
0254             CovM3(0, 0) = pow(muonRes.a(pt, eta), 2);
0255             CovM3(1, 1) = pow(muonRes.b(pt, eta), 2);
0256             CovM3(2, 2) = pow(muonRes.c(pt, eta), 2);
0257             return CovM3;
0258           case TopKinFitter::kEtEtaPhi:
0259             if (binsLep_.empty()) {
0260               CovM3(0, 0) = pow(muonRes.et(pt, eta), 2);
0261               CovM3(1, 1) = pow(muonRes.eta(pt, eta), 2);
0262               CovM3(2, 2) = pow(muonRes.phi(pt, eta), 2);
0263             } else {
0264               CovM3(0, 0) = pow(getResolution(object, objType, "et"), 2);
0265               CovM3(1, 1) = pow(getResolution(object, objType, "eta"), 2);
0266               CovM3(2, 2) = pow(getResolution(object, objType, "phi"), 2);
0267             }
0268             return CovM3;
0269           case TopKinFitter::kEtThetaPhi:
0270             CovM3(0, 0) = pow(muonRes.et(pt, eta), 2);
0271             CovM3(1, 1) = pow(muonRes.theta(pt, eta), 2);
0272             CovM3(2, 2) = pow(muonRes.phi(pt, eta), 2);
0273             return CovM3;
0274         }
0275       }
0276       break;
0277     case kElectron: {
0278       // if object is an electron
0279       res::HelperElectron elecRes;
0280       switch (param) {
0281         case TopKinFitter::kEMom:
0282           CovM3(0, 0) = pow(elecRes.a(pt, eta), 2);
0283           CovM3(1, 1) = pow(elecRes.b(pt, eta), 2);
0284           CovM3(2, 2) = pow(elecRes.c(pt, eta), 2);
0285           return CovM3;
0286         case TopKinFitter::kEtEtaPhi:
0287           if (binsLep_.empty()) {
0288             CovM3(0, 0) = pow(elecRes.et(pt, eta), 2);
0289             CovM3(1, 1) = pow(elecRes.eta(pt, eta), 2);
0290             CovM3(2, 2) = pow(elecRes.phi(pt, eta), 2);
0291           } else {
0292             CovM3(0, 0) = pow(getResolution(object, objType, "et"), 2);
0293             CovM3(1, 1) = pow(getResolution(object, objType, "eta"), 2);
0294             CovM3(2, 2) = pow(getResolution(object, objType, "phi"), 2);
0295           }
0296           return CovM3;
0297         case TopKinFitter::kEtThetaPhi:
0298           CovM3(0, 0) = pow(elecRes.et(pt, eta), 2);
0299           CovM3(1, 1) = pow(elecRes.theta(pt, eta), 2);
0300           CovM3(2, 2) = pow(elecRes.phi(pt, eta), 2);
0301           return CovM3;
0302       }
0303     } break;
0304     case kMet:
0305       // if object is met
0306       {
0307         res::HelperMET metRes;
0308         switch (param) {
0309           case TopKinFitter::kEMom:
0310             CovM3(0, 0) = pow(metRes.a(pt), 2);
0311             CovM3(1, 1) = pow(metRes.b(pt), 2);
0312             CovM3(2, 2) = pow(metRes.c(pt), 2);
0313             return CovM3;
0314           case TopKinFitter::kEtEtaPhi:
0315             if (binsMet_.empty()) {
0316               CovM3(0, 0) = pow(metRes.et(pt), 2);
0317               CovM3(1, 1) = pow(9999., 2);
0318               CovM3(2, 2) = pow(metRes.phi(pt), 2);
0319             } else {
0320               CovM3(0, 0) = pow(getResolution(object, objType, "et"), 2);
0321               CovM3(1, 1) = pow(getResolution(object, objType, "eta"), 2);
0322               CovM3(2, 2) = pow(getResolution(object, objType, "phi"), 2);
0323             }
0324             return CovM3;
0325           case TopKinFitter::kEtThetaPhi:
0326             CovM3(0, 0) = pow(metRes.et(pt), 2);
0327             CovM3(1, 1) = pow(9999., 2);
0328             CovM3(2, 2) = pow(metRes.phi(pt), 2);
0329             return CovM3;
0330         }
0331       }
0332       break;
0333   }
0334   cms::Exception("Logic") << "Something went wrong while trying to setup a covariance matrix!\n";
0335   return CovM4;  //should never get here
0336 }
0337 
0338 double CovarianceMatrix::getEtaDependentScaleFactor(const TLorentzVector& object) {
0339   double etaDependentScaleFactor = 1.;
0340   for (unsigned int i = 0; i < jetEnergyResolutionEtaBinning_.size(); i++) {
0341     if (std::abs(object.Eta()) >= jetEnergyResolutionEtaBinning_[i] && jetEnergyResolutionEtaBinning_[i] >= 0.) {
0342       if (i == jetEnergyResolutionEtaBinning_.size() - 1) {
0343         edm::LogWarning("CovarianceMatrix") << "object eta (" << std::abs(object.Eta()) << ") beyond last eta bin ("
0344                                             << jetEnergyResolutionEtaBinning_[i] << ") using scale factor 1.0!";
0345         etaDependentScaleFactor = 1.;
0346         break;
0347       }
0348       etaDependentScaleFactor = jetEnergyResolutionScaleFactors_[i];
0349     } else
0350       break;
0351   }
0352   return etaDependentScaleFactor;
0353 }
[ 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!