test/stubs/DTMeanTimerPlotter.cc

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  \author S. Bolognesi - INFN Torino
0005  */
0006
0007 #include "DTMeanTimerPlotter.h"
0008 //#include "CalibMuon/DTCalibration/src/vDriftHistos.h"
0009
0010 #include <iostream>
0011 #include <cstdio>
0012 #include <sstream>
0013 #include <string>
0014 #include <cmath>
0015
0016 #include "TFile.h"
0017 #include "TCanvas.h"
0018 #include "TCollection.h"
0019 #include "TSystem.h"
0020 #include "TF1.h"
0021 #include "TH1D.h"
0022 #include "TLegend.h"
0023 // #include "TProfile.h"
0024
0025 using namespace std;
0026
0027 // Constructor
0028 DTMeanTimerPlotter::DTMeanTimerPlotter(TFile* file) : theFile(file), theVerbosityLevel(1), theRebinning(1), color(0) {}
0029
0030 // Destructor
0031 DTMeanTimerPlotter::~DTMeanTimerPlotter() {}
0032
0033 void DTMeanTimerPlotter::plotMeanTimer(int wheel, int station, int sector, int sl, const TString& drawOptions) {
0034   TString histoName = getHistoNameSuffix(wheel, station, sector, sl);
0035
0036   if (drawOptions.Contains("SingleDeltaT0"))
0037     plotHistos(plotSingleTMaxDeltaT0(histoName), histoName, drawOptions);
0038   else if (drawOptions.Contains("SingleFormula"))
0039     plotHistos(plotSingleTMaxFormula(histoName), histoName, drawOptions);
0040   else {
0041     plotHistos(plotTotalTMax(histoName), histoName, drawOptions);
0042     cout << "use '(fit)SingleDeltaT0' or '(fit)SingleFormula' as option to fit each TMax histo separately" << endl;
0043   }
0044 }
0045
0046 // Set the verbosity of the output: 0 = silent, 1 = info, 2 = debug
0047 void DTMeanTimerPlotter::setVerbosity(unsigned int lvl) { theVerbosityLevel = lvl; }
0048
0049 void DTMeanTimerPlotter::setRebinning(unsigned int rebin) { theRebinning = rebin; }
0050
0051 void DTMeanTimerPlotter::resetColor() { color = 0; }
0052
0053 TString DTMeanTimerPlotter::getHistoNameSuffix(int wheel, int station, int sector, int sl) {
0054   TString N = (((((TString) "TMax" + (long)wheel) + (long)station) + (long)sector) + (long)sl);
0055   return N;
0056 }
0057
0058 void DTMeanTimerPlotter::plotHistos(vector<TH1D*> hTMaxes, TString& name, const TString& drawOptions) {
0059   TLegend* leg = new TLegend(0.5, 0.6, 0.7, 0.8);
0060   ;
0061   if (!drawOptions.Contains("same")) {
0062     new TCanvas(name, "Fit of Tmax histo");
0063     hTMaxes[0]->Draw();
0064   }
0065
0066   for (vector<TH1D*>::const_iterator ith = hTMaxes.begin(); ith != hTMaxes.end(); ++ith) {
0067     color++;
0068
0069     (*ith)->Rebin(theRebinning);
0070     (*ith)->SetLineColor(color);
0071     ((*ith)->GetXaxis())->SetRangeUser(200, 600);
0072     (*ith)->Draw("same");
0073
0074     leg->AddEntry((*ith), (*ith)->GetName(), "L");
0075   }
0076
0077   if (drawOptions.Contains("fit")) {
0078     int i = 0;
0079     vector<TF1*> functions = fitTMaxes(hTMaxes);
0080     for (vector<TF1*>::const_iterator funct = functions.begin(); funct != functions.end(); ++funct) {
0081       //     color++;
0082       (*funct)->SetLineColor(hTMaxes[i]->GetLineColor());
0083       (*funct)->Draw("same");
0084       i++;
0085     }
0086   }
0087   leg->SetFillColor(0);
0088   ;
0089   leg->Draw("same");
0090   hTMaxes[0]->SetMaximum(getMaximum(hTMaxes));
0091   setRebinning(1);
0092 }
0093
0094 vector<TH1D*> DTMeanTimerPlotter::plotSingleTMaxDeltaT0(TString& name) {
0095   // Retrieve histogram sets
0096   vector<TH1D*> hTMaxes;
0097   hTMaxes.push_back((TH1D*)theFile->Get(name + "_0"));
0098   hTMaxes.push_back((TH1D*)theFile->Get(name + "_t0"));
0099   hTMaxes.push_back((TH1D*)theFile->Get(name + "_2t0"));
0100   hTMaxes.push_back((TH1D*)theFile->Get(name + "_3t0"));
0101   if (theVerbosityLevel >= 1) {
0102     cout << "NOTE: these histos are not directly used to calibrate vdrift" << endl;
0103   }
0104   return hTMaxes;
0105 }
0106
0107 vector<TH1D*> DTMeanTimerPlotter::plotSingleTMaxFormula(TString& name) {
0108   // Retrieve histogram sets
0109   vector<TH1D*> hTMaxes;
0110   hTMaxes.push_back((TH1D*)theFile->Get(name + "_Tmax123"));
0111   hTMaxes.push_back((TH1D*)theFile->Get(name + "_Tmax124_s72"));
0112   hTMaxes.push_back((TH1D*)theFile->Get(name + "_Tmax124_s78"));
0113   hTMaxes.push_back((TH1D*)theFile->Get(name + "_Tmax134_s72"));
0114   hTMaxes.push_back((TH1D*)theFile->Get(name + "_Tmax134_s78"));
0115   hTMaxes.push_back((TH1D*)theFile->Get(name + "_Tmax234"));
0116   return hTMaxes;
0117 }
0118
0119 vector<TH1D*> DTMeanTimerPlotter::plotTotalTMax(TString& name) {
0120   TH1D* hTMax;  // histograms for <T_max> calculation
0121   hTMax = ((TH1D*)theFile->Get(name + "_Tmax123"));
0122   hTMax->Add((TH1D*)theFile->Get(name + "_Tmax124_s72"));
0123   hTMax->Add((TH1D*)theFile->Get(name + "_Tmax124_s78"));
0124   hTMax->Add((TH1D*)theFile->Get(name + "_Tmax134_s72"));
0125   hTMax->Add((TH1D*)theFile->Get(name + "_Tmax134_s78"));
0126   hTMax->Add((TH1D*)theFile->Get(name + "_Tmax234"));
0127   hTMax->SetName(name);
0128   hTMax->SetTitle(name);
0129
0130   vector<TH1D*> hTMaxes;
0131   hTMaxes.push_back(hTMax);
0132   if (theVerbosityLevel >= 1) {
0133     cout << "NOTE: this histo is not directly used to calibrate vdrift" << endl;
0134   }
0135   return hTMaxes;
0136 }
0137
0138 vector<TF1*> DTMeanTimerPlotter::fitTMaxes(vector<TH1D*> hTMaxes) {
0139   vector<TF1*> functions;
0140   for (vector<TH1D*>::const_iterator ith = hTMaxes.begin(); ith != hTMaxes.end(); ++ith) {
0141     // Find distribution peak and fit range
0142     Double_t peak = ((((((*ith)->GetXaxis())->GetXmax()) - (((*ith)->GetXaxis())->GetXmin())) / (*ith)->GetNbinsX()) *
0143                      ((*ith)->GetMaximumBin())) +
0144                     (((*ith)->GetXaxis())->GetXmin());
0145     if (theVerbosityLevel >= 1)
0146       cout << "Peak " << peak << " : "
0147            << "xmax " << (((*ith)->GetXaxis())->GetXmax()) << "            xmin " << (((*ith)->GetXaxis())->GetXmin())
0148            << "            nbin " << (*ith)->GetNbinsX() << "            bin with max " << ((*ith)->GetMaximumBin())
0149            << endl;
0150     Double_t range = 2. * (*ith)->GetRMS();
0151
0152     // Fit each Tmax (*ith)gram with a Gaussian in a restricted interval
0153     TF1* rGaus = new TF1("rGaus", "gaus", peak - range, peak + range);
0154     rGaus->SetMarkerSize();  //to stop gcc complain about unused var
0155     (*ith)->Fit("rGaus", "R0");
0156     functions.push_back((*ith)->GetFunction("rGaus"));
0157
0158     // Get mean, sigma and number of entries of each histogram
0159     cout << "Histo name " << (*ith)->GetName() << ": " << endl;
0160     cout << "mean  " << (((*ith)->GetFunction("rGaus"))->GetParameter(1)) << endl;
0161     cout << "sigma " << (((*ith)->GetFunction("rGaus"))->GetParameter(2)) << endl;
0162     cout << "count " << ((*ith)->GetEntries()) << endl << endl;
0163   }
0164   return functions;
0165 }
0166
0167 double DTMeanTimerPlotter::getMaximum(vector<TH1D*> hTMaxes) {
0168   double max = -pow(10.0, 10);
0169   for (vector<TH1D*>::const_iterator ith = hTMaxes.begin(); ith != hTMaxes.end(); ++ith) {
0170     double m = (*ith)->GetMaximum(pow(10.0, 10));
0171     if (m > max)
0172       max = m;
0173   }
0174   return max;
0175 }