-
Notifications
You must be signed in to change notification settings - Fork 1
/
energy_resolution.cxx
executable file
·147 lines (112 loc) · 4.6 KB
/
energy_resolution.cxx
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
#include "macros.h"
TString fFileName_spectrum;
std::vector<double> ffitrange_low;
std::vector<double> ffitrange_high;
std::vector<double> famp_st;
std::vector<double> fmean_st;
std::vector<double> fsigma_st;
std::vector<double> fconst_st;
std::vector<double> fslope_st;
// ----------------------------------------------------
// read parameters from user specified file
// ----------------------------------------------------
int read_parameters(TString FileName) {
double fitrange_low, fitrange_high, amp_st, mean_st, sigma_st, const_st, slope_st;
ifstream File;
File.open(FileName);
if (!File.is_open()) {
std::cout << "##### ERROR: could not open " << FileName << std::endl;
return 1;
}
std::string headerline;
getline(File, headerline);
File >> fFileName_spectrum;
getline(File, headerline);
getline(File, headerline);
while (true)
{
File >> fitrange_low >> fitrange_high >> amp_st >> mean_st >> sigma_st >> const_st >> slope_st;
ffitrange_low.push_back(fitrange_low);
ffitrange_high.push_back(fitrange_high);
famp_st.push_back(amp_st);
fmean_st.push_back(mean_st);
fsigma_st.push_back(sigma_st);
fconst_st.push_back(const_st);
fslope_st.push_back(slope_st);
if( File.eof() ) break;
}
File.close();
// print information
std::cout << "######################################" << std::endl;
std::cout << "#### Reading "<< FileName << " ..." << std::endl;
std::cout << "spectrum file name: " << fFileName_spectrum << std::endl;
std::cout << "fitrange low \t fitrange high \t counts \t mean \t sigma \t const. \t slope" << std::endl;
for (int i=0; i<ffitrange_low.size(); ++i)
{
std::cout << ffitrange_low[i] << "\t" << ffitrange_high[i] << "\t" << famp_st[i] << "\t" << fmean_st[i] << "\t" << fsigma_st[i]<< "\t" << fconst_st[i]<< "\t" << fslope_st[i]<< "\t" << std::endl;
}
std::cout << "######################################" << std::endl;
return 0;
}
// ----------------------------------------------------
// determine resolution
// ----------------------------------------------------
int main(int argc, char *argv[]) {
// argc should be 2 for correct execution
if ( argc != 2 ) {
// We print argv[0] assuming it is the program name
std::cout<<"usage: "<< argv[0] <<" <parameters_file.txt>" << std::endl;
return 1;
}
TString FileName_parameters = argv[1];
// read parameters from file
if (read_parameters(FileName_parameters)) {
return 1;
}
// get spectrum
TH1D* hist = getspectrum(fFileName_spectrum);
if (hist==0) {
return 1;
}
// draw spectrum
TCanvas* c1 = new TCanvas("c1");
hist->Draw();
TF1* fit = new TF1();
std::vector<double> peakpos, peakpos_err, sigma, sigma_err;
int npeaks = ffitrange_low.size();
// loop over all peaks
for (int i=0; i<npeaks; ++i) {
// fit peak with Gauss+Pol1
std::cout << "###################################################" << std::endl;
std::cout << "Fitting line " << i+1 << " in range " << ffitrange_low[i] << " - " << ffitrange_high[i] << std::endl;
fit = FitGaussPol1(hist,famp_st[i],fmean_st[i], fsigma_st[i], fconst_st[i], fslope_st[i], ffitrange_low[i], ffitrange_high[i]);
if (fit==0) {
return 1;
}
peakpos.push_back(fit->GetParameter(1));
peakpos_err.push_back(fit->GetParError(1));
sigma.push_back(fit->GetParameter(2));
sigma_err.push_back(fit->GetParError(2));
}
// save calibration fits
c1->SaveAs(fFileName_spectrum+"_resolution_fits.root");
c1->SaveAs(fFileName_spectrum+"_resolution_fits.pdf");
// plot peak position vs. energy
TGraphErrors* graph = new TGraphErrors(npeaks, &peakpos.at(0), &sigma.at(0), &peakpos_err.at(0), &sigma_err.at(0) );
// fit with sqrt
std::cout << "###################################################" << std::endl;
std::cout << "Fitting resolution curve ..." << std::endl;
TF1* fit_res = FitSqrt(graph,0.1,0.,0.);
if (fit_res==0) {
return 1;
}
TCanvas* c2 = new TCanvas("c2");
graph->SetTitle("Energy Resolution");
graph->SetMarkerStyle(2);
graph->GetXaxis()->SetTitle("Energy (keV)");
graph->GetYaxis()->SetTitle("Sigma (keV)");
graph->Draw("ap");
c2->SaveAs(fFileName_spectrum+"_resolution_function.root");
c2->SaveAs(fFileName_spectrum+"_resolution_function.pdf");
return 0;
}