MtdValidation/macros/makeTimeResPlots.C

0001 // =============================================================================
0002 //
0003 //  This macro reads the histograms of time residuals of the BTL and ETL
0004 //  UncalibratedRecHits from the DQM root file and produces time resolution
0005 //  plots vs the hit amplitude and the hit |eta|.
0006 //
0007 //  Usage:
0008 //     root -l make_resPlot.C\(\"DQM_filename.root\"\)
0009 //
0010 //  NB: In order to have the UncalibratedRecHits histograms filled,
0011 //      the MTD validation has to be run the flags:
0012 //
0013 //      process.btlLocalReco.UncalibRecHitsPlots = cms.bool(True);
0014 //      process.etlLocalReco.UncalibRecHitsPlots = cms.bool(True);
0015 //
0016 // =============================================================================
0017
0018 #include <iostream>
0019 #include "TStyle.h"
0020 #include "TFile.h"
0021 #include "TString.h"
0022 #include "TCanvas.h"
0023 #include "TH1F.h"
0024 #include "TFitResult.h"
0025 #include "TLegend.h"
0026
0027 // =============================================================================
0028 //  configuration
0029
0030 const Float_t minBinContent = 20.;
0031
0032 // --- BTL
0033
0034 const Int_t nBinsQ_BTL = 20;
0035 const Float_t binWidthQ_BTL = 30.;  // [pC]
0036 const Int_t nBinsQEta_BTL = 3;
0037 const Float_t binsQEta_BTL[4] = {0., 0.65, 1.15, 1.55};
0038
0039 const Int_t nBinsEta_BTL = 31;
0040 const Float_t binWidthEta_BTL = 0.05;
0041 const Int_t nBinsEtaQ_BTL = 7;
0042 const Float_t binsEtaQ_BTL[8] = {0., 30., 60., 90., 120., 150., 360., 600.};
0043
0044 // --- ETL
0045
0046 const Int_t nBinsQ_ETL = 20;
0047 const Float_t binWidthQ_ETL = 1.3;  // [MIP]
0048
0049 const Int_t nBinsEta_ETL = 26;
0050 const Float_t binWidthEta_ETL = 0.05;
0051 const Float_t etaMin_ETL = 1.65;
0052
0053 // =============================================================================
0054
0055 TCanvas* c[10];
0056
0057 TH1F* h_TimeResQ_BTL[nBinsQ_BTL];
0058 TH1F* h_TimeResQEta_BTL[nBinsQ_BTL][nBinsQEta_BTL];
0059 TH1F* h_TimeResEta_BTL[nBinsEta_BTL];
0060 TH1F* h_TimeResEtaQ_BTL[nBinsEta_BTL][nBinsEtaQ_BTL];
0061
0062 TH1F* g_TimeResQ_BTL;
0063 TH1F* g_TimeResQEta_BTL[nBinsQEta_BTL];
0064 TH1F* g_TimeResEta_BTL;
0065 TH1F* g_TimeResEtaQ_BTL[nBinsEtaQ_BTL];
0066
0067 TH1F* h_TimeResQ_ETL[2][nBinsQ_ETL];
0068 TH1F* h_TimeResEta_ETL[2][nBinsEta_ETL];
0069
0070 TH1F* g_TimeResQ_ETL[2];
0071 TH1F* g_TimeResEta_ETL[2];
0072
0073 void makeTimeResPlots(const TString DQMfilename = "DQM_V0001_UNKNOWN_R000000001.root") {
0074   gStyle->SetOptStat(kFALSE);
0075
0076   // --- Histograms booking
0077
0078   g_TimeResQ_BTL = new TH1F("g_TimeResQ_BTL", "BTL time resolution vs Q", nBinsQ_BTL, 0., nBinsQ_BTL * binWidthQ_BTL);
0079
0080   for (Int_t ih = 0; ih < nBinsQEta_BTL; ++ih) {
0081     TString hname = Form("g_TimeResQEta_BTL_%d", ih);
0082     TString htitle =
0083         Form("BTL time resolution vs Q (%4.2f < |#eta_{hit}| < %4.2f)", binsQEta_BTL[ih], binsQEta_BTL[ih + 1]);
0084     g_TimeResQEta_BTL[ih] = new TH1F(hname, htitle, nBinsQ_BTL, 0., nBinsQ_BTL * binWidthQ_BTL);
0085   }
0086
0087   g_TimeResEta_BTL = new TH1F(
0088       "g_TimeResEta_BTL", "BTL time resolution vs |#eta_{hit}|", nBinsEta_BTL, 0., nBinsEta_BTL * binWidthEta_BTL);
0089
0090   for (Int_t ih = 0; ih < nBinsEtaQ_BTL; ++ih) {
0091     TString hname = Form("g_TimeResEtaQ_BTL_%d", ih);
0092     TString htitle =
0093         Form("BTL time resolution vs |#eta_{hit}| (%4.2f < Q < %4.2f)", binsEtaQ_BTL[ih], binsEtaQ_BTL[ih + 1]);
0094     g_TimeResEtaQ_BTL[ih] = new TH1F(hname, htitle, nBinsEta_BTL, 0., nBinsEta_BTL * binWidthEta_BTL);
0095   }
0096
0097   for (Int_t iside = 0; iside < 2; ++iside) {
0098     TString hname = Form("g_TimeResQ_ETL_%d", iside);
0099     g_TimeResQ_ETL[iside] =
0100         new TH1F(hname, "ETL time resolution vs amplitude", nBinsQ_ETL, 0., nBinsQ_ETL * binWidthQ_ETL);
0101
0102     hname = Form("g_TimeResEta_ETL_%d", iside);
0103     g_TimeResEta_ETL[iside] = new TH1F(hname,
0104                                        "ETL time resolution vs |#eta_{hit}|",
0105                                        nBinsEta_ETL,
0106                                        etaMin_ETL,
0107                                        etaMin_ETL + nBinsEta_ETL * binWidthEta_ETL);
0108   }
0109
0110   std::cout << " Processing file " << DQMfilename.Data() << " ... " << std::endl;
0111   TFile* input_file = new TFile(DQMfilename.Data());
0112
0113   // ---------------------------------------------------------------------------
0114   // BTL time resolution vs Q
0115
0116   for (Int_t iq = 0; iq < nBinsQ_BTL; ++iq) {
0117     TString hname = Form("DQMData/Run 1/MTD/Run summary/BTL/LocalReco/TimeResQ_%d", iq);
0118     h_TimeResQ_BTL[iq] = (TH1F*)input_file->Get(hname);
0119
0120     if (h_TimeResQ_BTL[iq]->GetEntries() < minBinContent)
0121       continue;
0122
0123     Float_t low_limit = h_TimeResQ_BTL[iq]->GetMean() - 2. * h_TimeResQ_BTL[iq]->GetRMS();
0124     Float_t up_limit = h_TimeResQ_BTL[iq]->GetMean() + 3 * h_TimeResQ_BTL[iq]->GetRMS();
0125
0126     TFitResultPtr fit_res = h_TimeResQ_BTL[iq]->Fit("gaus", "SERQ0", "", low_limit, up_limit);
0127
0128     if (fit_res->Status() != 0) {
0129       std::cout << " *** ERROR: fit failed for the histogram" << h_TimeResQ_BTL[iq]->GetName() << std::endl;
0130       continue;
0131     }
0132
0133     g_TimeResQ_BTL->SetBinContent(iq + 1, fit_res->GetParams()[2]);
0134     g_TimeResQ_BTL->SetBinError(iq + 1, fit_res->GetErrors()[2]);
0135
0136   }  // iq loop
0137
0138   // ---------------------------------------------------------------------------
0139   // BTL time resolution vs Q in bind of |eta|
0140
0141   for (Int_t iq = 0; iq < nBinsQ_BTL; ++iq) {
0142     for (Int_t ieta = 0; ieta < nBinsQEta_BTL; ++ieta) {
0143       TString hname = Form("DQMData/Run 1/MTD/Run summary/BTL/LocalReco/TimeResQvsEta_%d_%d", iq, ieta);
0144       h_TimeResQEta_BTL[iq][ieta] = (TH1F*)input_file->Get(hname);
0145
0146       if (h_TimeResQEta_BTL[iq][ieta]->GetEntries() < minBinContent)
0147         continue;
0148
0149       Float_t low_limit = h_TimeResQEta_BTL[iq][ieta]->GetMean() - 2. * h_TimeResQEta_BTL[iq][ieta]->GetRMS();
0150       Float_t up_limit = h_TimeResQEta_BTL[iq][ieta]->GetMean() + 3 * h_TimeResQEta_BTL[iq][ieta]->GetRMS();
0151
0152       TFitResultPtr fit_res = h_TimeResQEta_BTL[iq][ieta]->Fit("gaus", "SERQ0", "", low_limit, up_limit);
0153
0154       if (fit_res->Status() != 0) {
0155         std::cout << " *** ERROR: fit failed for the histogram" << h_TimeResQEta_BTL[iq][ieta]->GetName() << std::endl;
0156         continue;
0157       }
0158
0159       g_TimeResQEta_BTL[ieta]->SetBinContent(iq + 1, fit_res->GetParams()[2]);
0160       g_TimeResQEta_BTL[ieta]->SetBinError(iq + 1, fit_res->GetErrors()[2]);
0161
0162     }  // ieta loop
0163
0164   }  // iq loop
0165
0166   // ---------------------------------------------------------------------------
0167   // BTL time resolution vs |eta|
0168
0169   for (Int_t ieta = 0; ieta < nBinsEta_BTL; ++ieta) {
0170     TString hname = Form("DQMData/Run 1/MTD/Run summary/BTL/LocalReco/TimeResEta_%d", ieta);
0171     h_TimeResEta_BTL[ieta] = (TH1F*)input_file->Get(hname);
0172
0173     if (h_TimeResEta_BTL[ieta]->GetEntries() < minBinContent)
0174       continue;
0175
0176     Float_t low_limit = h_TimeResEta_BTL[ieta]->GetMean() - 2. * h_TimeResEta_BTL[ieta]->GetRMS();
0177     Float_t up_limit = h_TimeResEta_BTL[ieta]->GetMean() + 3 * h_TimeResEta_BTL[ieta]->GetRMS();
0178
0179     TFitResultPtr fit_res = h_TimeResEta_BTL[ieta]->Fit("gaus", "SERQ0", "", low_limit, up_limit);
0180
0181     if (fit_res->Status() != 0) {
0182       std::cout << " *** ERROR: fit failed for the histogram" << h_TimeResEta_BTL[ieta]->GetName() << std::endl;
0183       continue;
0184     }
0185
0186     g_TimeResEta_BTL->SetBinContent(ieta + 1, fit_res->GetParams()[2]);
0187     g_TimeResEta_BTL->SetBinError(ieta + 1, fit_res->GetErrors()[2]);
0188
0189   }  // ieta loop
0190
0191   // ---------------------------------------------------------------------------
0192   // BTL time resolution vs |eta| in bind of Q
0193
0194   for (Int_t ieta = 0; ieta < nBinsEta_BTL; ++ieta) {
0195     for (Int_t iq = 0; iq < nBinsEtaQ_BTL; ++iq) {
0196       TString hname = Form("DQMData/Run 1/MTD/Run summary/BTL/LocalReco/TimeResEtavsQ_%d_%d", ieta, iq);
0197       h_TimeResEtaQ_BTL[ieta][iq] = (TH1F*)input_file->Get(hname);
0198
0199       if (h_TimeResEtaQ_BTL[ieta][iq]->GetEntries() < minBinContent)
0200         continue;
0201
0202       Float_t low_limit = h_TimeResEtaQ_BTL[ieta][iq]->GetMean() - 2. * h_TimeResEtaQ_BTL[ieta][iq]->GetRMS();
0203       Float_t up_limit = h_TimeResEtaQ_BTL[ieta][iq]->GetMean() + 3 * h_TimeResEtaQ_BTL[ieta][iq]->GetRMS();
0204
0205       TFitResultPtr fit_res = h_TimeResEtaQ_BTL[ieta][iq]->Fit("gaus", "SERQ0", "", low_limit, up_limit);
0206
0207       if (fit_res->Status() != 0) {
0208         std::cout << " *** ERROR: fit failed for the histogram" << h_TimeResEtaQ_BTL[ieta][iq]->GetName() << std::endl;
0209         continue;
0210       }
0211
0212       g_TimeResEtaQ_BTL[iq]->SetBinContent(ieta + 1, fit_res->GetParams()[2]);
0213       g_TimeResEtaQ_BTL[iq]->SetBinError(ieta + 1, fit_res->GetErrors()[2]);
0214
0215     }  // iq loop
0216
0217   }  // ieta loop
0218
0219   // ---------------------------------------------------------------------------
0220   // ETL time resolution vs amplitude and |eta|
0221
0222   for (Int_t iside = 0; iside < 2; ++iside) {
0223     for (Int_t iq = 0; iq < nBinsQ_ETL; ++iq) {
0224       TString hname = Form("DQMData/Run 1/MTD/Run summary/ETL/LocalReco/TimeResQ_%d_%d", iside, iq);
0225       h_TimeResQ_ETL[iside][iq] = (TH1F*)input_file->Get(hname);
0226
0227       if (h_TimeResQ_ETL[iside][iq]->GetEntries() < minBinContent)
0228         continue;
0229
0230       Float_t low_limit = h_TimeResQ_ETL[iside][iq]->GetMean() - 2. * h_TimeResQ_ETL[iside][iq]->GetRMS();
0231       Float_t up_limit = h_TimeResQ_ETL[iside][iq]->GetMean() + 3 * h_TimeResQ_ETL[iside][iq]->GetRMS();
0232
0233       TFitResultPtr fit_res = h_TimeResQ_ETL[iside][iq]->Fit("gaus", "SERQ0", "", low_limit, up_limit);
0234
0235       if (fit_res->Status() != 0) {
0236         std::cout << " *** ERROR: fit failed for the histogram" << h_TimeResQ_ETL[iside][iq]->GetName() << std::endl;
0237         continue;
0238       }
0239
0240       g_TimeResQ_ETL[iside]->SetBinContent(iq + 1, fit_res->GetParams()[2]);
0241       g_TimeResQ_ETL[iside]->SetBinError(iq + 1, fit_res->GetErrors()[2]);
0242
0243     }  // iq loop
0244
0245     for (Int_t ieta = 0; ieta < nBinsEta_ETL; ++ieta) {
0246       TString hname = Form("DQMData/Run 1/MTD/Run summary/ETL/LocalReco/TimeResEta_%d_%d", iside, ieta);
0247       h_TimeResEta_ETL[iside][ieta] = (TH1F*)input_file->Get(hname);
0248
0249       if (h_TimeResEta_ETL[iside][ieta]->GetEntries() < minBinContent)
0250         continue;
0251
0252       Float_t low_limit = h_TimeResEta_ETL[iside][ieta]->GetMean() - 2. * h_TimeResEta_ETL[iside][ieta]->GetRMS();
0253       Float_t up_limit = h_TimeResEta_ETL[iside][ieta]->GetMean() + 3 * h_TimeResEta_ETL[iside][ieta]->GetRMS();
0254
0255       TFitResultPtr fit_res = h_TimeResEta_ETL[iside][ieta]->Fit("gaus", "SERQ0", "", low_limit, up_limit);
0256
0257       if (fit_res->Status() != 0) {
0258         std::cout << " *** ERROR: fit failed for the histogram" << h_TimeResEta_ETL[iside][ieta]->GetName()
0259                   << std::endl;
0260         continue;
0261       }
0262
0263       g_TimeResEta_ETL[iside]->SetBinContent(ieta + 1, fit_res->GetParams()[2]);
0264       g_TimeResEta_ETL[iside]->SetBinError(ieta + 1, fit_res->GetErrors()[2]);
0265
0266     }  // ieta loop
0267
0268   }  // iside loop
0269
0270   // =============================================================================
0271   // Draw the histograms
0272
0273   // --- BTL
0274
0275   c[0] = new TCanvas("c_0", "BTL time resolution vs Q");
0276   c[0]->SetGridy(kTRUE);
0277   g_TimeResQ_BTL->SetTitle("BTL UncalibratedRecHits");
0278   g_TimeResQ_BTL->SetMarkerStyle(20);
0279   g_TimeResQ_BTL->SetMarkerColor(4);
0280   g_TimeResQ_BTL->SetXTitle("hit amplitude [pC]");
0281   g_TimeResQ_BTL->SetYTitle("#sigma_{T} [ns]");
0282   g_TimeResQ_BTL->GetYaxis()->SetRangeUser(0., 0.075);
0283   g_TimeResQ_BTL->Draw("EP");
0284
0285   c[1] = new TCanvas("c_1", "BTL time resolution vs Q");
0286   c[1]->SetGridy(kTRUE);
0287   g_TimeResQEta_BTL[0]->SetTitle("BTL UncalibratedRecHits");
0288   g_TimeResQEta_BTL[0]->SetXTitle("hit amplitude [pC]");
0289   g_TimeResQEta_BTL[0]->SetYTitle("#sigma_{T} [ns]");
0290   g_TimeResQEta_BTL[0]->GetYaxis()->SetRangeUser(0., 0.075);
0291   g_TimeResQEta_BTL[0]->SetMarkerStyle(20);
0292   g_TimeResQEta_BTL[0]->SetMarkerColor(1);
0293   g_TimeResQEta_BTL[0]->Draw("EP");
0294   g_TimeResQEta_BTL[1]->SetMarkerStyle(20);
0295   g_TimeResQEta_BTL[1]->SetMarkerColor(4);
0296   g_TimeResQEta_BTL[1]->Draw("EPSAME");
0297   g_TimeResQEta_BTL[2]->SetMarkerStyle(20);
0298   g_TimeResQEta_BTL[2]->SetMarkerColor(2);
0299   g_TimeResQEta_BTL[2]->Draw("EPSAME");
0300
0301   TLegend* legend_1 = new TLegend(0.510, 0.634, 0.862, 0.847);
0302   legend_1->SetBorderSize(0.);
0303   legend_1->SetFillStyle(0);
0304   for (Int_t ih = 0; ih < nBinsQEta_BTL; ++ih)
0305     legend_1->AddEntry(
0306         g_TimeResQEta_BTL[ih], Form("%4.2f #leq |#eta_{hit}| < %4.2f", binsQEta_BTL[ih], binsQEta_BTL[ih + 1]), "LP");
0307   legend_1->Draw();
0308
0309   c[2] = new TCanvas("c_2", "BTL time resolution vs eta");
0310   c[2]->SetGridy(kTRUE);
0311   g_TimeResEta_BTL->SetTitle("BTL UncalibratedRecHits");
0312   g_TimeResEta_BTL->SetMarkerStyle(20);
0313   g_TimeResEta_BTL->SetMarkerColor(4);
0314   g_TimeResEta_BTL->SetXTitle("|#eta_{hit}|");
0315   g_TimeResEta_BTL->SetYTitle("#sigma_{T} [ns]");
0316   g_TimeResEta_BTL->GetYaxis()->SetRangeUser(0., 0.075);
0317   g_TimeResEta_BTL->Draw("EP");
0318
0319   c[3] = new TCanvas("c_3", "BTL time resolution vs eta");
0320   c[3]->SetGridy(kTRUE);
0321   g_TimeResEtaQ_BTL[1]->SetTitle("BTL UncalibratedRecHits");
0322   g_TimeResEtaQ_BTL[1]->SetXTitle("|#eta_{hit}|");
0323   g_TimeResEtaQ_BTL[1]->SetYTitle("#sigma_{T} [ns]");
0324   g_TimeResEtaQ_BTL[1]->GetYaxis()->SetRangeUser(0., 0.075);
0325   g_TimeResEtaQ_BTL[1]->SetMarkerStyle(20);
0326   g_TimeResEtaQ_BTL[1]->SetMarkerColor(1);
0327   g_TimeResEtaQ_BTL[1]->Draw("EP");
0328
0329   TLegend* legend_3 = new TLegend(0.649, 0.124, 0.862, 0.328);
0330   legend_3->SetBorderSize(0.);
0331   legend_3->SetFillStyle(0);
0332   for (Int_t ih = 1; ih < nBinsEtaQ_BTL; ++ih) {
0333     if (ih > 1) {
0334       g_TimeResEtaQ_BTL[ih]->SetMarkerStyle(20);
0335       g_TimeResEtaQ_BTL[ih]->SetMarkerColor(ih);
0336       g_TimeResEtaQ_BTL[ih]->Draw("EPSAME");
0337     }
0338
0339     legend_3->AddEntry(
0340         g_TimeResEtaQ_BTL[ih], Form("%4.0f #leq Q_{hit} < %4.0f pC", binsEtaQ_BTL[ih], binsEtaQ_BTL[ih + 1]), "LP");
0341   }
0342
0343   legend_3->Draw();
0344
0345   // --- ETL
0346
0347   c[4] = new TCanvas("c_4", "ETL time resolution vs amplitude");
0348   c[4]->SetGridy(kTRUE);
0349   g_TimeResQ_ETL[0]->SetTitle("ETL UncalibratedRecHits");
0350   g_TimeResQ_ETL[0]->SetXTitle("hit amplitude [MIP]");
0351   g_TimeResQ_ETL[0]->SetYTitle("#sigma_{T} [ns]");
0352   g_TimeResQ_ETL[0]->GetYaxis()->SetRangeUser(0., 0.075);
0353   g_TimeResQ_ETL[0]->SetMarkerStyle(20);
0354   g_TimeResQ_ETL[0]->SetMarkerColor(4);
0355   g_TimeResQ_ETL[0]->Draw("EP");
0356   g_TimeResQ_ETL[1]->SetMarkerStyle(20);
0357   g_TimeResQ_ETL[1]->SetMarkerColor(2);
0358   g_TimeResQ_ETL[1]->Draw("EPSAME");
0359
0360   TLegend* legend_4 = new TLegend(0.673, 0.655, 0.872, 0.819);
0361   legend_4->SetBorderSize(0.);
0362   legend_4->SetFillStyle(0);
0363   legend_4->AddEntry(g_TimeResQ_ETL[0], "ETL-", "LP");
0364   legend_4->AddEntry(g_TimeResQ_ETL[1], "ETL+", "LP");
0365   legend_4->DrawClone();
0366
0367   c[5] = new TCanvas("c_5", "ETL time resolution vs |eta|");
0368   c[5]->SetGridy(kTRUE);
0369   g_TimeResEta_ETL[0]->SetTitle("ETL UncalibratedRecHits");
0370   g_TimeResEta_ETL[0]->SetXTitle("|#eta_{hit}|");
0371   g_TimeResEta_ETL[0]->SetYTitle("#sigma_{T} [ns]");
0372   g_TimeResEta_ETL[0]->GetYaxis()->SetRangeUser(0., 0.075);
0373   g_TimeResEta_ETL[0]->SetMarkerStyle(20);
0374   g_TimeResEta_ETL[0]->SetMarkerColor(4);
0375   g_TimeResEta_ETL[0]->Draw("EP");
0376   g_TimeResEta_ETL[1]->SetMarkerStyle(20);
0377   g_TimeResEta_ETL[1]->SetMarkerColor(2);
0378   g_TimeResEta_ETL[1]->Draw("EPSAME");
0379
0380   legend_4->DrawClone();
0381 }