MtdValidation/macros/Pt_residuals_fit.C

0001 // -*- C -*-
0002 //
0003 // 'Resolution in track Pt' macro
0004 //
0005 // \author 12/2023 - Raffaele Delli Gatti
0006
0007 #include "Riostream.h"
0008 #include "TFile.h"
0009 #include "TDirectoryFile.h"
0010 #include "TTree.h"
0011 #include "TCanvas.h"
0012 #include "TH1D.h"
0013 #include "TStyle.h"
0014 #include "TLegend.h"
0015 #include "TLatex.h"
0016 #include "TString.h"
0017 #include "RooRealVar.h"
0018 #include "RooCrystalBall.h"
0019 #include "RooAddPdf.h"
0020 #include "RooDataHist.h"
0021 #include "RooPlot.h"
0022 #include "RooFitResult.h"
0023 using namespace std;
0024 using namespace RooFit;
0025
0026 // Funtion for fitting to data with roofit library (defined below)
0027 // ---------------------------------------------------------------
0028 void fit_to_data(TH1D* histogram, TString name_file);
0029
0030 // Main function
0031 //--------------
0032 void residuals_fit() {
0033   // Open the root file to read
0034   // --------------------------
0035
0036   TFile* file_DQM = TFile::Open("./DQM_V0001_R000000001__Global__CMSSW_X_Y_Z__RECO.root", "READ");
0037   // and take its directories
0038   TDirectoryFile* dir_DQMData = (TDirectoryFile*)file_DQM->Get("DQMData");
0039   if (!dir_DQMData)
0040     cout << "Cannot find dir_DQMData" << endl;
0041   TDirectoryFile* dir_Run1 = (TDirectoryFile*)dir_DQMData->Get("Run 1");
0042   if (!dir_Run1)
0043     cout << "Cannot find dir_Run1" << endl;
0044   TDirectoryFile* dir_MTD = (TDirectoryFile*)dir_Run1->Get("MTD");
0045   if (!dir_MTD)
0046     cout << "Cannot find dir_MTD" << endl;
0047   TDirectoryFile* dir_Runsum = (TDirectoryFile*)dir_MTD->Get("Run summary");
0048   if (!dir_Runsum)
0049     cout << "Cannot find dir_Runsum" << endl;
0050   TDirectoryFile* dir_Tracks = (TDirectoryFile*)dir_Runsum->Get("Tracks");
0051   if (!dir_Tracks)
0052     cout << "Cannot find dir_Tracks" << endl;
0053
0054   // Take the trees with the method Get()
0055   TH1D* h_TrackMatchedTPBTLPtRatioGen = (TH1D*)dir_Tracks->Get("TrackMatchedTPBTLPtRatioGen");
0056   TH1D* h_TrackMatchedTPBTLPtRatioMtd = (TH1D*)dir_Tracks->Get("TrackMatchedTPBTLPtRatioMtd");
0057   TH1D* h_TrackMatchedTPBTLPtResMtd = (TH1D*)dir_Tracks->Get("TrackMatchedTPBTLPtResMtd");
0058
0059   TH1D* h_TrackMatchedTPETLPtRatioGen = (TH1D*)dir_Tracks->Get("TrackMatchedTPETLPtRatioGen");
0060   TH1D* h_TrackMatchedTPETLPtRatioMtd = (TH1D*)dir_Tracks->Get("TrackMatchedTPETLPtRatioMtd");
0061   TH1D* h_TrackMatchedTPETLPtResMtd = (TH1D*)dir_Tracks->Get("TrackMatchedTPETLPtResMtd");
0062
0063   TH1D* h_TrackMatchedTPETL2PtRatioGen = (TH1D*)dir_Tracks->Get("TrackMatchedTPETL2PtRatioGen");
0064   TH1D* h_TrackMatchedTPETL2PtRatioMtd = (TH1D*)dir_Tracks->Get("TrackMatchedTPETL2PtRatioMtd");
0065   TH1D* h_TrackMatchedTPETL2PtResMtd = (TH1D*)dir_Tracks->Get("TrackMatchedTPETL2PtResMtd");
0066
0067   // Fit to data with the function fit_to_data
0068   //------------------------------------------
0069   fit_to_data(h_TrackMatchedTPBTLPtRatioGen, "BTLPtRatioGen_fit.pdf");
0070   fit_to_data(h_TrackMatchedTPETLPtRatioGen, "ETLPtRatioGen_fit.pdf");
0071   fit_to_data(h_TrackMatchedTPETL2PtRatioGen, "ETL2PtRatioGen_fit.pdf");
0072
0073   fit_to_data(h_TrackMatchedTPBTLPtRatioMtd, "BTLPtRatioMtd_fit.pdf");
0074   fit_to_data(h_TrackMatchedTPETLPtRatioMtd, "ETLPtRatioMtd_fit.pdf");
0075   fit_to_data(h_TrackMatchedTPETL2PtRatioMtd, "ETL2PtRatioMtd_fit.pdf");
0076
0077   fit_to_data(h_TrackMatchedTPBTLPtResMtd, "BTLPtResMtd_fit.pdf");
0078   fit_to_data(h_TrackMatchedTPETLPtResMtd, "ETLPtResMtd_fit.pdf");
0079   fit_to_data(h_TrackMatchedTPETL2PtResMtd, "ETL2PtResMtd_fit.pdf");
0080
0081   return;
0082 }
0083
0084 void fit_to_data(TH1D* h_TrackMatchedTP, TString str) {
0085   // Fit to data using roofit library
0086   // --------------------------------
0087
0088   Double_t bin_width = h_TrackMatchedTP->GetBinWidth(0);
0089   Double_t range_min = h_TrackMatchedTP->GetXaxis()->GetBinLowEdge(1);
0090   Double_t range_max =
0091       h_TrackMatchedTP->GetXaxis()->GetBinLowEdge(h_TrackMatchedTP->GetNbinsX()) + h_TrackMatchedTP->GetBinWidth(0);
0092
0093   // Observable
0094   RooRealVar x_res("x res", "", range_min, range_max);
0095
0096   // Import data from histogram
0097   RooDataHist* h_ = new RooDataHist("h_", "h_", x_res, h_TrackMatchedTP);
0098
0099   // Parameters
0100   RooRealVar mean("mean", "mean", h_TrackMatchedTP->GetMean(), range_min, range_max);
0101   RooRealVar sigmaL("sigmaL", "sigmaL", 0.05, 0., 5.);
0102   RooRealVar sigmaR("sigmaR", "sigmaR", 0.05, 0., 5.);
0103   RooRealVar alphaL("alphaL", "alphaL", 1., 0., 5.);
0104   RooRealVar alphaR("alphaR", "alphaR", 1., 0., 5.);
0105   RooRealVar nL("NL", "NL", 5., 0., 100.);
0106   RooRealVar nR("NR", "NR", 5., 0., 100.);
0107
0108   // Build a double sided crystall ball PDF
0109   RooCrystalBall* pdf = new RooCrystalBall("pdf", "pdf", x_res, mean, sigmaL, sigmaR, alphaL, nL, alphaR, nR);
0110   // Construct a signal PDF
0111   RooRealVar nsig("nsig", "#signal events", h_TrackMatchedTP->GetEntries(), 0., h_TrackMatchedTP->GetEntries() * 2);
0112   RooAddPdf* model = new RooAddPdf("model", "model", {*pdf}, {nsig});
0113
0114   // The PDF fit to that data set using an un-binned maximum likelihood fit
0115   // Then the data are visualized with the PDF overlaid
0116
0117   // Perform extended ML fit of PDF to data
0118   model->fitTo(*h_);
0119
0120   // Retrieve values from the fit
0121   Double_t mean_fit = mean.getVal();
0122   Double_t err_mean = mean.getError();
0123   Double_t sigmaR_fit = sigmaR.getVal();
0124   Double_t err_sigmaR = sigmaR.getError();
0125   Double_t sigmaL_fit = sigmaL.getVal();
0126   Double_t err_sigmaL = sigmaL.getError();
0127
0128   // Compute resolution as half-width of the interval containing 68% of all entries (including overflows), centered around the MPV of the residuals
0129   // ----------------------------------------------------------------------------------------------------------------------------------------------
0130
0131   Double_t res = 0;
0132   Double_t min = 0;
0133   double_t integral = 0;
0134   for (Int_t i = 1; i < h_TrackMatchedTP->GetNbinsX() / 2; i++) {
0135     Int_t bin_mean = (mean_fit - range_min) / bin_width;
0136     double_t int_norm = h_TrackMatchedTP->Integral(bin_mean - i, bin_mean + i) / h_TrackMatchedTP->Integral();
0137     if (int_norm - 0.68 < min)
0138       res = i * bin_width;
0139   }
0140   cout << "Resolution = " << res << " +- " << bin_width << endl;
0141
0142   // Create a RooPlot to draw on
0143   //----------------------------
0144   // We don't manage the memory of the returned pointer
0145   // Instead we let it leak such that the plot still exists at the end of the macro and we can take a look at it
0146   RooPlot* xresframe = x_res.frame();
0147
0148   // Plot data and PDF overlaid
0149   h_->plotOn(xresframe, MarkerSize(0.8), Name("histogram"));
0150   model->plotOn(xresframe, LineColor(kBlue), LineWidth(3), Name("model"));
0151   // In the previous lines, name is needed for the legend
0152
0153   auto legend_res = new TLegend(0.65, 0.8, 0.8, 0.9);
0154   gStyle->SetLegendTextSize(0.033);
0155   TLatex* header = new TLatex();
0156   header->SetTextSize(0.035);
0157
0158   TCanvas* c1 = new TCanvas("c1", "c1", 600, 500);
0159   c1->cd();
0160   xresframe->GetXaxis()->SetTitle(h_TrackMatchedTP->GetXaxis()->GetTitle());
0161   xresframe->Draw();
0162   if (str.Contains("BTL") && str.Contains("Mtd"))
0163     legend_res->AddEntry("histogram", "BTL", "PLerr");
0164   else if (str.Contains("BTL") && str.Contains("Gen"))
0165     legend_res->AddEntry("histogram", "GenTrack (barrel)", "PLerr");
0166   else if (str.Contains("ETLPt") && str.Contains("Mtd"))
0167     legend_res->AddEntry("histogram", "ETL (one hit)", "PLerr");
0168   else if (str.Contains("ETLPt") && str.Contains("Gen"))
0169     legend_res->AddEntry("histogram", "GenTrack (end-caps)", "PLerr");
0170   else if (str.Contains("ETL2") && str.Contains("Mtd"))
0171     legend_res->AddEntry("histogram", "ETL (2 hits)", "PLerr");
0172   else if (str.Contains("ETL2") && str.Contains("Gen"))
0173     legend_res->AddEntry("histogram", "GenTrack (end-caps)", "PLerr");
0174   legend_res->AddEntry("model", "DSCB Fit", "L");
0175   legend_res->Draw("same");
0176   header->DrawLatexNDC(0.12, 0.96, "MC Simulation");
0177   header->DrawLatexNDC(0.81, 0.96, "#sqrt{s} = 14 TeV");
0178   header->DrawLatexNDC(0.66, 0.75, TString::Format("#mu = %.4f #pm %.4f", mean_fit, err_mean));
0179   header->DrawLatexNDC(0.66, 0.71, TString::Format("#sigma_{R} = %.4f #pm %.4f", sigmaR_fit, err_sigmaR));
0180   header->DrawLatexNDC(0.66, 0.67, TString::Format("#sigma_{L} = %.4f #pm %.4f", sigmaL_fit, err_sigmaL));
0181   if (str.Contains("Ratio"))
0182     header->DrawLatexNDC(0.66, 0.63, TString::Format("#sigma (0.68) = %.3f #pm %.3f", res, bin_width));
0183   if (str.Contains("Ratio"))
0184     header->DrawLatexNDC(0.66, 0.59, TString::Format("#chi^{2}/ndf = %.2f", xresframe->chiSquare()));
0185   if (str.Contains("Res"))
0186     header->DrawLatexNDC(0.66, 0.63, TString::Format("#chi^{2}/ndf = %.2f", xresframe->chiSquare()));
0187   c1->Print(str);
0188 }