From 4202a56ccbd00668a19e5df723ec5c3517f3a079 Mon Sep 17 00:00:00 2001 From: mmusich Date: Fri, 17 May 2024 16:16:32 +0200 Subject: [PATCH] add methods to allow auto-scaling of the mass range in the profiles (off by default) --- .../macros/DiMuonMassProfiles.C | 163 +++++++++++++----- 1 file changed, 120 insertions(+), 43 deletions(-) diff --git a/Alignment/OfflineValidation/macros/DiMuonMassProfiles.C b/Alignment/OfflineValidation/macros/DiMuonMassProfiles.C index 38ea202c16167..2601c64099e45 100644 --- a/Alignment/OfflineValidation/macros/DiMuonMassProfiles.C +++ b/Alignment/OfflineValidation/macros/DiMuonMassProfiles.C @@ -29,10 +29,10 @@ #include #include -// style +// CMSSW classes / style +#include "Alignment/OfflineValidation/interface/FitWithRooFit.h" #include "Alignment/OfflineValidation/macros/CMS_lumi.h" #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h" -#include "Alignment/OfflineValidation/interface/FitWithRooFit.h" bool debugMode_{false}; Int_t def_markers[9] = {kFullSquare, @@ -47,6 +47,90 @@ Int_t def_markers[9] = {kFullSquare, Int_t def_colors[9] = {kBlack, kRed, kBlue, kMagenta, kGreen, kCyan, kViolet, kOrange, kGreen + 2}; +namespace diMuonMassBias { + + struct fitOutputs { + public: + fitOutputs(const Measurement1D& bias, const Measurement1D& width) : m_bias(bias), m_width(width) {} + + // getters + const Measurement1D getBias() { return m_bias; } + const Measurement1D getWidth() { return m_width; } + bool isInvalid() const { + return (m_bias.value() == 0.f && m_bias.error() == 0.f && m_width.value() == 0.f && m_width.error() == 0.f); + } + + private: + Measurement1D m_bias; + Measurement1D m_width; + }; + + static constexpr int minimumHits = 10; + + using histoMap = std::map; + using histo2DMap = std::map; + using histo3DMap = std::map; + + using extrema = std::pair; + using extremaMap = std::map; + +} // namespace diMuonMassBias + +/*--------------------------------------------------------------------*/ +template +std::map> transformMaps(const std::vector>& vecOfMaps) +/*--------------------------------------------------------------------*/ +{ + std::map> result; + + // Lambda to insert each key-value pair into the result map + auto insert_into_result = [&result](const std::map& map) { + for (const auto& pair : map) { + result[pair.first].push_back(pair.second); + } + }; + + // Apply the lambda to each map in the vector + std::for_each(vecOfMaps.begin(), vecOfMaps.end(), insert_into_result); + + return result; +} + +/*--------------------------------------------------------------------*/ +template +diMuonMassBias::extremaMap getExtrema(const T& inputColl, const float sigma) +/*--------------------------------------------------------------------*/ +{ + diMuonMassBias::extremaMap result; // output + + // first transform the map vector of maps with a map of vectors + const auto& mapOfVecs = transformMaps(inputColl); + + for (const auto& [key, vec] : mapOfVecs) { + TObjArray* array = new TObjArray(vec.size()); + + for (const auto& histo : vec) { + array->Add(histo); + } + + Double_t theMaximum = (static_cast(array->At(0)))->GetMaximum(); + Double_t theMinimum = (static_cast(array->At(0)))->GetMinimum(); + for (Int_t i = 0; i < array->GetSize(); i++) { + if ((static_cast(array->At(i)))->GetMaximum() > theMaximum) { + theMaximum = (static_cast(array->At(i)))->GetMaximum(); + } + if ((static_cast(array->At(i)))->GetMinimum() < theMinimum) { + theMinimum = (static_cast(array->At(i)))->GetMinimum(); + } + } + delete array; + + const auto& delta = theMaximum - theMinimum; + result.insert({key, std::make_pair(theMinimum - (sigma * delta), theMaximum + (sigma * delta))}); + } + return result; +} + /*--------------------------------------------------------------------*/ std::pair getClosestFactors(int input) /*--------------------------------------------------------------------*/ @@ -99,32 +183,6 @@ void MakeNicePlotStyle(T* hist) hist->GetXaxis()->SetLabelSize(.05); } -namespace diMuonMassBias { - - struct fitOutputs { - public: - fitOutputs(const Measurement1D& bias, const Measurement1D& width) : m_bias(bias), m_width(width) {} - - // getters - const Measurement1D getBias() { return m_bias; } - const Measurement1D getWidth() { return m_width; } - bool isInvalid() const { - return (m_bias.value() == 0.f && m_bias.error() == 0.f && m_width.value() == 0.f && m_width.error() == 0.f); - } - - private: - Measurement1D m_bias; - Measurement1D m_width; - }; - - static constexpr int minimumHits = 10; - - using histoMap = std::map; - using histo2DMap = std::map; - using histo3DMap = std::map; - -} // namespace diMuonMassBias - //----------------------------------------------------------------------------------- diMuonMassBias::fitOutputs fitBWTimesCB(TH1* hist) //----------------------------------------------------------------------------------- @@ -520,6 +578,7 @@ void producePlots(const std::vector& inputMap, const std::vector& MEtoHarvest, const std::vector& labels, const TString& Rlabel, + const bool useAutoLimits, const bool isWidth) /************************************************/ { @@ -537,6 +596,9 @@ void producePlots(const std::vector& inputMap, infoBox->SetFillColor(kWhite); infoBox->SetTextSize(0.035); + // get the extrema + diMuonMassBias::extremaMap extrema = getExtrema(inputMap, 3.f); + for (const auto& var : MEtoHarvest) { TCanvas* c = new TCanvas( ((isWidth ? "width_" : "mean_") + var).c_str(), ((isWidth ? "width_" : "mean_") + var).c_str(), W, H); @@ -574,12 +636,18 @@ void producePlots(const std::vector& inputMap, histoMap.at(var)->SetMarkerColor(def_colors[count]); histoMap.at(var)->SetMarkerStyle(def_markers[count]); histoMap.at(var)->SetMarkerSize(1.5); - if (isWidth) { - // for width resolution between 0.5 and 2.8 - histoMap.at(var)->GetYaxis()->SetRangeUser(0.5, 2.85); + + // set the limits + if (!useAutoLimits) { + if (isWidth) { + // for width resolution between 0.5 and 2.8 + histoMap.at(var)->GetYaxis()->SetRangeUser(0.5, 2.85); + } else { + // for mass between 90.5 and 91.5 + histoMap.at(var)->GetYaxis()->SetRangeUser(90.5, 91.5); + } } else { - // for mass between 90.5 and 91.5 - histoMap.at(var)->GetYaxis()->SetRangeUser(90.5, 91.5); + histoMap.at(var)->GetYaxis()->SetRangeUser(extrema.at(var).first, extrema.at(var).second); } MakeNicePlotStyle(histoMap.at(var)); @@ -621,6 +689,7 @@ void produceMaps(const std::vector& inputMap, const std::vector& MEtoHarvest, const std::vector& labels, const TString& Rlabel, + const bool useAutoLimits, const bool isWidth) /************************************************/ { @@ -644,6 +713,9 @@ void produceMaps(const std::vector& inputMap, infoBox->SetFillColor(kWhite); infoBox->SetTextSize(0.035); + // get the extrema + diMuonMassBias::extremaMap extrema = getExtrema(inputMap, 0.f); + for (const auto& var : MEtoHarvest) { TCanvas* c = new TCanvas( ((isWidth ? "width_" : "mean_") + var).c_str(), ((isWidth ? "width_" : "mean_") + var).c_str(), W, H); @@ -676,12 +748,17 @@ void produceMaps(const std::vector& inputMap, } } - if (isWidth) { - // for width resolution between 0.5 and 2.8 - histoMap.at(var)->GetZaxis()->SetRangeUser(0.5, 2.85); + // set the limits + if (!useAutoLimits) { + if (isWidth) { + // for width resolution between 0.5 and 2.8 + histoMap.at(var)->GetZaxis()->SetRangeUser(0.5, 2.85); + } else { + // for mass between 90.5 and 91.5 + histoMap.at(var)->GetZaxis()->SetRangeUser(90., 92.); + } } else { - // for mass between 90.5 and 91.5 - histoMap.at(var)->GetZaxis()->SetRangeUser(90., 92.); + histoMap.at(var)->GetZaxis()->SetRangeUser(extrema.at(var).first, extrema.at(var).second); } MakeNicePlotStyle(histoMap.at(var)); @@ -716,7 +793,7 @@ void produceMaps(const std::vector& inputMap, } /************************************************/ -void DiMuonMassProfiles(TString namesandlabels, const TString& Rlabel = "") +void DiMuonMassProfiles(TString namesandlabels, const TString& Rlabel = "", const bool useAutoLimits = false) /************************************************/ { gStyle->SetOptStat(0); @@ -802,13 +879,13 @@ void DiMuonMassProfiles(TString namesandlabels, const TString& Rlabel = "") "TkTkMassVsPhiPlusInEtaBins/th2d_mass_PhiPlus_backward-backward", "TkTkMassVsPhiPlusInEtaBins/th2d_mass_PhiPlus_forward-backward"}; - producePlots(v_meanHistos, MEtoHarvest, labels, Rlabel, false); - producePlots(v_widthHistos, MEtoHarvest, labels, Rlabel, true); + producePlots(v_meanHistos, MEtoHarvest, labels, Rlabel, useAutoLimits, false); + producePlots(v_widthHistos, MEtoHarvest, labels, Rlabel, useAutoLimits, true); std::vector MEtoHarvest3D = {"th3d_mass_vs_eta_phi_plus", "th3d_mass_vs_eta_phi_minus"}; - produceMaps(v_meanMaps, MEtoHarvest3D, labels, Rlabel, false); - produceMaps(v_widthMaps, MEtoHarvest3D, labels, Rlabel, true); + produceMaps(v_meanMaps, MEtoHarvest3D, labels, Rlabel, useAutoLimits, false); + produceMaps(v_widthMaps, MEtoHarvest3D, labels, Rlabel, useAutoLimits, true); // finally close the file for (const auto& file : sourceFiles) {