Back to home page

Project CMSSW displayed by LXR

 
 

    


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 // constructor with path; default should not be used
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   // find etabin values
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 // destructor
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.);  // Total freedom
0177     obj.setResolutionPhi(network_[5]->Evaluate(0, v));
0178     obj.setResolutionEt(network_[6]->Evaluate(0, v));
0179     obj.setResolutionEta(1000000.);  // Total freedom
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.);  // Total freedom
0186     obj.setResolutionPhi(this->obsRes(5, 0, obj.et()));
0187     obj.setResolutionEt(this->obsRes(6, 0, obj.et()));
0188     obj.setResolutionEta(1000000.);  // Total freedom
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