OfflineValidation/macros/DiMuonMassProfiles.C

0001 // ROOT includes
0002 #include "RooAddPdf.h"
0003 #include "RooCBShape.h"
0004 #include "RooDataHist.h"
0005 #include "RooExponential.h"
0006 #include "RooGaussian.h"
0007 #include "RooPlot.h"
0008 #include "RooRealVar.h"
0009 #include "RooVoigtian.h"
0010 #include "TCanvas.h"
0011 #include "TClass.h"
0012 #include "TDirectory.h"
0013 #include "TFile.h"
0014 #include "TGaxis.h"
0015 #include "TH1.h"
0016 #include "TH2.h"
0017 #include "TH3.h"
0018 #include "TKey.h"
0019 #include "TLegend.h"
0020 #include "TObjString.h"
0021 #include "TObject.h"
0022 #include "TProfile.h"
0023 #include "TRatioPlot.h"
0024 #include "TStyle.h"
0025
0026 // standard includes
0027 #include <iomanip>
0028 #include <iostream>
0029 #include <map>
0030 #include <fmt/core.h>
0031
0032 // CMSSW classes / style
0033 #include "Alignment/OfflineValidation/interface/FitWithRooFit.h"
0034 #include "Alignment/OfflineValidation/macros/CMS_lumi.h"
0035 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
0036
0037 bool debugMode_{false};
0038 Int_t def_markers[9] = {kFullSquare,
0039                         kFullCircle,
0040                         kFullTriangleDown,
0041                         kOpenSquare,
0042                         kDot,
0043                         kOpenCircle,
0044                         kFullTriangleDown,
0045                         kFullTriangleUp,
0046                         kOpenTriangleDown};
0047
0048 Int_t def_colors[9] = {kBlack, kRed, kBlue, kMagenta, kGreen, kCyan, kViolet, kOrange, kGreen + 2};
0049
0050 namespace diMuonMassBias {
0051
0052   struct fitOutputs {
0053   public:
0054     fitOutputs(const Measurement1D& bias, const Measurement1D& width) : m_bias(bias), m_width(width) {}
0055
0056     // getters
0057     const Measurement1D getBias() { return m_bias; }
0058     const Measurement1D getWidth() { return m_width; }
0059     bool isInvalid() const {
0060       return (m_bias.value() == 0.f && m_bias.error() == 0.f && m_width.value() == 0.f && m_width.error() == 0.f);
0061     }
0062
0063   private:
0064     Measurement1D m_bias;
0065     Measurement1D m_width;
0066   };
0067
0068   static constexpr int minimumHits = 10;
0069
0070   using histoMap = std::map<std::string, TH1F*>;
0071   using histo2DMap = std::map<std::string, TH2F*>;
0072   using histo3DMap = std::map<std::string, TH3F*>;
0073
0074   using extrema = std::pair<Double_t, Double_t>;
0075   using extremaMap = std::map<std::string, diMuonMassBias::extrema>;
0076
0077 }  // namespace diMuonMassBias
0078
0079 /*--------------------------------------------------------------------*/
0080 template <typename T>
0081 std::map<std::string, std::vector<T>> transformMaps(const std::vector<std::map<std::string, T>>& vecOfMaps)
0082 /*--------------------------------------------------------------------*/
0083 {
0084   std::map<std::string, std::vector<T>> result;
0085
0086   // Lambda to insert each key-value pair into the result map
0087   auto insert_into_result = [&result](const std::map<std::string, T>& map) {
0088     for (const auto& pair : map) {
0089       result[pair.first].push_back(pair.second);
0090     }
0091   };
0092
0093   // Apply the lambda to each map in the vector
0094   std::for_each(vecOfMaps.begin(), vecOfMaps.end(), insert_into_result);
0095
0096   return result;
0097 }
0098
0099 /*--------------------------------------------------------------------*/
0100 template <typename T>
0101 diMuonMassBias::extremaMap getExtrema(const T& inputColl, const float sigma)
0102 /*--------------------------------------------------------------------*/
0103 {
0104   diMuonMassBias::extremaMap result;  // output
0105
0106   // first transform the map vector of maps with a map of vectors
0107   const auto& mapOfVecs = transformMaps(inputColl);
0108
0109   for (const auto& [key, vec] : mapOfVecs) {
0110     TObjArray* array = new TObjArray(vec.size());
0111
0112     for (const auto& histo : vec) {
0113       array->Add(histo);
0114     }
0115
0116     Double_t theMaximum = (static_cast<TH1*>(array->At(0)))->GetMaximum();
0117     Double_t theMinimum = (static_cast<TH1*>(array->At(0)))->GetMinimum();
0118     for (Int_t i = 0; i < array->GetSize(); i++) {
0119       if ((static_cast<TH1*>(array->At(i)))->GetMaximum() > theMaximum) {
0120         theMaximum = (static_cast<TH1*>(array->At(i)))->GetMaximum();
0121       }
0122       if ((static_cast<TH1*>(array->At(i)))->GetMinimum() < theMinimum) {
0123         theMinimum = (static_cast<TH1*>(array->At(i)))->GetMinimum();
0124       }
0125     }
0126     delete array;
0127
0128     const auto& delta = theMaximum - theMinimum;
0129     result.insert({key, std::make_pair(theMinimum - (sigma * delta), theMaximum + (sigma * delta))});
0130   }
0131   return result;
0132 }
0133
0134 /*--------------------------------------------------------------------*/
0135 std::pair<int, int> getClosestFactors(int input)
0136 /*--------------------------------------------------------------------*/
0137 {
0138   if ((input % 2 != 0) && input > 1) {
0139     input += 1;
0140   }
0141
0142   int testNum = (int)sqrt(input);
0143   while (input % testNum != 0) {
0144     testNum--;
0145   }
0146   return std::make_pair(testNum, input / testNum);
0147 }
0148
0149 /*--------------------------------------------------------------------*/
0150 template <typename T>
0151 void MakeNicePlotStyle(T* hist)
0152 /*--------------------------------------------------------------------*/
0153 {
0154   //hist->SetStats(kFALSE);
0155   hist->SetLineWidth(2);
0156   hist->GetXaxis()->SetNdivisions(505);
0157   hist->GetXaxis()->CenterTitle(true);
0158   hist->GetYaxis()->CenterTitle(true);
0159   hist->GetXaxis()->SetTitleFont(42);
0160   hist->GetYaxis()->SetTitleFont(42);
0161   hist->GetXaxis()->SetTitleSize(0.05);
0162   hist->GetYaxis()->SetTitleSize(0.05);
0163   if (((TObject*)hist)->IsA()->InheritsFrom("TH2")) {
0164     hist->GetZaxis()->CenterTitle(true);
0165     hist->GetZaxis()->SetTitleSize(0.04);
0166     hist->GetZaxis()->SetLabelFont(42);
0167     hist->GetYaxis()->SetLabelSize(.05);
0168     hist->GetXaxis()->SetTitleOffset(1.0);
0169     hist->GetYaxis()->SetTitleOffset(1.0);
0170     hist->GetZaxis()->SetTitleOffset(1.15);
0171   } else {
0172     hist->GetXaxis()->SetTitleOffset(0.9);
0173     hist->GetYaxis()->SetTitleOffset(1.7);
0174   }
0175   hist->GetXaxis()->SetLabelFont(42);
0176   hist->GetYaxis()->SetLabelFont(42);
0177   if (((TObject*)hist)->IsA()->InheritsFrom("TGraph")) {
0178     hist->GetYaxis()->SetLabelSize(.025);
0179     //hist->GetYaxis()->SetNdivisions(505);
0180   } else {
0181     hist->GetYaxis()->SetLabelSize(.05);
0182   }
0183   hist->GetXaxis()->SetLabelSize(.05);
0184 }
0185
0186 //-----------------------------------------------------------------------------------
0187 diMuonMassBias::fitOutputs fitBWTimesCB(TH1* hist)
0188 //-----------------------------------------------------------------------------------
0189 {
0190   if (hist->GetEntries() < diMuonMassBias::minimumHits) {
0191     std::cout << " Input histogram:" << hist->GetName() << " has not enough entries (" << hist->GetEntries()
0192               << ") for a meaningful Voigtian fit!\n"
0193               << "Skipping!";
0194
0195     return diMuonMassBias::fitOutputs(Measurement1D(0., 0.), Measurement1D(0., 0.));
0196   }
0197
0198   TCanvas* c1 = new TCanvas();
0199   if (debugMode_) {
0200     c1->Clear();
0201     c1->SetLeftMargin(0.15);
0202     c1->SetRightMargin(0.10);
0203   }
0204
0205   // silence messages
0206   RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
0207
0208   Double_t xMean = 91.1876;
0209   Double_t xMin = hist->GetXaxis()->GetXmin();
0210   Double_t xMax = hist->GetXaxis()->GetXmax();
0211
0212   if (debugMode_) {
0213     std::cout << " fitting range: (" << xMin << "-" << xMax << ")" << std::endl;
0214   }
0215
0216   double sigma(2.);
0217   double sigmaMin(0.1);
0218   double sigmaMax(10.);
0219
0220   double sigma2(0.1);
0221   double sigma2Min(0.);
0222   double sigma2Max(10.);
0223
0224   std::unique_ptr<FitWithRooFit> fitter = std::make_unique<FitWithRooFit>();
0225
0226   bool useChi2(false);
0227
0228   fitter->useChi2_ = useChi2;
0229   fitter->initMean(xMean, xMin, xMax);
0230   fitter->initSigma(sigma, sigmaMin, sigmaMax);
0231   fitter->initSigma2(sigma2, sigma2Min, sigma2Max);
0232   fitter->initAlpha(1.5, 0.05, 10.);
0233   fitter->initN(1, 0.01, 100.);
0234   fitter->initFGCB(0.4, 0., 1.);
0235   fitter->initGamma(2.4952, 0., 10.);
0236   fitter->gamma()->setConstant(kTRUE);
0237   fitter->initMean2(0., -20., 20.);
0238   fitter->mean2()->setConstant(kTRUE);
0239   fitter->initSigma(1.2, 0., 5.);
0240   fitter->initAlpha(1.5, 0.05, 10.);
0241   fitter->initN(1, 0.01, 100.);
0242   fitter->initExpCoeffA0(-1., -10., 10.);
0243   fitter->initExpCoeffA1(0., -10., 10.);
0244   fitter->initExpCoeffA2(0., -2., 2.);
0245   fitter->initFsig(0.9, 0., 1.);
0246   fitter->initA0(0., -10., 10.);
0247   fitter->initA1(0., -10., 10.);
0248   fitter->initA2(0., -10., 10.);
0249   fitter->initA3(0., -10., 10.);
0250   fitter->initA4(0., -10., 10.);
0251   fitter->initA5(0., -10., 10.);
0252   fitter->initA6(0., -10., 10.);
0253
0254   fitter->fit(hist, "breitWignerTimesCB", "exponential", xMin, xMax, false);
0255
0256   TString histName = hist->GetName();
0257
0258   if (debugMode_) {
0259     c1->Print("fit_debug" + histName + ".pdf");
0260   }
0261   delete c1;
0262
0263   Measurement1D resultM(fitter->mean()->getValV(), fitter->mean()->getError());
0264   Measurement1D resultW(fitter->sigma()->getValV(), fitter->sigma()->getError());
0265
0266   return diMuonMassBias::fitOutputs(resultM, resultW);
0267 }
0268
0269 //-----------------------------------------------------------------------------------
0270 void fitAndFillHisto(std::pair<std::string, TH2F*> toHarvest,
0271                      diMuonMassBias::histoMap& meanHistos_,
0272                      diMuonMassBias::histoMap& widthHistos_)
0273 //-----------------------------------------------------------------------------------
0274 {
0275   const auto& key = toHarvest.first;
0276   const auto& ME = toHarvest.second;
0277
0278   if (debugMode_)
0279     std::cout << "dealing with key: " << key << std::endl;
0280
0281   if (ME == nullptr) {
0282     std::cout << "could not find MonitorElement for key: " << key << std::endl;
0283     return;
0284   }
0285
0286   for (int bin = 1; bin <= ME->GetNbinsY(); bin++) {
0287     const auto& yaxis = ME->GetYaxis();
0288     const auto& low_edge = yaxis->GetBinLowEdge(bin);
0289     const auto& high_edge = yaxis->GetBinUpEdge(bin);
0290
0291     if (debugMode_)
0292       std::cout << "dealing with bin: " << bin << " range: (" << low_edge << "," << high_edge << ")";
0293     TH1D* Proj = ME->ProjectionX(Form("%s_proj_%i", key.c_str(), bin), bin, bin);
0294     Proj->SetTitle(Form("%s #in (%.2f,%.2f), bin: %i", Proj->GetTitle(), low_edge, high_edge, bin));
0295
0296     diMuonMassBias::fitOutputs results = fitBWTimesCB(Proj);
0297
0298     if (results.isInvalid()) {
0299       std::cout << "the current bin has invalid data" << std::endl;
0300       continue;
0301     }
0302
0303     // fill the mean profiles
0304     const Measurement1D& bias = results.getBias();
0305     meanHistos_[key]->SetBinContent(bin, bias.value());
0306     meanHistos_[key]->SetBinError(bin, bias.error());
0307
0308     // fill the width profiles
0309     const Measurement1D& width = results.getWidth();
0310     widthHistos_[key]->SetBinContent(bin, width.value());
0311     widthHistos_[key]->SetBinError(bin, width.error());
0312
0313     if (debugMode_) {
0314       std::cout << "key: " << key << " bin: " << bin << " bias: " << bias.value()
0315                 << " (in histo: " << meanHistos_[key]->GetBinContent(bin) << ") width: " << width.value()
0316                 << " (in histo: " << widthHistos_[key]->GetBinContent(bin) << ")" << std::endl;
0317     }
0318   }
0319 }
0320
0321 //-----------------------------------------------------------------------------------
0322 void bookMaps(const diMuonMassBias::histo3DMap& harvestTargets_,
0323               diMuonMassBias::histo2DMap& meanMaps_,
0324               diMuonMassBias::histo2DMap& widthMaps_,
0325               const unsigned int counter)
0326 //-----------------------------------------------------------------------------------
0327 {
0328   for (const auto& [key, ME] : harvestTargets_) {
0329     if (ME == nullptr) {
0330       std::cout << "could not find MonitorElement for key: " << key << std::endl;
0331       continue;
0332     }
0333
0334     const auto& title = ME->GetTitle();
0335     const auto& xtitle = ME->GetYaxis()->GetTitle();
0336     const auto& ytitle = ME->GetZaxis()->GetTitle();
0337
0338     const auto& nxbins = ME->GetNbinsY();
0339     const auto& xmin = ME->GetYaxis()->GetXmin();
0340     const auto& xmax = ME->GetYaxis()->GetXmax();
0341
0342     const auto& nybins = ME->GetNbinsZ();
0343     const auto& ymin = ME->GetZaxis()->GetXmin();
0344     const auto& ymax = ME->GetZaxis()->GetXmax();
0345
0346     if (debugMode_) {
0347       std::cout << "Booking " << key << std::endl;
0348     }
0349
0350     TH2F* meanToBook =
0351         new TH2F(fmt::format("Mean_{}_{}", counter, key).c_str(),
0352                  fmt::format("{};{};{};#LT M_{{#mu^{{-}}#mu^{{+}}}} #GT [GeV]", title, xtitle, ytitle).c_str(),
0353                  nxbins,
0354                  xmin,
0355                  xmax,
0356                  nybins,
0357                  ymin,
0358                  ymax);
0359
0360     if (debugMode_) {
0361       std::cout << "after creating mean" << key << std::endl;
0362     }
0363
0364     meanMaps_.insert({key, meanToBook});
0365
0366     if (debugMode_) {
0367       std::cout << "after inserting mean" << key << std::endl;
0368     }
0369
0370     TH2F* sigmaToBook =
0371         new TH2F(fmt::format("Sigma_{}_{}", counter, key).c_str(),
0372                  fmt::format("{};{};#sigma of  M_{{#mu^{{-}}#mu^{{+}}}} [GeV]", title, xtitle, ytitle).c_str(),
0373                  nxbins,
0374                  xmin,
0375                  xmax,
0376                  nybins,
0377                  ymin,
0378                  ymax);
0379
0380     if (debugMode_) {
0381       std::cout << "after creating sigma" << key << std::endl;
0382     }
0383
0384     widthMaps_.insert({key, sigmaToBook});
0385
0386     if (debugMode_) {
0387       std::cout << "after inserting sigma" << key << std::endl;
0388     }
0389   }
0390 }
0391
0392 //-----------------------------------------------------------------------------------
0393 void fitAndFillMap(std::pair<std::string, TH3F*> toHarvest,
0394                    diMuonMassBias::histo2DMap& meanMaps_,
0395                    diMuonMassBias::histo2DMap& widthMaps_)
0396 //-----------------------------------------------------------------------------------
0397 {
0398   const auto& key = toHarvest.first;
0399   const auto& ME = toHarvest.second;
0400
0401   if (debugMode_)
0402     std::cout << "dealing with key: " << key << std::endl;
0403
0404   if (ME == nullptr) {
0405     std::cout << "could not find MonitorElement for key: " << key << std::endl;
0406     return;
0407   }
0408
0409   for (int binY = 1; binY <= ME->GetNbinsY(); binY++) {
0410     const auto& yaxis = ME->GetYaxis();
0411     const auto& y_low_edge = yaxis->GetBinLowEdge(binY);
0412     const auto& y_high_edge = yaxis->GetBinUpEdge(binY);
0413
0414     for (int binZ = 1; binZ <= ME->GetNbinsZ(); binZ++) {
0415       const auto& zaxis = ME->GetZaxis();
0416       const auto& z_low_edge = zaxis->GetBinLowEdge(binZ);
0417       const auto& z_high_edge = zaxis->GetBinUpEdge(binZ);
0418
0419       if (debugMode_) {
0420         std::cout << "dealing with y bin: " << binY << " range: (" << y_low_edge << "," << y_high_edge << ")";
0421         std::cout << " dealing with z bin: " << binZ << " range: (" << z_low_edge << "," << z_high_edge << ")"
0422                   << std::endl;
0423       }
0424
0425       const auto& ProjYZ = ME->ProjectionX(Form("%s_proj_%i_proj_%i", key.c_str(), binY, binZ), binY, binY, binZ, binZ);
0426       ProjYZ->SetTitle(Form("%s #in (%.2f,%.2f) - (%.2f,%.2f), bin: %i-%i ",
0427                             ProjYZ->GetTitle(),
0428                             y_low_edge,
0429                             y_high_edge,
0430                             z_low_edge,
0431                             z_high_edge,
0432                             binY,
0433                             binZ));
0434
0435       diMuonMassBias::fitOutputs results = fitBWTimesCB(ProjYZ);
0436
0437       if (results.isInvalid()) {
0438         std::cout << "the current bin has invalid data" << std::endl;
0439         continue;
0440       }
0441
0442       // fill the mean profiles
0443       const Measurement1D& bias = results.getBias();
0444       meanMaps_[key]->SetBinContent(binY, binZ, bias.value());
0445       meanMaps_[key]->SetBinError(binY, binZ, bias.error());
0446
0447       // fill the width profiles
0448       const Measurement1D& width = results.getWidth();
0449       widthMaps_[key]->SetBinContent(binY, binZ, width.value());
0450       widthMaps_[key]->SetBinError(binY, binZ, width.error());
0451
0452       if (debugMode_) {
0453         std::cout << "key: " << key << " bin: (" << binY << "," << binZ << ")  bias: " << bias.value()
0454                   << " (in histo: " << meanMaps_[key]->GetBinContent(binY, binZ) << ") width: " << width.value()
0455                   << " (in histo: " << widthMaps_[key]->GetBinContent(binY, binZ) << ")" << std::endl;
0456       }
0457     }
0458   }
0459 }
0460
0461 //-----------------------------------------------------------------------------------
0462 void bookHistos(const diMuonMassBias::histo2DMap& harvestTargets_,
0463                 diMuonMassBias::histoMap& meanHistos_,
0464                 diMuonMassBias::histoMap& widthHistos_,
0465                 const unsigned int counter)
0466 //-----------------------------------------------------------------------------------
0467 {
0468   for (const auto& [key, ME] : harvestTargets_) {
0469     if (ME == nullptr) {
0470       std::cout << "could not find MonitorElement for key: " << key << std::endl;
0471       continue;
0472     }
0473
0474     const auto& title = ME->GetTitle();
0475     const auto& xtitle = ME->GetYaxis()->GetTitle();
0476     const auto& ytitle = ME->GetXaxis()->GetTitle();
0477
0478     const auto& nxbins = ME->GetNbinsY();
0479     const auto& xmin = ME->GetYaxis()->GetXmin();
0480     const auto& xmax = ME->GetYaxis()->GetXmax();
0481
0482     if (debugMode_) {
0483       std::cout << "Booking " << key << std::endl;
0484     }
0485
0486     TH1F* meanToBook = new TH1F(fmt::format("Mean_{}_{}", counter, key).c_str(),
0487                                 fmt::format("{};{};#LT M_{{#mu^{{-}}#mu^{{+}}}} #GT [GeV]", title, xtitle).c_str(),
0488                                 nxbins,
0489                                 xmin,
0490                                 xmax);
0491
0492     if (debugMode_) {
0493       std::cout << "after creating mean" << key << std::endl;
0494     }
0495
0496     meanHistos_.insert({key, meanToBook});
0497
0498     if (debugMode_) {
0499       std::cout << "after inserting mean" << key << std::endl;
0500     }
0501
0502     TH1F* sigmaToBook = new TH1F(fmt::format("Sigma_{}_{}", counter, key).c_str(),
0503                                  fmt::format("{};{};#sigma of  M_{{#mu^{{-}}#mu^{{+}}}} [GeV]", title, xtitle).c_str(),
0504                                  nxbins,
0505                                  xmin,
0506                                  xmax);
0507
0508     if (debugMode_) {
0509       std::cout << "after creating sigma" << key << std::endl;
0510     }
0511
0512     widthHistos_.insert({key, sigmaToBook});
0513
0514     if (debugMode_) {
0515       std::cout << "after inserting sigma" << key << std::endl;
0516     }
0517   }
0518 }
0519
0520 //-----------------------------------------------------------------------------------
0521 void getMEsToHarvest(diMuonMassBias::histo2DMap& harvestTargets_,
0522                      diMuonMassBias::histo3DMap& harvestTargets3D_,
0523                      TFile* file)
0524 //-----------------------------------------------------------------------------------
0525 {
0526   std::string inFolder = "DiMuonMassValidation";
0527
0528   // list of 2D histograms to flatten
0529   std::vector<std::string> MEtoHarvest_ = {"th2d_mass_CosThetaCS",
0530                                            "th2d_mass_DeltaEta",
0531                                            "th2d_mass_EtaMinus",
0532                                            "th2d_mass_EtaPlus",
0533                                            "th2d_mass_PhiCS",
0534                                            "th2d_mass_PhiMinus",
0535                                            "th2d_mass_PhiPlus",
0536                                            "TkTkMassVsPhiMinusInEtaBins/th2d_mass_PhiMinus_barrel-barrel",
0537                                            "TkTkMassVsPhiMinusInEtaBins/th2d_mass_PhiMinus_barrel-forward",
0538                                            "TkTkMassVsPhiMinusInEtaBins/th2d_mass_PhiMinus_barrel-backward",
0539                                            "TkTkMassVsPhiMinusInEtaBins/th2d_mass_PhiMinus_forward-forward",
0540                                            "TkTkMassVsPhiMinusInEtaBins/th2d_mass_PhiMinus_backward-backward",
0541                                            "TkTkMassVsPhiMinusInEtaBins/th2d_mass_PhiMinus_forward-backward",
0542                                            "TkTkMassVsPhiPlusInEtaBins/th2d_mass_PhiPlus_barrel-barrel",
0543                                            "TkTkMassVsPhiPlusInEtaBins/th2d_mass_PhiPlus_barrel-forward",
0544                                            "TkTkMassVsPhiPlusInEtaBins/th2d_mass_PhiPlus_barrel-backward",
0545                                            "TkTkMassVsPhiPlusInEtaBins/th2d_mass_PhiPlus_forward-forward",
0546                                            "TkTkMassVsPhiPlusInEtaBins/th2d_mass_PhiPlus_backward-backward",
0547                                            "TkTkMassVsPhiPlusInEtaBins/th2d_mass_PhiPlus_forward-backward"};
0548
0549   //loop on the list of histograms to harvest
0550   for (const auto& hname : MEtoHarvest_) {
0551     std::cout << "trying to get: " << hname << std::endl;
0552     TH2F* toHarvest = (TH2F*)file->Get((inFolder + "/" + hname).c_str());
0553
0554     if (toHarvest == nullptr) {
0555       std::cout << "could not find input MonitorElement: " << inFolder + "/" + hname << std::endl;
0556       continue;
0557     }
0558     harvestTargets_.insert({hname, toHarvest});
0559   }
0560
0561   // list of 3D histograms to flatten
0562   std::vector<std::string> MEtoHarvest3D_ = {"th3d_mass_vs_eta_phi_plus", "th3d_mass_vs_eta_phi_minus"};
0563
0564   for (const auto& hname : MEtoHarvest3D_) {
0565     std::cout << "trying to get: " << hname << std::endl;
0566     TH3F* toHarvest = (TH3F*)file->Get((inFolder + "/" + hname).c_str());
0567
0568     if (toHarvest == nullptr) {
0569       std::cout << "could not find input MonitorElement: " << inFolder + "/" + hname << std::endl;
0570       continue;
0571     }
0572     harvestTargets3D_.insert({hname, toHarvest});
0573   }
0574 }
0575
0576 /************************************************/
0577 void producePlots(const std::vector<diMuonMassBias::histoMap>& inputMap,
0578                   const std::vector<std::string>& MEtoHarvest,
0579                   const std::vector<TString>& labels,
0580                   const TString& Rlabel,
0581                   const bool useAutoLimits,
0582                   const bool isWidth)
0583 /************************************************/
0584 {
0585   int W = 800;
0586   int H = 800;
0587   // references for T, B, L, R
0588   float T = 0.08 * H;
0589   float B = 0.12 * H;
0590   float L = 0.12 * W;
0591   float R = 0.04 * W;
0592
0593   // Draw the legend
0594   TLegend* infoBox = new TLegend(0.65, 0.75, 0.95, 0.90, "");
0595   infoBox->SetShadowColor(0);  // 0 = transparent
0596   infoBox->SetFillColor(kWhite);
0597   infoBox->SetTextSize(0.035);
0598
0599   // get the extrema
0600   diMuonMassBias::extremaMap extrema = getExtrema(inputMap, 3.f);
0601
0602   for (const auto& var : MEtoHarvest) {
0603     TCanvas* c = new TCanvas(
0604         ((isWidth ? "width_" : "mean_") + var).c_str(), ((isWidth ? "width_" : "mean_") + var).c_str(), W, H);
0605     c->SetFillColor(0);
0606     c->SetBorderMode(0);
0607     c->SetFrameFillStyle(0);
0608     c->SetFrameBorderMode(0);
0609     c->SetLeftMargin(L / W + 0.05);
0610     c->SetRightMargin(R / W);
0611     c->SetTopMargin(T / H);
0612     c->SetBottomMargin(B / H);
0613     c->SetTickx(0);
0614     c->SetTicky(0);
0615     c->SetGrid();
0616
0617     unsigned int count{0};
0618
0619     for (const auto& histoMap : inputMap) {
0620       if (debugMode_) {
0621         std::cout << var << "  n.bins: " << histoMap.at(var)->GetNbinsX()
0622                   << " entries: " << histoMap.at(var)->GetEntries() << "title: " << histoMap.at(var)->GetTitle()
0623                   << " x-axis title: " << histoMap.at(var)->GetXaxis()->GetTitle() << std::endl;
0624       }
0625
0626       if (debugMode_) {
0627         for (int bin = 1; bin <= histoMap.at(var)->GetNbinsX(); bin++) {
0628           std::cout << var << " bin " << bin << " : " << histoMap.at(var)->GetBinContent(bin) << " +/-"
0629                     << histoMap.at(var)->GetBinError(bin) << std::endl;
0630         }
0631       }
0632
0633       //histoMap.at(var)->SaveAs((var+".root").c_str());
0634
0635       histoMap.at(var)->SetLineColor(def_colors[count]);
0636       histoMap.at(var)->SetMarkerColor(def_colors[count]);
0637       histoMap.at(var)->SetMarkerStyle(def_markers[count]);
0638       histoMap.at(var)->SetMarkerSize(1.5);
0639
0640       // set the limits
0641       if (!useAutoLimits) {
0642         if (isWidth) {
0643           // for width resolution between 0.5 and 2.8
0644           histoMap.at(var)->GetYaxis()->SetRangeUser(0.5, 2.85);
0645         } else {
0646           // for mass between 90.5 and 91.5
0647           histoMap.at(var)->GetYaxis()->SetRangeUser(90.5, 91.5);
0648         }
0649       } else {
0650         histoMap.at(var)->GetYaxis()->SetRangeUser(extrema.at(var).first, extrema.at(var).second);
0651       }
0652
0653       MakeNicePlotStyle<TH1>(histoMap.at(var));
0654
0655       c->cd();
0656       if (count == 0) {
0657         histoMap.at(var)->Draw("E1");
0658       } else {
0659         histoMap.at(var)->Draw("E1same");
0660       }
0661
0662       // fill the legend only if that's the first element in the vector of variables
0663       if (var == MEtoHarvest[0]) {
0664         infoBox->AddEntry(histoMap.at(var), labels[count], "LP");
0665       }
0666       infoBox->Draw("same");
0667       count++;
0668     }
0669
0670     CMS_lumi(c, 0, 3, Rlabel);
0671
0672     // Find the position of the first '/'
0673     size_t pos = var.find('/');
0674     std::string outputName{var};
0675
0676     // Check if '/' is found
0677     if (pos != std::string::npos) {
0678       // Erase the substring before the '/' (including the '/')
0679       outputName.erase(0, pos + 1);
0680     }
0681
0682     c->SaveAs(((isWidth ? "width_" : "mean_") + outputName + ".png").c_str());
0683     c->SaveAs(((isWidth ? "width_" : "mean_") + outputName + ".pdf").c_str());
0684   }
0685 }
0686
0687 /************************************************/
0688 void produceMaps(const std::vector<diMuonMassBias::histo2DMap>& inputMap,
0689                  const std::vector<std::string>& MEtoHarvest,
0690                  const std::vector<TString>& labels,
0691                  const TString& Rlabel,
0692                  const bool useAutoLimits,
0693                  const bool isWidth)
0694 /************************************************/
0695 {
0696   const auto& sides = getClosestFactors(labels.size());
0697
0698   if (debugMode_) {
0699     std::cout << "SIDES:" << sides.second << " :" << sides.first << std::endl;
0700   }
0701
0702   int W = 800 * sides.second;
0703   int H = 600 * sides.first;
0704   // references for T, B, L, R
0705   float T = 0.08 * H;
0706   float B = 0.12 * H;
0707   float L = 0.12 * W;
0708   float R = 0.04 * W;
0709
0710   // Draw the legend
0711   TLegend* infoBox = new TLegend(0.65, 0.75, 0.95, 0.90, "");
0712   infoBox->SetShadowColor(0);  // 0 = transparent
0713   infoBox->SetFillColor(kWhite);
0714   infoBox->SetTextSize(0.035);
0715
0716   // get the extrema
0717   diMuonMassBias::extremaMap extrema = getExtrema(inputMap, 0.f);
0718
0719   for (const auto& var : MEtoHarvest) {
0720     TCanvas* c = new TCanvas(
0721         ((isWidth ? "width_" : "mean_") + var).c_str(), ((isWidth ? "width_" : "mean_") + var).c_str(), W, H);
0722
0723     c->SetFillColor(0);
0724     c->SetBorderMode(0);
0725     c->SetFrameFillStyle(0);
0726     c->SetFrameBorderMode(0);
0727     c->SetLeftMargin(L / W);
0728     c->SetRightMargin(R / W);
0729     c->SetTopMargin(T / H);
0730     c->SetBottomMargin(B / H);
0731     c->SetTickx(0);
0732     c->SetTicky(0);
0733
0734     c->Divide(sides.second, sides.first);
0735
0736     unsigned int count{0};
0737     for (const auto& histoMap : inputMap) {
0738       if (debugMode_) {
0739         std::cout << var << "  n.bins: " << histoMap.at(var)->GetNbinsX()
0740                   << " entries: " << histoMap.at(var)->GetEntries() << "title: " << histoMap.at(var)->GetTitle()
0741                   << " x-axis title: " << histoMap.at(var)->GetXaxis()->GetTitle() << std::endl;
0742       }
0743
0744       if (debugMode_) {
0745         for (int bin = 1; bin <= histoMap.at(var)->GetNbinsX(); bin++) {
0746           std::cout << var << " bin " << bin << " : " << histoMap.at(var)->GetBinContent(bin) << " +/-"
0747                     << histoMap.at(var)->GetBinError(bin) << std::endl;
0748         }
0749       }
0750
0751       // set the limits
0752       if (!useAutoLimits) {
0753         if (isWidth) {
0754           // for width resolution between 0.5 and 2.8
0755           histoMap.at(var)->GetZaxis()->SetRangeUser(0.5, 2.85);
0756         } else {
0757           // for mass between 90.5 and 91.5
0758           histoMap.at(var)->GetZaxis()->SetRangeUser(90., 92.);
0759         }
0760       } else {
0761         histoMap.at(var)->GetZaxis()->SetRangeUser(extrema.at(var).first, extrema.at(var).second);
0762       }
0763
0764       MakeNicePlotStyle<TH2>(histoMap.at(var));
0765
0766       c->cd(count + 1);
0767       c->cd(count + 1)->SetRightMargin(0.15);
0768       histoMap.at(var)->Draw("colz");
0769
0770       // fill the legend only if that's the first element in the vector of variables
0771       if (var == MEtoHarvest[0]) {
0772         infoBox->AddEntry(histoMap.at(var), labels[count], "LP");
0773       }
0774       //infoBox->Draw("same");
0775
0776       CMS_lumi(dynamic_cast<TPad*>(gPad), 0, 3, labels[count]);
0777       count++;
0778     }
0779
0780     // Find the position of the first '/'
0781     size_t pos = var.find('/');
0782     std::string outputName{var};
0783
0784     // Check if '/' is found
0785     if (pos != std::string::npos) {
0786       // Erase the substring before the '/' (including the '/')
0787       outputName.erase(0, pos + 1);
0788     }
0789
0790     c->SaveAs(((isWidth ? "width_" : "mean_") + outputName + ".png").c_str());
0791     c->SaveAs(((isWidth ? "width_" : "mean_") + outputName + ".pdf").c_str());
0792   }
0793 }
0794
0795 /************************************************/
0796 void DiMuonMassProfiles(TString namesandlabels, const TString& Rlabel = "", const bool useAutoLimits = false)
0797 /************************************************/
0798 {
0799   gStyle->SetOptStat(0);
0800   gStyle->SetOptTitle(0);
0801
0802   std::vector<TString> labels;
0803   std::vector<TFile*> sourceFiles;
0804
0805   namesandlabels.Remove(TString::kTrailing, ',');
0806   TObjArray* nameandlabelpairs = namesandlabels.Tokenize(",");
0807   for (Int_t i = 0; i < nameandlabelpairs->GetEntries(); ++i) {
0808     TObjArray* aFileLegPair = TString(nameandlabelpairs->At(i)->GetName()).Tokenize("=");
0809     if (aFileLegPair->GetEntries() == 2) {
0810       sourceFiles.push_back(TFile::Open(aFileLegPair->At(0)->GetName(), "READ"));
0811       TObjString* s_label = (TObjString*)aFileLegPair->At(1);
0812       labels.push_back(s_label->String());
0813     } else {
0814       std::cout << "Please give file name and legend entry in the following form:\n"
0815                 << " filename1=legendentry1,filename2=legendentry2\n";
0816       return;
0817     }
0818   }
0819
0820   // for the bias plots
0821   std::vector<diMuonMassBias::histoMap> v_meanHistos;
0822   std::vector<diMuonMassBias::histoMap> v_widthHistos;
0823
0824   // for the maps
0825   std::vector<diMuonMassBias::histo2DMap> v_meanMaps;
0826   std::vector<diMuonMassBias::histo2DMap> v_widthMaps;
0827
0828   unsigned int countFiles{0};
0829   for (const auto& file : sourceFiles) {
0830     diMuonMassBias::histo2DMap harvestTargets;
0831     diMuonMassBias::histo3DMap harvestTargets3D;
0832
0833     getMEsToHarvest(harvestTargets, harvestTargets3D, file);
0834
0835     diMuonMassBias::histoMap meanHistos;
0836     diMuonMassBias::histoMap widthHistos;
0837
0838     bookHistos(harvestTargets, meanHistos, widthHistos, countFiles);
0839
0840     diMuonMassBias::histo2DMap meanMaps;
0841     diMuonMassBias::histo2DMap widthMaps;
0842
0843     bookMaps(harvestTargets3D, meanMaps, widthMaps, countFiles);
0844
0845     for (const auto& element : harvestTargets) {
0846       fitAndFillHisto(element, meanHistos, widthHistos);
0847     }
0848
0849     v_meanHistos.push_back(meanHistos);
0850     v_widthHistos.push_back(widthHistos);
0851
0852     for (const auto& element : harvestTargets3D) {
0853       fitAndFillMap(element, meanMaps, widthMaps);
0854     }
0855
0856     v_meanMaps.push_back(meanMaps);
0857     v_widthMaps.push_back(widthMaps);
0858
0859     countFiles++;
0860   }
0861
0862   // now do the plotting
0863   std::vector<std::string> MEtoHarvest = {"th2d_mass_CosThetaCS",
0864                                           "th2d_mass_DeltaEta",
0865                                           "th2d_mass_EtaMinus",
0866                                           "th2d_mass_EtaPlus",
0867                                           "th2d_mass_PhiCS",
0868                                           "th2d_mass_PhiMinus",
0869                                           "th2d_mass_PhiPlus",
0870                                           "TkTkMassVsPhiMinusInEtaBins/th2d_mass_PhiMinus_barrel-barrel",
0871                                           "TkTkMassVsPhiMinusInEtaBins/th2d_mass_PhiMinus_barrel-forward",
0872                                           "TkTkMassVsPhiMinusInEtaBins/th2d_mass_PhiMinus_barrel-backward",
0873                                           "TkTkMassVsPhiMinusInEtaBins/th2d_mass_PhiMinus_forward-forward",
0874                                           "TkTkMassVsPhiMinusInEtaBins/th2d_mass_PhiMinus_backward-backward",
0875                                           "TkTkMassVsPhiPlusInEtaBins/th2d_mass_PhiPlus_barrel-barrel",
0876                                           "TkTkMassVsPhiPlusInEtaBins/th2d_mass_PhiPlus_barrel-forward",
0877                                           "TkTkMassVsPhiPlusInEtaBins/th2d_mass_PhiPlus_barrel-backward",
0878                                           "TkTkMassVsPhiPlusInEtaBins/th2d_mass_PhiPlus_forward-forward",
0879                                           "TkTkMassVsPhiPlusInEtaBins/th2d_mass_PhiPlus_backward-backward",
0880                                           "TkTkMassVsPhiPlusInEtaBins/th2d_mass_PhiPlus_forward-backward"};
0881
0882   producePlots(v_meanHistos, MEtoHarvest, labels, Rlabel, useAutoLimits, false);
0883   producePlots(v_widthHistos, MEtoHarvest, labels, Rlabel, useAutoLimits, true);
0884
0885   std::vector<std::string> MEtoHarvest3D = {"th3d_mass_vs_eta_phi_plus", "th3d_mass_vs_eta_phi_minus"};
0886
0887   produceMaps(v_meanMaps, MEtoHarvest3D, labels, Rlabel, useAutoLimits, false);
0888   produceMaps(v_widthMaps, MEtoHarvest3D, labels, Rlabel, useAutoLimits, true);
0889
0890   // finally close the file
0891   for (const auto& file : sourceFiles) {
0892     file->Close();
0893   }
0894 }