Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-06-20 01:53:36 

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