Phase2L1ParticleFlow/src/corrector.cc

0001 #include "L1Trigger/Phase2L1ParticleFlow/interface/corrector.h"
0002
0003 #include <iostream>
0004 #include <sstream>
0005 #include <cstdio>
0006 #include <cstdlib>
0007 #include <cassert>
0008 #include <unordered_map>
0009 #include <TFile.h>
0010 #include <TKey.h>
0011 #include <TH1.h>
0012 #include <TH2.h>
0013 #include <TAxis.h>
0014
0015 #ifdef CMSSW_GIT_HASH
0016 #include "FWCore/Utilities/interface/CPUTimer.h"
0017 #include "FWCore/Utilities/interface/Exception.h"
0018 #include "FWCore/ParameterSet/interface/FileInPath.h"
0019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0020
0021 #include "DataFormats/L1TParticleFlow/interface/PFCluster.h"
0022 #else
0023 #include <filesystem>
0024 #include <sstream>
0025 #include <stdexcept>
0026 #endif
0027
0028 /* ---
0029  * Note: #ifdef CMSSW_GIT_HASH is used to determine whether compilation is in CMSSW context or not
0030  * since this same implementation is used in the standalone tests of corrections in comparison to
0031  * L1T firmware
0032  * --- */
0033
0034 l1tpf::corrector::corrector(const std::string &filename, float emfMax, bool debug, bool emulate)
0035     : emfMax_(emfMax), emulate_(emulate) {
0036   if (!filename.empty())
0037     init_(filename, "", debug, emulate);
0038 }
0039
0040 l1tpf::corrector::corrector(
0041     const std::string &filename, const std::string &directory, float emfMax, bool debug, bool emulate)
0042     : emfMax_(emfMax), emulate_(emulate) {
0043   if (!filename.empty())
0044     init_(filename, directory, debug, emulate);
0045 }
0046
0047 l1tpf::corrector::corrector(TDirectory *src, float emfMax, bool debug, bool emulate)
0048     : emfMax_(emfMax), emulate_(emulate) {
0049   init_(src, debug);
0050 }
0051
0052 void l1tpf::corrector::init_(const std::string &filename, const std::string &directory, bool debug, bool emulate) {
0053   std::string resolvedFileName = filename;
0054   if (filename[0] != '/') {
0055 #ifdef CMSSW_GIT_HASH
0056     resolvedFileName = edm::FileInPath(filename).fullPath();
0057 #else
0058     resolvedFileName = std::filesystem::absolute(filename);
0059 #endif
0060   }
0061   TFile *lFile = TFile::Open(resolvedFileName.c_str());
0062   if (!lFile || lFile->IsZombie()) {
0063 #ifdef CMSSW_GIT_HASH
0064     throw cms::Exception("Configuration", "cannot read file " + filename);
0065 #else
0066     throw std::runtime_error("Cannot read file " + filename);
0067 #endif
0068   }
0069
0070   TDirectory *dir = directory.empty() ? lFile : lFile->GetDirectory(directory.c_str());
0071   if (!dir) {
0072 #ifdef CMSSW_GIT_HASH
0073     throw cms::Exception("Configuration", "cannot find directory '" + directory + "' in file " + filename);
0074 #else
0075     throw std::runtime_error("Cannot find directory '" + directory + "' in file " + filename);
0076 #endif
0077   }
0078   init_(dir, debug);
0079   if (emulate)
0080     initEmulation_(dir, debug);
0081
0082   lFile->Close();
0083 }
0084
0085 void l1tpf::corrector::init_(TDirectory *lFile, bool debug) {
0086   TH1 *index = (TH1 *)lFile->Get("INDEX");
0087   if (!index) {
0088 #ifdef CMSSW_GIT_HASH
0089     throw cms::Exception("Configuration")
0090         << "invalid input file " << lFile->GetPath() << ": INDEX histogram not found.\n";
0091 #else
0092     std::stringstream ss;
0093     ss << "invalid input file " << lFile->GetPath() << ": INDEX histogram nit found.\n";
0094     throw std::runtime_error(ss.str());
0095 #endif
0096   }
0097   index_.reset((TH1 *)index->Clone());
0098   index_->SetDirectory(nullptr);
0099
0100   is2d_ = index_->InheritsFrom("TH2");
0101
0102   std::unordered_map<std::string, TGraph *> graphs;
0103   TKey *key;
0104   TIter nextkey(lFile->GetListOfKeys());
0105   while ((key = (TKey *)nextkey())) {
0106     if (strncmp(key->GetName(), "eta_", 4) == 0) {
0107       TGraph *gr = (TGraph *)key->ReadObj();
0108       if (!gr->TestBit(TGraph::kIsSortedX))
0109         gr->Sort();
0110       graphs[key->GetName()] = gr;
0111     }
0112   }
0113
0114   neta_ = index_->GetNbinsX();
0115   nemf_ = (is2d_ ? index_->GetNbinsY() : 1);
0116   corrections_.resize(neta_ * nemf_);
0117   std::fill(corrections_.begin(), corrections_.end(), nullptr);
0118   char buff[32];
0119   for (unsigned int ieta = 0; ieta < neta_; ++ieta) {
0120     for (unsigned int iemf = 0; iemf < nemf_; ++iemf) {
0121       if (is2d_) {
0122         snprintf(buff, 31, "eta_bin%d_emf_bin%d", ieta + 1, iemf + 1);
0123       } else {
0124         snprintf(buff, 31, "eta_bin%d", ieta + 1);
0125       }
0126       TGraph *graph = graphs[buff];
0127       if (debug) {
0128 #ifdef CMSSW_GIT_HASH
0129         edm::LogPrint("corrector") << "   eta bin " << ieta << " emf bin " << iemf << " graph " << buff
0130                                    << (graph ? " <valid>" : " <nil>") << "\n";
0131 #else
0132         std::cout << "   eta bin " << ieta << " emf bin " << iemf << " graph " << buff
0133                   << (graph ? " <valid>" : " <nil>") << "\n";
0134 #endif
0135       }
0136       if (graph) {
0137         corrections_[ieta * nemf_ + iemf] = (TGraph *)graph->Clone();
0138       }
0139       if (std::abs(index_->GetXaxis()->GetBinCenter(ieta + 1)) > 3.0) {
0140         break;  // no EMF bins beyond eta = 3
0141       }
0142     }
0143   }
0144 }
0145
0146 void l1tpf::corrector::initEmulation_(TDirectory *lFile, bool debug) {
0147   std::unordered_map<std::string, TH1 *> hists;
0148   TKey *key;
0149   TIter nextkey(lFile->GetListOfKeys());
0150   while ((key = (TKey *)nextkey())) {
0151     if (strncmp(key->GetName(), "emul_eta", 8) == 0) {
0152       TH1 *hist = (TH1 *)key->ReadObj();
0153       hists[key->GetName()] = hist;
0154     }
0155   }
0156
0157   neta_ = index_->GetNbinsX();
0158   correctionsEmulated_.resize(neta_);
0159   std::fill(correctionsEmulated_.begin(), correctionsEmulated_.end(), nullptr);
0160   char buff[32];
0161   for (unsigned int ieta = 0; ieta < neta_; ++ieta) {
0162     snprintf(buff, 31, "emul_eta_bin%d", ieta + 1);
0163     TH1 *hist = hists[buff];
0164     if (debug)
0165 #ifdef CMSSW_GIT_HASH
0166       edm::LogPrint("corrector") << "   eta bin " << ieta << " hist " << buff << (hist ? " <valid>" : " <nil>") << "\n";
0167 #else
0168       std::cout << "   eta bin " << ieta << " hist " << buff << (hist ? " <valid>" : " <nil>") << "\n";
0169 #endif
0170     if (hist) {
0171       correctionsEmulated_[ieta] = (TH1 *)hist->Clone();
0172       correctionsEmulated_[ieta]->SetDirectory(nullptr);
0173     }
0174   }
0175 }
0176
0177 l1tpf::corrector::corrector(const TH1 *index, float emfMax)
0178     : index_((TH1 *)index->Clone("INDEX")),
0179       is2d_(index->InheritsFrom("TH2")),
0180       neta_(index->GetNbinsX()),
0181       nemf_(is2d_ ? index->GetNbinsY() : 1),
0182       emfMax_(emfMax) {
0183   index_->SetDirectory(nullptr);
0184   corrections_.resize(neta_ * nemf_);
0185   std::fill(corrections_.begin(), corrections_.end(), nullptr);
0186 }
0187
0188 l1tpf::corrector::~corrector() {
0189   for (TGraph *&p : corrections_) {
0190     delete p;
0191     p = nullptr;
0192   }
0193   corrections_.clear();
0194   for (TH1 *&p : correctionsEmulated_) {
0195     delete p;
0196     p = nullptr;
0197   }
0198   correctionsEmulated_.clear();
0199 }
0200
0201 l1tpf::corrector::corrector(corrector &&corr)
0202     : index_(std::move(corr.index_)),
0203       corrections_(std::move(corr.corrections_)),
0204       correctionsEmulated_(std::move(corr.correctionsEmulated_)),
0205       is2d_(corr.is2d_),
0206       neta_(corr.neta_),
0207       nemf_(corr.nemf_),
0208       emfMax_(corr.emfMax_),
0209       emulate_(corr.emulate_) {}
0210
0211 l1tpf::corrector &l1tpf::corrector::operator=(corrector &&corr) {
0212   std::swap(is2d_, corr.is2d_);
0213   std::swap(neta_, corr.neta_);
0214   std::swap(nemf_, corr.nemf_);
0215   std::swap(emfMax_, corr.emfMax_);
0216   std::swap(emulate_, corr.emulate_);
0217
0218   index_.swap(corr.index_);
0219   corrections_.swap(corr.corrections_);
0220   correctionsEmulated_.swap(corr.correctionsEmulated_);
0221   return *this;
0222 }
0223
0224 float l1tpf::corrector::correctedPt(float pt, float emPt, float eta) const {
0225   unsigned int ipt, ieta;
0226   float total = std::max(pt, emPt), abseta = std::abs(eta);
0227   float emf = emPt / total;
0228   if (emfMax_ > 0 && emf > emfMax_)
0229     return total;  // no correction
0230   ieta = std::min(std::max<unsigned>(1, index_->GetXaxis()->FindBin(abseta)), neta_) - 1;
0231   unsigned int iemf =
0232       is2d_ && abseta < 3.0 ? std::min(std::max<unsigned>(1, index_->GetYaxis()->FindBin(emf)), nemf_) - 1 : 0;
0233   float ptcorr = 0;
0234   if (!emulate_) {  // not emulation - read from the TGraph as normal
0235     TGraph *graph = corrections_[ieta * nemf_ + iemf];
0236     if (!graph) {
0237 #ifdef CMSSW_GIT_HASH
0238       throw cms::Exception("RuntimeError") << "Error trying to read calibration for eta " << eta << " emf " << emf
0239                                            << " which are not available." << std::endl;
0240 #else
0241       std::stringstream ss;
0242       ss << "Error trying to read calibration for eta " << eta << " emf " << emf << " which are not available."
0243          << std::endl;
0244       throw std::runtime_error(ss.str());
0245 #endif
0246     }
0247     ptcorr = std::min<float>(graph->Eval(total), 4 * total);
0248   } else {  // emulation - read from the pt binned histogram
0249     TH1 *hist = correctionsEmulated_[ieta];
0250     if (!hist) {
0251 #ifdef CMSSW_GIT_HASH
0252       throw cms::Exception("RuntimeError")
0253           << "Error trying to read emulated calibration for eta " << eta << " which are not available." << std::endl;
0254 #else
0255       std::stringstream ss;
0256       ss << "Error trying to read emulated calibration for eta " << eta << " which are not available." << std::endl;
0257       throw std::runtime_error(ss.str());
0258 #endif
0259     }
0260     ipt = hist->GetXaxis()->FindBin(pt);
0261     ptcorr = hist->GetBinContent(ipt);
0262   }
0263   return ptcorr;
0264 }
0265
0266 #ifdef CMSSW_GIT_HASH
0267 void l1tpf::corrector::correctPt(l1t::PFCluster &cluster, float preserveEmEt) const {
0268   float newpt = correctedPt(cluster.pt(), cluster.emEt(), cluster.eta());
0269   cluster.calibratePt(newpt, preserveEmEt);
0270 }
0271 #endif
0272
0273 void l1tpf::corrector::setGraph(const TGraph &graph, int ieta, int iemf) {
0274   char buff[32];
0275   if (is2d_) {
0276     snprintf(buff, 31, "eta_bin%d_emf_bin%d", (unsigned int)(ieta + 1), (unsigned int)(iemf + 1));
0277   } else {
0278     snprintf(buff, 31, "eta_bin%d", (unsigned int)(ieta + 1));
0279   }
0280   TGraph *gclone = (TGraph *)graph.Clone(buff);
0281   delete corrections_[ieta * nemf_ + iemf];
0282   corrections_[ieta * nemf_ + iemf] = gclone;
0283 }
0284
0285 void l1tpf::corrector::writeToFile(const std::string &filename, const std::string &directory) const {
0286   TFile *lFile = TFile::Open(filename.c_str(), "RECREATE");
0287   TDirectory *dir = directory.empty() ? lFile : lFile->mkdir(directory.c_str());
0288   writeToFile(dir);
0289   lFile->Close();
0290 }
0291
0292 void l1tpf::corrector::writeToFile(TDirectory *dest) const {
0293   TH1 *index = (TH1 *)index_->Clone();
0294   index->SetDirectory(dest);
0295   dest->WriteTObject(index);
0296
0297   for (const TGraph *p : corrections_) {
0298     if (p != nullptr) {
0299       dest->WriteTObject((TGraph *)p->Clone(), p->GetName());
0300     }
0301   }
0302
0303   for (const TH1 *p : correctionsEmulated_) {
0304     if (p != nullptr) {
0305       dest->WriteTObject((TH1 *)p->Clone(), p->GetName());
0306     }
0307   }
0308 }