File indexing completed on 2024-05-22 04:02:32
0001
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
0027 #include <iomanip>
0028 #include <iostream>
0029 #include <map>
0030 #include <fmt/core.h>
0031
0032
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
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 }
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
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
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;
0105
0106
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
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
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
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
0304 const Measurement1D& bias = results.getBias();
0305 meanHistos_[key]->SetBinContent(bin, bias.value());
0306 meanHistos_[key]->SetBinError(bin, bias.error());
0307
0308
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
0443 const Measurement1D& bias = results.getBias();
0444 meanMaps_[key]->SetBinContent(binY, binZ, bias.value());
0445 meanMaps_[key]->SetBinError(binY, binZ, bias.error());
0446
0447
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
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
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
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
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
0594 TLegend* infoBox = new TLegend(0.65, 0.75, 0.95, 0.90, "");
0595 infoBox->SetShadowColor(0);
0596 infoBox->SetFillColor(kWhite);
0597 infoBox->SetTextSize(0.035);
0598
0599
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
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
0641 if (!useAutoLimits) {
0642 if (isWidth) {
0643
0644 histoMap.at(var)->GetYaxis()->SetRangeUser(0.5, 2.85);
0645 } else {
0646
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
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
0673 size_t pos = var.find('/');
0674 std::string outputName{var};
0675
0676
0677 if (pos != std::string::npos) {
0678
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
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
0711 TLegend* infoBox = new TLegend(0.65, 0.75, 0.95, 0.90, "");
0712 infoBox->SetShadowColor(0);
0713 infoBox->SetFillColor(kWhite);
0714 infoBox->SetTextSize(0.035);
0715
0716
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
0752 if (!useAutoLimits) {
0753 if (isWidth) {
0754
0755 histoMap.at(var)->GetZaxis()->SetRangeUser(0.5, 2.85);
0756 } else {
0757
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
0771 if (var == MEtoHarvest[0]) {
0772 infoBox->AddEntry(histoMap.at(var), labels[count], "LP");
0773 }
0774
0775
0776 CMS_lumi(dynamic_cast<TPad*>(gPad), 0, 3, labels[count]);
0777 count++;
0778 }
0779
0780
0781 size_t pos = var.find('/');
0782 std::string outputName{var};
0783
0784
0785 if (pos != std::string::npos) {
0786
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
0821 std::vector<diMuonMassBias::histoMap> v_meanHistos;
0822 std::vector<diMuonMassBias::histoMap> v_widthHistos;
0823
0824
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
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
0891 for (const auto& file : sourceFiles) {
0892 file->Close();
0893 }
0894 }