File indexing completed on 2024-05-22 04:02:32
0001 #ifndef Alignment_OfflineValidation_DiLeptonVertexHelpers_h
0002 #define Alignment_OfflineValidation_DiLeptonVertexHelpers_h
0003
0004 #include <vector>
0005 #include <string>
0006 #include <fmt/printf.h>
0007 #include "TH2F.h"
0008 #include "TLorentzVector.h"
0009 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012
0013 namespace DiLeptonHelp {
0014
0015
0016
0017
0018 struct Counts {
0019 unsigned int eventsTotal;
0020 unsigned int eventsAfterMult;
0021 unsigned int eventsAfterPt;
0022 unsigned int eventsAfterEta;
0023 unsigned int eventsAfterVtx;
0024 unsigned int eventsAfterDist;
0025
0026 public:
0027 void printCounts() {
0028 edm::LogInfo("DiLeptonHelpCounts") << " Total Events: " << eventsTotal << "\n"
0029 << " After multiplicity: " << eventsAfterMult << "\n"
0030 << " After pT cut: " << eventsAfterPt << "\n"
0031 << " After eta cut: " << eventsAfterEta << "\n"
0032 << " After Vtx: " << eventsAfterVtx << "\n"
0033 << " After VtxDist: " << eventsAfterDist << std::endl;
0034 }
0035
0036 void zeroAll() {
0037 eventsTotal = 0;
0038 eventsAfterMult = 0;
0039 eventsAfterPt = 0;
0040 eventsAfterEta = 0;
0041 eventsAfterVtx = 0;
0042 eventsAfterDist = 0;
0043 }
0044 };
0045
0046 enum flavour { MM = 0, EE = 1, UNDEF = -1 };
0047
0048 enum etaRegion { BARBAR, BARFWD, BARBWD, FWDFWD, BWDBWD, FWDBWD, END };
0049
0050
0051
0052
0053
0054 class PlotsVsDiLeptonRegion {
0055 public:
0056 PlotsVsDiLeptonRegion(const float etaBoundary) : m_etaBoundary(etaBoundary) {}
0057 ~PlotsVsDiLeptonRegion() = default;
0058
0059
0060 inline void bookSet(const TFileDirectory& fs, const TH1* histo) {
0061 const std::string name = histo->GetName();
0062 const std::string title = histo->GetTitle();
0063 const std::string xTitle = histo->GetXaxis()->GetTitle();
0064 const std::string yTitle = histo->GetYaxis()->GetTitle();
0065 std::string zTitle = "";
0066 if (((TObject*)histo)->InheritsFrom("TH2")) {
0067 zTitle = histo->GetZaxis()->GetTitle();
0068 }
0069
0070 for (const auto& etaReg : m_etaRegions) {
0071 if (etaReg == etaRegion::END)
0072 continue;
0073
0074 const auto& toSub = m_etaRegionNames[etaReg];
0075
0076 if (((TObject*)histo)->InheritsFrom("TH2")) {
0077 m_h2_map[etaReg] = fs.make<TH2F>(
0078 (name + "_" + toSub).c_str(),
0079 (title + " (" + toSub + ");" + xTitle + " (" + toSub + ") ;" + yTitle + " (" + toSub + ");" + zTitle)
0080 .c_str(),
0081 histo->GetNbinsX(),
0082 histo->GetXaxis()->GetXmin(),
0083 histo->GetXaxis()->GetXmax(),
0084 histo->GetNbinsY(),
0085 histo->GetYaxis()->GetXmin(),
0086 histo->GetYaxis()->GetXmax());
0087 } else {
0088 m_h1_map[etaReg] = fs.make<TH1F>(
0089 (name + "_" + toSub).c_str(),
0090 (title + " (" + toSub + ");" + xTitle + " (" + toSub + ") ;" + yTitle + " (" + toSub + ")").c_str(),
0091 histo->GetNbinsX(),
0092 histo->GetXaxis()->GetXmin(),
0093 histo->GetXaxis()->GetXmax());
0094 }
0095 }
0096
0097
0098 m_isBooked = true;
0099 }
0100
0101
0102
0103 etaRegion getEtaRegion(const double eta1, const double eta2) {
0104 bool isEta1Barrel = std::abs(eta1) <= m_etaBoundary;
0105 bool isEta2Barrel = std::abs(eta2) <= m_etaBoundary;
0106
0107 if (isEta1Barrel && isEta2Barrel) {
0108 return etaRegion::BARBAR;
0109 } else if ((isEta1Barrel && eta2 > m_etaBoundary) || (isEta2Barrel && eta1 > m_etaBoundary)) {
0110 return etaRegion::BARFWD;
0111 } else if ((isEta1Barrel && eta2 < -m_etaBoundary) || (isEta2Barrel && eta1 < -m_etaBoundary)) {
0112 return etaRegion::BARBWD;
0113 } else if (eta1 > m_etaBoundary && eta2 > m_etaBoundary) {
0114 return etaRegion::FWDFWD;
0115 } else if (eta1 < -m_etaBoundary && eta2 < -m_etaBoundary) {
0116 return etaRegion::BWDBWD;
0117 } else if ((eta1 > m_etaBoundary && eta2 < -m_etaBoundary) || (eta2 > m_etaBoundary && eta1 < -m_etaBoundary)) {
0118 return etaRegion::FWDBWD;
0119 }
0120
0121
0122 return etaRegion::END;
0123 }
0124
0125
0126 inline void fillTH1Plots(const float val, const std::pair<TLorentzVector, TLorentzVector>& momenta) {
0127 if (!m_isBooked) {
0128 edm::LogError("PlotsVsDiLeptonRegion")
0129 << "In" << __FUNCTION__ << "," << __LINE__ << "trying to fill a plot not booked!" << std::endl;
0130 return;
0131 }
0132
0133 etaRegion region = getEtaRegion(momenta.first.Eta(), momenta.second.Eta());
0134 if (region == etaRegion::END) {
0135 edm::LogError("PlotsVsDiLeptonRegion") << "undefined di-muon kinematics" << std::endl;
0136 }
0137 m_h1_map[region]->Fill(val);
0138 }
0139
0140
0141 inline void fillTH2Plots(const float valX,
0142 const float valY,
0143 const std::pair<TLorentzVector, TLorentzVector>& momenta) {
0144 if (!m_isBooked) {
0145 edm::LogError("PlotsVsDiLeptonRegion")
0146 << "In" << __FUNCTION__ << "," << __LINE__ << "trying to fill a plot not booked!" << std::endl;
0147 return;
0148 }
0149
0150 etaRegion region = getEtaRegion(momenta.first.Eta(), momenta.second.Eta());
0151 if (region == etaRegion::END) {
0152 edm::LogError("PlotsVsDiLeptonRegion") << "undefined di-muon kinematics" << std::endl;
0153 }
0154 m_h2_map[region]->Fill(valX, valY);
0155 }
0156
0157 private:
0158 const std::vector<etaRegion> m_etaRegions = {etaRegion::BARBAR,
0159 etaRegion::BARFWD,
0160 etaRegion::BARBWD,
0161 etaRegion::FWDFWD,
0162 etaRegion::BWDBWD,
0163 etaRegion::FWDBWD};
0164
0165 const std::vector<std::string> m_etaRegionNames = {"barrel-barrel",
0166 "barrel-forward",
0167 "barrel-backward",
0168 "forward-forward",
0169 "backward-backward",
0170 "forward-backward"};
0171 const float m_etaBoundary;
0172 bool m_isBooked;
0173 std::map<etaRegion, TH1F*> m_h1_map;
0174 std::map<etaRegion, TH2F*> m_h2_map;
0175 };
0176
0177
0178
0179
0180 class PlotsVsKinematics {
0181 public:
0182 PlotsVsKinematics(flavour FLAV) : m_name(""), m_title(""), m_ytitle(""), m_isBooked(false), m_flav(FLAV) {}
0183
0184
0185
0186 PlotsVsKinematics(flavour FLAV, const std::string& name, const std::string& tt, const std::string& ytt)
0187 : m_name(name), m_title(tt), m_ytitle(ytt), m_isBooked(false), m_flav(FLAV) {
0188 if (m_flav < 0) {
0189 edm::LogError("PlotsVsKinematics") << "The initialization flavour is not correct!" << std::endl;
0190 }
0191 }
0192
0193 ~PlotsVsKinematics() = default;
0194
0195
0196 inline void bookFromPSet(const TFileDirectory& fs, const edm::ParameterSet& hpar) {
0197 std::string namePostfix;
0198 std::string titlePostfix;
0199 float xmin, xmax;
0200
0201 std::string sed = (m_flav ? "e" : "#mu");
0202
0203 for (const auto& xAx : axisChoices) {
0204 switch (xAx) {
0205 case xAxis::Z_PHI:
0206 xmin = -M_PI;
0207 xmax = M_PI;
0208 namePostfix = m_flav ? "EEPhi" : "MMPhi";
0209 titlePostfix = fmt::sprintf("%s%s pair #phi;%s^{+}%s^{-} #phi", sed, sed, sed, sed);
0210 break;
0211 case xAxis::Z_ETA:
0212 xmin = -3.5;
0213 xmax = 3.5;
0214 namePostfix = m_flav ? "EEEta" : "MuMuEta";
0215 titlePostfix = fmt::sprintf("%s%s pair #eta;%s^{+}%s^{-} #eta", sed, sed, sed, sed);
0216 break;
0217 case xAxis::LP_PHI:
0218 xmin = -M_PI;
0219 xmax = M_PI;
0220 namePostfix = m_flav ? "EPlusPhi" : "MuPlusPhi";
0221 titlePostfix = fmt::sprintf("%s^{+} #phi;%s^{+} #phi [rad]", sed, sed);
0222 break;
0223 case xAxis::LP_ETA:
0224 xmin = -2.4;
0225 xmax = 2.4;
0226 namePostfix = m_flav ? "EPlusEta" : "MuPlusEta";
0227 titlePostfix = fmt::sprintf("%s^{+} #eta;%s^{+} #eta", sed, sed);
0228 break;
0229 case xAxis::LM_PHI:
0230 xmin = -M_PI;
0231 xmax = M_PI;
0232 namePostfix = m_flav ? "EMinusPhi" : "MuMinusPhi";
0233 titlePostfix = fmt::sprintf("%s^{-} #phi;%s^{-} #phi [rad]", sed, sed);
0234 break;
0235 case xAxis::LM_ETA:
0236 xmin = -2.4;
0237 xmax = 2.4;
0238 namePostfix = m_flav ? "EMinusEta" : "MuMinusEta";
0239 titlePostfix = fmt::sprintf("%s^{-} #eta;%s^{+} #eta", sed, sed);
0240 break;
0241 default:
0242 throw cms::Exception("LogicalError") << " there is not such Axis choice as " << xAx;
0243 }
0244
0245 const auto& h2name = fmt::sprintf("%sVs%s", hpar.getParameter<std::string>("name"), namePostfix);
0246 const auto& h2title = fmt::sprintf("%s vs %s;%s% s",
0247 hpar.getParameter<std::string>("title"),
0248 titlePostfix,
0249 hpar.getParameter<std::string>("title"),
0250 hpar.getParameter<std::string>("yUnits"));
0251
0252 m_h2_map[xAx] = fs.make<TH2F>(h2name.c_str(),
0253 h2title.c_str(),
0254 hpar.getParameter<int32_t>("NxBins"),
0255 xmin,
0256 xmax,
0257 hpar.getParameter<int32_t>("NyBins"),
0258 hpar.getParameter<double>("ymin"),
0259 hpar.getParameter<double>("ymax"));
0260 }
0261
0262
0263 m_isBooked = true;
0264 }
0265
0266
0267 inline void bookPlots(
0268 TFileDirectory& fs, const float valmin, const float valmax, const int nxbins, const int nybins) {
0269 if (m_name.empty() && m_title.empty() && m_ytitle.empty()) {
0270 edm::LogError("PlotsVsKinematics")
0271 << "In" << __FUNCTION__ << "," << __LINE__
0272 << "trying to book plots without the right constructor being called!" << std::endl;
0273 return;
0274 }
0275
0276 std::string dilep = (m_flav ? "e^{+}e^{-}" : "#mu^{+}#mu^{-}");
0277 std::string lep = (m_flav ? "e^{+}" : "#mu^{+}");
0278 std::string lem = (m_flav ? "e^{-}" : "#mu^{-}");
0279
0280 static constexpr float maxMuEta = 2.4;
0281 static constexpr float maxMuMuEta = 3.5;
0282 TH1F::SetDefaultSumw2(kTRUE);
0283
0284
0285 m_h2_map[xAxis::Z_ETA] = fs.make<TH2F>(fmt::sprintf("%sVsMuMuEta", m_name).c_str(),
0286 fmt::sprintf("%s vs %s pair #eta;%s #eta;%s", m_title, dilep, dilep, m_ytitle).c_str(),
0287 nxbins, -M_PI, M_PI,
0288 nybins, valmin, valmax);
0289
0290 m_h2_map[xAxis::Z_PHI] = fs.make<TH2F>(fmt::sprintf("%sVsMuMuPhi", m_name).c_str(),
0291 fmt::sprintf("%s vs %s pair #phi;%s #phi [rad];%s", m_title, dilep, dilep, m_ytitle).c_str(),
0292 nxbins, -maxMuMuEta, maxMuMuEta,
0293 nybins, valmin, valmax);
0294
0295 m_h2_map[xAxis::LP_ETA] = fs.make<TH2F>(fmt::sprintf("%sVsMuPlusEta", m_name).c_str(),
0296 fmt::sprintf("%s vs %s #eta;%s #eta;%s", m_title, lep, lep, m_ytitle).c_str(),
0297 nxbins, -maxMuEta, maxMuEta,
0298 nybins, valmin, valmax);
0299
0300 m_h2_map[xAxis::LP_PHI] = fs.make<TH2F>(fmt::sprintf("%sVsMuPlusPhi", m_name).c_str(),
0301 fmt::sprintf("%s vs %s #phi;%s #phi [rad];%s", m_title, lep, lep, m_ytitle).c_str(),
0302 nxbins, -M_PI, M_PI,
0303 nybins, valmin, valmax);
0304
0305 m_h2_map[xAxis::LM_ETA] = fs.make<TH2F>(fmt::sprintf("%sVsMuMinusEta", m_name).c_str(),
0306 fmt::sprintf("%s vs %s #eta;%s #eta;%s", m_title, lem, lem, m_ytitle).c_str(),
0307 nxbins, -maxMuEta, maxMuEta,
0308 nybins, valmin, valmax);
0309
0310 m_h2_map[xAxis::LM_PHI] = fs.make<TH2F>(fmt::sprintf("%sVsMuMinusPhi", m_name).c_str(),
0311 fmt::sprintf("%s vs %s #phi;%s #phi [rad];%s", m_title, lem, lem, m_ytitle).c_str(),
0312 nxbins, -M_PI, M_PI,
0313 nybins, valmin, valmax);
0314
0315
0316
0317 m_isBooked = true;
0318 }
0319
0320
0321 inline void fillPlots(const float val, const std::pair<TLorentzVector, TLorentzVector>& momenta) {
0322 if (!m_isBooked) {
0323 edm::LogError("PlotsVsKinematics")
0324 << "In" << __FUNCTION__ << "," << __LINE__ << "trying to fill a plot not booked!" << std::endl;
0325 return;
0326 }
0327
0328 m_h2_map[xAxis::Z_ETA]->Fill((momenta.first + momenta.second).Eta(), val);
0329 m_h2_map[xAxis::Z_PHI]->Fill((momenta.first + momenta.second).Phi(), val);
0330 m_h2_map[xAxis::LP_ETA]->Fill((momenta.first).Eta(), val);
0331 m_h2_map[xAxis::LP_PHI]->Fill((momenta.first).Phi(), val);
0332 m_h2_map[xAxis::LM_ETA]->Fill((momenta.second).Eta(), val);
0333 m_h2_map[xAxis::LM_PHI]->Fill((momenta.second).Phi(), val);
0334 }
0335
0336 private:
0337 enum xAxis { Z_PHI, Z_ETA, LP_PHI, LP_ETA, LM_PHI, LM_ETA };
0338 const std::vector<xAxis> axisChoices = {
0339 xAxis::Z_PHI, xAxis::Z_ETA, xAxis::LP_PHI, xAxis::LP_ETA, xAxis::LM_PHI, xAxis::LM_ETA};
0340
0341 const std::string m_name;
0342 const std::string m_title;
0343 const std::string m_ytitle;
0344
0345 bool m_isBooked;
0346 flavour m_flav;
0347
0348 std::map<xAxis, TH2F*> m_h2_map;
0349 };
0350 }
0351 #endif
