File indexing completed on 2024-04-06 12:24:06
0001
0002
0003
0004 #include "PhysicsTools/PatUtils/interface/ObjectResolutionCalc.h"
0005
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "FWCore/Utilities/interface/Exception.h"
0008
0009 #include "TH1.h"
0010 #include "TKey.h"
0011
0012 using namespace pat;
0013
0014
0015 ObjectResolutionCalc::ObjectResolutionCalc(const TString& resopath, bool useNN = false) : useNN_(useNN) {
0016 edm::LogVerbatim("ObjectResolutionCalc")
0017 << ("ObjectResolutionCalc") << "=== Constructing a TopObjectResolutionCalc...";
0018 resoFile_ = new TFile(resopath);
0019 if (!resoFile_)
0020 throw edm::Exception(edm::errors::LogicError)
0021 << "ObjectResolutionCalc: no resolutions fits for this file available: " << resopath << "...";
0022
0023 TString resObsName[8] = {"_ares", "_bres", "_cres", "_dres", "_thres", "_phres", "_etres", "_etares"};
0024 TList* keys = resoFile_->GetListOfKeys();
0025 TIter nextitem(keys);
0026 TKey* key = nullptr;
0027 while ((key = (TKey*)nextitem())) {
0028 TString name = key->GetName();
0029 if (useNN_) {
0030 for (Int_t ro = 0; ro < 8; ro++) {
0031 TString obsName = resObsName[ro];
0032 obsName += "_NN";
0033 if (name.Contains(obsName)) {
0034 network_[ro] = (TMultiLayerPerceptron*)resoFile_->GetKey(name)->ReadObj();
0035 }
0036 }
0037 } else {
0038 if (name.Contains("etabin") && (!name.Contains("etbin"))) {
0039 for (int p = 0; p < 8; p++) {
0040 if (name.Contains(resObsName[p])) {
0041 TString etabin = name;
0042 etabin.Remove(0, etabin.Index("_") + 1);
0043 etabin.Remove(0, etabin.Index("_") + 7);
0044 int etaBin = etabin.Atoi();
0045 TH1F* tmp = (TH1F*)(resoFile_->GetKey(name)->ReadObj());
0046 fResVsEt_[p][etaBin] = (TF1)(*(tmp->GetFunction("F_" + name)));
0047 }
0048 }
0049 }
0050 }
0051 }
0052
0053 TH1F* tmpEta = (TH1F*)(resoFile_->GetKey("hEtaBins")->ReadObj());
0054 for (int b = 1; b <= tmpEta->GetNbinsX(); b++)
0055 etaBinVals_.push_back(tmpEta->GetXaxis()->GetBinLowEdge(b));
0056 etaBinVals_.push_back(tmpEta->GetXaxis()->GetBinUpEdge(tmpEta->GetNbinsX()));
0057 edm::LogVerbatim("ObjectResolutionCalc") << "Found " << etaBinVals_.size() - 1 << " eta-bins with edges: ( ";
0058 for (size_t u = 0; u < etaBinVals_.size(); u++)
0059 edm::LogVerbatim("ObjectResolutionCalc") << etaBinVals_[u] << ", ";
0060 edm::LogVerbatim("ObjectResolutionCalc") << "\b\b )" << std::endl;
0061
0062 edm::LogVerbatim("ObjectResolutionCalc") << "=== done." << std::endl;
0063 }
0064
0065
0066 ObjectResolutionCalc::~ObjectResolutionCalc() { delete resoFile_; }
0067
0068 float ObjectResolutionCalc::obsRes(int obs, int eta, float eT) {
0069 if (useNN_)
0070 throw edm::Exception(edm::errors::LogicError,
0071 "TopObjectResolutionCalc::obsRes should never be called when using a NN for resolutions.");
0072 float res = fResVsEt_[obs][eta].Eval(eT);
0073 return res;
0074 }
0075
0076 int ObjectResolutionCalc::etaBin(float eta) {
0077 int nrEtaBins = etaBinVals_.size() - 1;
0078 int bin = nrEtaBins - 1;
0079 for (int i = 0; i < nrEtaBins; i++) {
0080 if (fabs(eta) > etaBinVals_[i] && fabs(eta) < etaBinVals_[i + 1]) {
0081 bin = i;
0082 break;
0083 }
0084 }
0085 return bin;
0086 }
0087
0088 #if OBSOLETE
0089 void ObjectResolutionCalc::operator()(Electron& obj) {
0090 if (useNN_) {
0091 double v[2];
0092 v[0] = obj.et();
0093 v[1] = obj.eta();
0094 obj.setResolutionA(network_[0]->Evaluate(0, v));
0095 obj.setResolutionB(network_[1]->Evaluate(0, v));
0096 obj.setResolutionC(network_[2]->Evaluate(0, v));
0097 obj.setResolutionD(network_[3]->Evaluate(0, v));
0098 obj.setResolutionTheta(network_[4]->Evaluate(0, v));
0099 obj.setResolutionPhi(network_[5]->Evaluate(0, v));
0100 obj.setResolutionEt(network_[6]->Evaluate(0, v));
0101 obj.setResolutionEta(network_[7]->Evaluate(0, v));
0102 } else {
0103 int bin = this->etaBin(obj.eta());
0104 obj.setResolutionA(this->obsRes(0, bin, obj.et()));
0105 obj.setResolutionB(this->obsRes(1, bin, obj.et()));
0106 obj.setResolutionC(this->obsRes(2, bin, obj.et()));
0107 obj.setResolutionD(this->obsRes(3, bin, obj.et()));
0108 obj.setResolutionTheta(this->obsRes(4, bin, obj.et()));
0109 obj.setResolutionPhi(this->obsRes(5, bin, obj.et()));
0110 obj.setResolutionEt(this->obsRes(6, bin, obj.et()));
0111 obj.setResolutionEta(this->obsRes(7, bin, obj.et()));
0112 }
0113 }
0114
0115 void ObjectResolutionCalc::operator()(Muon& obj) {
0116 if (useNN_) {
0117 double v[2];
0118 v[0] = obj.et();
0119 v[1] = obj.eta();
0120 obj.setResolutionA(network_[0]->Evaluate(0, v));
0121 obj.setResolutionB(network_[1]->Evaluate(0, v));
0122 obj.setResolutionC(network_[2]->Evaluate(0, v));
0123 obj.setResolutionD(network_[3]->Evaluate(0, v));
0124 obj.setResolutionTheta(network_[4]->Evaluate(0, v));
0125 obj.setResolutionPhi(network_[5]->Evaluate(0, v));
0126 obj.setResolutionEt(network_[6]->Evaluate(0, v));
0127 obj.setResolutionEta(network_[7]->Evaluate(0, v));
0128 } else {
0129 int bin = this->etaBin(obj.eta());
0130 obj.setResolutionA(this->obsRes(0, bin, obj.et()));
0131 obj.setResolutionB(this->obsRes(1, bin, obj.et()));
0132 obj.setResolutionC(this->obsRes(2, bin, obj.et()));
0133 obj.setResolutionD(this->obsRes(3, bin, obj.et()));
0134 obj.setResolutionTheta(this->obsRes(4, bin, obj.et()));
0135 obj.setResolutionPhi(this->obsRes(5, bin, obj.et()));
0136 obj.setResolutionEt(this->obsRes(6, bin, obj.et()));
0137 obj.setResolutionEta(this->obsRes(7, bin, obj.et()));
0138 }
0139 }
0140
0141 void ObjectResolutionCalc::operator()(Jet& obj) {
0142 if (useNN_) {
0143 double v[2];
0144 v[0] = obj.et();
0145 v[1] = obj.eta();
0146 obj.setResolutionA(network_[0]->Evaluate(0, v));
0147 obj.setResolutionB(network_[1]->Evaluate(0, v));
0148 obj.setResolutionC(network_[2]->Evaluate(0, v));
0149 obj.setResolutionD(network_[3]->Evaluate(0, v));
0150 obj.setResolutionTheta(network_[4]->Evaluate(0, v));
0151 obj.setResolutionPhi(network_[5]->Evaluate(0, v));
0152 obj.setResolutionEt(network_[6]->Evaluate(0, v));
0153 obj.setResolutionEta(network_[7]->Evaluate(0, v));
0154 } else {
0155 int bin = this->etaBin(obj.eta());
0156 obj.setResolutionA(this->obsRes(0, bin, obj.et()));
0157 obj.setResolutionB(this->obsRes(1, bin, obj.et()));
0158 obj.setResolutionC(this->obsRes(2, bin, obj.et()));
0159 obj.setResolutionD(this->obsRes(3, bin, obj.et()));
0160 obj.setResolutionTheta(this->obsRes(4, bin, obj.et()));
0161 obj.setResolutionPhi(this->obsRes(5, bin, obj.et()));
0162 obj.setResolutionEt(this->obsRes(6, bin, obj.et()));
0163 obj.setResolutionEta(this->obsRes(7, bin, obj.et()));
0164 }
0165 }
0166
0167 void ObjectResolutionCalc::operator()(MET& obj) {
0168 if (useNN_) {
0169 double v[2];
0170 v[0] = obj.et();
0171 v[1] = obj.eta();
0172 obj.setResolutionA(network_[0]->Evaluate(0, v));
0173 obj.setResolutionB(network_[1]->Evaluate(0, v));
0174 obj.setResolutionC(network_[2]->Evaluate(0, v));
0175 obj.setResolutionD(network_[3]->Evaluate(0, v));
0176 obj.setResolutionTheta(1000000.);
0177 obj.setResolutionPhi(network_[5]->Evaluate(0, v));
0178 obj.setResolutionEt(network_[6]->Evaluate(0, v));
0179 obj.setResolutionEta(1000000.);
0180 } else {
0181 obj.setResolutionA(this->obsRes(0, 0, obj.et()));
0182 obj.setResolutionC(this->obsRes(1, 0, obj.et()));
0183 obj.setResolutionB(this->obsRes(2, 0, obj.et()));
0184 obj.setResolutionD(this->obsRes(3, 0, obj.et()));
0185 obj.setResolutionTheta(1000000.);
0186 obj.setResolutionPhi(this->obsRes(5, 0, obj.et()));
0187 obj.setResolutionEt(this->obsRes(6, 0, obj.et()));
0188 obj.setResolutionEta(1000000.);
0189 }
0190 }
0191
0192 void ObjectResolutionCalc::operator()(Tau& obj) {
0193 if (useNN_) {
0194 double v[2];
0195 v[0] = obj.et();
0196 v[1] = obj.eta();
0197 obj.setResolutionA(network_[0]->Evaluate(0, v));
0198 obj.setResolutionB(network_[1]->Evaluate(0, v));
0199 obj.setResolutionC(network_[2]->Evaluate(0, v));
0200 obj.setResolutionD(network_[3]->Evaluate(0, v));
0201 obj.setResolutionTheta(network_[4]->Evaluate(0, v));
0202 obj.setResolutionPhi(network_[5]->Evaluate(0, v));
0203 obj.setResolutionEt(network_[6]->Evaluate(0, v));
0204 obj.setResolutionEta(network_[7]->Evaluate(0, v));
0205 } else {
0206 int bin = this->etaBin(obj.eta());
0207 obj.setResolutionA(this->obsRes(0, bin, obj.et()));
0208 obj.setResolutionB(this->obsRes(1, bin, obj.et()));
0209 obj.setResolutionC(this->obsRes(2, bin, obj.et()));
0210 obj.setResolutionD(this->obsRes(3, bin, obj.et()));
0211 obj.setResolutionTheta(this->obsRes(4, bin, obj.et()));
0212 obj.setResolutionPhi(this->obsRes(5, bin, obj.et()));
0213 obj.setResolutionEt(this->obsRes(6, bin, obj.et()));
0214 obj.setResolutionEta(this->obsRes(7, bin, obj.et()));
0215 }
0216 }
0217 #endif
