|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
|
/*
* See header file for a description of this class.
*
* \author S. Bolognesi - INFN Torino
*/
#include "DTMeanTimerPlotter.h"
//#include "CalibMuon/DTCalibration/src/vDriftHistos.h"
#include <iostream>
#include <cstdio>
#include <sstream>
#include <string>
#include <cmath>
#include "TFile.h"
#include "TCanvas.h"
#include "TCollection.h"
#include "TSystem.h"
#include "TF1.h"
#include "TH1D.h"
#include "TLegend.h"
// #include "TProfile.h"
using namespace std;
// Constructor
DTMeanTimerPlotter::DTMeanTimerPlotter(TFile* file) : theFile(file), theVerbosityLevel(1), theRebinning(1), color(0) {}
// Destructor
DTMeanTimerPlotter::~DTMeanTimerPlotter() {}
void DTMeanTimerPlotter::plotMeanTimer(int wheel, int station, int sector, int sl, const TString& drawOptions) {
TString histoName = getHistoNameSuffix(wheel, station, sector, sl);
if (drawOptions.Contains("SingleDeltaT0"))
plotHistos(plotSingleTMaxDeltaT0(histoName), histoName, drawOptions);
else if (drawOptions.Contains("SingleFormula"))
plotHistos(plotSingleTMaxFormula(histoName), histoName, drawOptions);
else {
plotHistos(plotTotalTMax(histoName), histoName, drawOptions);
cout << "use '(fit)SingleDeltaT0' or '(fit)SingleFormula' as option to fit each TMax histo separately" << endl;
}
}
// Set the verbosity of the output: 0 = silent, 1 = info, 2 = debug
void DTMeanTimerPlotter::setVerbosity(unsigned int lvl) { theVerbosityLevel = lvl; }
void DTMeanTimerPlotter::setRebinning(unsigned int rebin) { theRebinning = rebin; }
void DTMeanTimerPlotter::resetColor() { color = 0; }
TString DTMeanTimerPlotter::getHistoNameSuffix(int wheel, int station, int sector, int sl) {
TString N = (((((TString) "TMax" + (long)wheel) + (long)station) + (long)sector) + (long)sl);
return N;
}
void DTMeanTimerPlotter::plotHistos(vector<TH1D*> hTMaxes, TString& name, const TString& drawOptions) {
TLegend* leg = new TLegend(0.5, 0.6, 0.7, 0.8);
;
if (!drawOptions.Contains("same")) {
new TCanvas(name, "Fit of Tmax histo");
hTMaxes[0]->Draw();
}
for (vector<TH1D*>::const_iterator ith = hTMaxes.begin(); ith != hTMaxes.end(); ++ith) {
color++;
(*ith)->Rebin(theRebinning);
(*ith)->SetLineColor(color);
((*ith)->GetXaxis())->SetRangeUser(200, 600);
(*ith)->Draw("same");
leg->AddEntry((*ith), (*ith)->GetName(), "L");
}
if (drawOptions.Contains("fit")) {
int i = 0;
vector<TF1*> functions = fitTMaxes(hTMaxes);
for (vector<TF1*>::const_iterator funct = functions.begin(); funct != functions.end(); ++funct) {
// color++;
(*funct)->SetLineColor(hTMaxes[i]->GetLineColor());
(*funct)->Draw("same");
i++;
}
}
leg->SetFillColor(0);
;
leg->Draw("same");
hTMaxes[0]->SetMaximum(getMaximum(hTMaxes));
setRebinning(1);
}
vector<TH1D*> DTMeanTimerPlotter::plotSingleTMaxDeltaT0(TString& name) {
// Retrieve histogram sets
vector<TH1D*> hTMaxes;
hTMaxes.push_back((TH1D*)theFile->Get(name + "_0"));
hTMaxes.push_back((TH1D*)theFile->Get(name + "_t0"));
hTMaxes.push_back((TH1D*)theFile->Get(name + "_2t0"));
hTMaxes.push_back((TH1D*)theFile->Get(name + "_3t0"));
if (theVerbosityLevel >= 1) {
cout << "NOTE: these histos are not directly used to calibrate vdrift" << endl;
}
return hTMaxes;
}
vector<TH1D*> DTMeanTimerPlotter::plotSingleTMaxFormula(TString& name) {
// Retrieve histogram sets
vector<TH1D*> hTMaxes;
hTMaxes.push_back((TH1D*)theFile->Get(name + "_Tmax123"));
hTMaxes.push_back((TH1D*)theFile->Get(name + "_Tmax124_s72"));
hTMaxes.push_back((TH1D*)theFile->Get(name + "_Tmax124_s78"));
hTMaxes.push_back((TH1D*)theFile->Get(name + "_Tmax134_s72"));
hTMaxes.push_back((TH1D*)theFile->Get(name + "_Tmax134_s78"));
hTMaxes.push_back((TH1D*)theFile->Get(name + "_Tmax234"));
return hTMaxes;
}
vector<TH1D*> DTMeanTimerPlotter::plotTotalTMax(TString& name) {
TH1D* hTMax; // histograms for <T_max> calculation
hTMax = ((TH1D*)theFile->Get(name + "_Tmax123"));
hTMax->Add((TH1D*)theFile->Get(name + "_Tmax124_s72"));
hTMax->Add((TH1D*)theFile->Get(name + "_Tmax124_s78"));
hTMax->Add((TH1D*)theFile->Get(name + "_Tmax134_s72"));
hTMax->Add((TH1D*)theFile->Get(name + "_Tmax134_s78"));
hTMax->Add((TH1D*)theFile->Get(name + "_Tmax234"));
hTMax->SetName(name);
hTMax->SetTitle(name);
vector<TH1D*> hTMaxes;
hTMaxes.push_back(hTMax);
if (theVerbosityLevel >= 1) {
cout << "NOTE: this histo is not directly used to calibrate vdrift" << endl;
}
return hTMaxes;
}
vector<TF1*> DTMeanTimerPlotter::fitTMaxes(vector<TH1D*> hTMaxes) {
vector<TF1*> functions;
for (vector<TH1D*>::const_iterator ith = hTMaxes.begin(); ith != hTMaxes.end(); ++ith) {
// Find distribution peak and fit range
Double_t peak = ((((((*ith)->GetXaxis())->GetXmax()) - (((*ith)->GetXaxis())->GetXmin())) / (*ith)->GetNbinsX()) *
((*ith)->GetMaximumBin())) +
(((*ith)->GetXaxis())->GetXmin());
if (theVerbosityLevel >= 1)
cout << "Peak " << peak << " : "
<< "xmax " << (((*ith)->GetXaxis())->GetXmax()) << " xmin " << (((*ith)->GetXaxis())->GetXmin())
<< " nbin " << (*ith)->GetNbinsX() << " bin with max " << ((*ith)->GetMaximumBin())
<< endl;
Double_t range = 2. * (*ith)->GetRMS();
// Fit each Tmax (*ith)gram with a Gaussian in a restricted interval
TF1* rGaus = new TF1("rGaus", "gaus", peak - range, peak + range);
rGaus->SetMarkerSize(); //to stop gcc complain about unused var
(*ith)->Fit("rGaus", "R0");
functions.push_back((*ith)->GetFunction("rGaus"));
// Get mean, sigma and number of entries of each histogram
cout << "Histo name " << (*ith)->GetName() << ": " << endl;
cout << "mean " << (((*ith)->GetFunction("rGaus"))->GetParameter(1)) << endl;
cout << "sigma " << (((*ith)->GetFunction("rGaus"))->GetParameter(2)) << endl;
cout << "count " << ((*ith)->GetEntries()) << endl << endl;
}
return functions;
}
double DTMeanTimerPlotter::getMaximum(vector<TH1D*> hTMaxes) {
double max = -pow(10.0, 10);
for (vector<TH1D*>::const_iterator ith = hTMaxes.begin(); ith != hTMaxes.end(); ++ith) {
double m = (*ith)->GetMaximum(pow(10.0, 10));
if (m > max)
max = m;
}
return max;
}
|