Phase2L1ParticleFlow/src/corrector.cc

0001 #include "L1Trigger/Phase2L1ParticleFlow/src/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 #include "FWCore/Utilities/interface/CPUTimer.h"
0015 #include "FWCore/Utilities/interface/Exception.h"
0016 #include "FWCore/ParameterSet/interface/FileInPath.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018
0019 #include "DataFormats/L1TParticleFlow/interface/PFCluster.h"
0020
0021 l1tpf::corrector::corrector(const std::string &filename, float emfMax, bool debug) : emfMax_(emfMax) {
0022   if (!filename.empty())
0023     init_(filename, "", debug);
0024 }
0025 l1tpf::corrector::corrector(const std::string &filename, const std::string &directory, float emfMax, bool debug)
0026     : emfMax_(emfMax) {
0027   if (!filename.empty())
0028     init_(filename, directory, debug);
0029 }
0030
0031 l1tpf::corrector::corrector(TDirectory *src, float emfMax, bool debug) : emfMax_(emfMax) { init_(src, debug); }
0032
0033 void l1tpf::corrector::init_(const std::string &filename, const std::string &directory, bool debug) {
0034   std::string resolvedFileName = filename;
0035   if (filename[0] != '/')
0036     resolvedFileName = edm::FileInPath(filename).fullPath();
0037   TFile *lFile = TFile::Open(resolvedFileName.c_str());
0038   if (!lFile || lFile->IsZombie())
0039     throw cms::Exception("Configuration", "cannot read file " + filename);
0040
0041   TDirectory *dir = directory.empty() ? lFile : lFile->GetDirectory(directory.c_str());
0042   if (!dir)
0043     throw cms::Exception("Configuration", "cannot find directory '" + directory + "' in file " + filename);
0044   init_(dir, debug);
0045
0046   lFile->Close();
0047 }
0048
0049 void l1tpf::corrector::init_(TDirectory *lFile, bool debug) {
0050   TH1 *index = (TH1 *)lFile->Get("INDEX");
0051   if (!index)
0052     throw cms::Exception("Configuration")
0053         << "invalid input file " << lFile->GetPath() << ": INDEX histogram not found.\n";
0054   index_.reset((TH1 *)index->Clone());
0055   index_->SetDirectory(nullptr);
0056
0057   is2d_ = index_->InheritsFrom("TH2");
0058
0059   std::unordered_map<std::string, TGraph *> graphs;
0060   TKey *key;
0061   TIter nextkey(lFile->GetListOfKeys());
0062   while ((key = (TKey *)nextkey())) {
0063     if (strncmp(key->GetName(), "eta_", 4) == 0) {
0064       TGraph *gr = (TGraph *)key->ReadObj();
0065       if (!gr->TestBit(TGraph::kIsSortedX))
0066         gr->Sort();
0067       graphs[key->GetName()] = gr;
0068     }
0069   }
0070
0071   neta_ = index_->GetNbinsX();
0072   nemf_ = (is2d_ ? index_->GetNbinsY() : 1);
0073   corrections_.resize(neta_ * nemf_);
0074   std::fill(corrections_.begin(), corrections_.end(), nullptr);
0075   char buff[32];
0076   int ngraphs = 0;
0077   for (unsigned int ieta = 0; ieta < neta_; ++ieta) {
0078     for (unsigned int iemf = 0; iemf < nemf_; ++iemf) {
0079       if (is2d_) {
0080         snprintf(buff, 31, "eta_bin%d_emf_bin%d", ieta + 1, iemf + 1);
0081       } else {
0082         snprintf(buff, 31, "eta_bin%d", ieta + 1);
0083       }
0084       TGraph *graph = graphs[buff];
0085       if (debug)
0086         edm::LogPrint("corrector") << "   eta bin " << ieta << " emf bin " << iemf << " graph " << buff
0087                                    << (graph ? " <valid>" : " <nil>") << "\n";
0088       if (graph) {
0089         ngraphs++;
0090         corrections_[ieta * nemf_ + iemf] = (TGraph *)graph->Clone();
0091       }
0092       if (std::abs(index_->GetXaxis()->GetBinCenter(ieta + 1)) > 3.0) {
0093         break;  // no EMF bins beyond eta = 3
0094       }
0095     }
0096   }
0097 }
0098
0099 l1tpf::corrector::corrector(const TH1 *index, float emfMax)
0100     : index_((TH1 *)index->Clone("INDEX")),
0101       is2d_(index->InheritsFrom("TH2")),
0102       neta_(index->GetNbinsX()),
0103       nemf_(is2d_ ? index->GetNbinsY() : 1),
0104       emfMax_(emfMax) {
0105   index_->SetDirectory(nullptr);
0106   corrections_.resize(neta_ * nemf_);
0107   std::fill(corrections_.begin(), corrections_.end(), nullptr);
0108 }
0109
0110 l1tpf::corrector::~corrector() {
0111   for (TGraph *&p : corrections_) {
0112     delete p;
0113     p = nullptr;
0114   }
0115   corrections_.clear();
0116 }
0117
0118 l1tpf::corrector::corrector(corrector &&corr)
0119     : index_(std::move(corr.index_)),
0120       corrections_(std::move(corr.corrections_)),
0121       is2d_(corr.is2d_),
0122       neta_(corr.neta_),
0123       nemf_(corr.nemf_),
0124       emfMax_(corr.emfMax_) {}
0125
0126 l1tpf::corrector &l1tpf::corrector::operator=(corrector &&corr) {
0127   std::swap(is2d_, corr.is2d_);
0128   std::swap(neta_, corr.neta_);
0129   std::swap(nemf_, corr.nemf_);
0130   std::swap(emfMax_, corr.emfMax_);
0131
0132   index_.swap(corr.index_);
0133   corrections_.swap(corr.corrections_);
0134   return *this;
0135 }
0136
0137 float l1tpf::corrector::correctedPt(float pt, float emPt, float eta) const {
0138   float total = std::max(pt, emPt), abseta = std::abs(eta);
0139   float emf = emPt / total;
0140   if (emfMax_ > 0 && emf > emfMax_)
0141     return total;  // no correction
0142   unsigned int ieta = std::min(std::max<unsigned>(1, index_->GetXaxis()->FindBin(abseta)), neta_) - 1;
0143   unsigned int iemf =
0144       is2d_ && abseta < 3.0 ? std::min(std::max<unsigned>(1, index_->GetYaxis()->FindBin(emf)), nemf_) - 1 : 0;
0145   TGraph *graph = corrections_[ieta * nemf_ + iemf];
0146   if (!graph) {
0147     throw cms::Exception("RuntimeError") << "Error trying to read calibration for eta " << eta << " emf " << emf
0148                                          << " which are not available." << std::endl;
0149   }
0150   float ptcorr = std::min<float>(graph->Eval(total), 4 * total);
0151   return ptcorr;
0152 }
0153
0154 void l1tpf::corrector::correctPt(l1t::PFCluster &cluster, float preserveEmEt) const {
0155   float newpt = correctedPt(cluster.pt(), cluster.emEt(), cluster.eta());
0156   cluster.calibratePt(newpt, preserveEmEt);
0157 }
0158
0159 void l1tpf::corrector::setGraph(const TGraph &graph, int ieta, int iemf) {
0160   char buff[32];
0161   if (is2d_) {
0162     snprintf(buff, 31, "eta_bin%d_emf_bin%d", (unsigned int)(ieta + 1), (unsigned int)(iemf + 1));
0163   } else {
0164     snprintf(buff, 31, "eta_bin%d", (unsigned int)(ieta + 1));
0165   }
0166   TGraph *gclone = (TGraph *)graph.Clone(buff);
0167   delete corrections_[ieta * nemf_ + iemf];
0168   corrections_[ieta * nemf_ + iemf] = gclone;
0169 }
0170
0171 void l1tpf::corrector::writeToFile(const std::string &filename, const std::string &directory) const {
0172   TFile *lFile = TFile::Open(filename.c_str(), "RECREATE");
0173   TDirectory *dir = directory.empty() ? lFile : lFile->mkdir(directory.c_str());
0174   writeToFile(dir);
0175   lFile->Close();
0176 }
0177
0178 void l1tpf::corrector::writeToFile(TDirectory *dest) const {
0179   TH1 *index = (TH1 *)index_->Clone();
0180   index->SetDirectory(dest);
0181   dest->WriteTObject(index);
0182
0183   for (const TGraph *p : corrections_) {
0184     if (p != nullptr) {
0185       dest->WriteTObject((TGraph *)p->Clone(), p->GetName());
0186     }
0187   }
0188 }