GRSISort "v4.0.0.5"
An extension of the ROOT analysis Framework
Loading...
Searching...
No Matches
TPeakFitter.cxx
Go to the documentation of this file.
1#include "TPeakFitter.h"
2#include "Math/Minimizer.h"
3#include "TH1.h"
4#include "TClass.h"
5
6#include "Globals.h"
7
9
10TPeakFitter::TPeakFitter(const Double_t& rangeLow, const Double_t& rangeHigh)
11 : fRangeLow(rangeLow), fRangeHigh(rangeHigh)
12{
13 fBGToFit = new TF1("fbg", this, &TPeakFitter::DefaultBackgroundFunction, fRangeLow, fRangeHigh, 4, "TPeakFitter", "DefaultBackgroundFunction");
14 fBGToFit->FixParameter(3, 0);
15 fBGToFit->SetLineColor(static_cast<Color_t>(kRed + fColorIndex));
16 if(fVerboseLevel >= EVerbosity::kSubroutines) { std::cout << __PRETTY_FUNCTION__ << ": constructed peak fitter " << this << " using range " << fRangeLow << " - " << fRangeHigh << std::endl; } // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay)
17}
18
20{
21 if(fVerboseLevel >= EVerbosity::kSubroutines) { std::cout << __PRETTY_FUNCTION__ << ": destroying peak fitter " << this << " for range " << fRangeLow << " - " << fRangeHigh << std::endl; } // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay)
22}
23
24void TPeakFitter::Print(Option_t* opt) const
25{
26 /// Print information from the fit, opt is passed along to each individual peaks Print function.
27 if(fTotalFitFunction != nullptr) {
28 double xMin = 0.;
29 double xMax = 0.;
30 fTotalFitFunction->GetRange(xMin, xMax);
31 std::cout << "Fitting range: " << xMin << " to " << xMax << std::endl;
32 } else {
33 std::cout << "Range: " << fRangeLow << " to " << fRangeHigh << std::endl;
34 }
35 Int_t counter = 0;
36 if(!fPeaksToFit.empty()) {
37 std::cout << "Peaks: " << std::endl;
38 for(auto* peak : fPeaksToFit) {
39 std::cout << "Peak #" << counter++ << " - ";
40 peak->Print(opt);
41 }
42 std::cout << "Chi2/Ndf = " << fPeaksToFit.front()->GetChi2() << "/" << fPeaksToFit.front()->GetNDF() << " = " << fPeaksToFit.front()->GetReducedChi2() << std::endl;
43 } else {
44 std::cout << "Could not find peak to print" << std::endl;
45 }
46 if(fBGToFit != nullptr) {
47 std::cout << "BG: " << std::endl;
48 fBGToFit->Print(opt);
49 } else {
50 std::cout << "Could not find a BG to print" << std::endl;
51 }
52}
53
55{
56 /// Print the range of the fit and the parameters of each peak on a single line.
57 std::cout << "Range: " << fRangeLow << " to " << fRangeHigh << " - ";
58 int i = 0;
59 for(auto* peak : fPeaksToFit) {
60 std::cout << peak->IsA()->GetName() << " #" << i++ << " ";
61 peak->PrintParameters();
62 }
63 std::cout << std::endl;
64}
65
67{
68 Int_t n_par = 0;
69 for(auto* peak : fPeaksToFit) {
70 n_par += peak->GetNParameters();
71 }
72 if(fBGToFit != nullptr) {
73 n_par += fBGToFit->GetNpar();
74 }
75
76 return n_par;
77}
78
79void TPeakFitter::SetRange(const Double_t& low, const Double_t& high)
80{
81 fRangeLow = low;
82 fRangeHigh = high;
83 if(fTotalFitFunction != nullptr) {
84 fTotalFitFunction->SetTitle(Form("total_fit_%.0f_%.0f", fRangeLow, fRangeHigh));
85 }
86}
87
88TFitResultPtr TPeakFitter::Fit(TH1* fit_hist, Option_t* opt)
89{
90 /// Fit the histogram. Recognized options are "q" for a quiet fit, "retryfit" to retry a fit without parameter limits
91 /// if one of the parameters got close to its limit. All options (apart from "retryfit") are passed on to the actual
92 /// call to TH1::Fit.
93 if(fPeaksToFit.empty()) {
94 std::cout << "No peaks provided!" << std::endl;
95 return {};
96 }
97
98 TString options = opt;
99 options.ToLower();
100 bool quiet = options.Contains("q");
101 bool verbose = options.Contains("v");
102 if(quiet && verbose) {
103 std::cout << "Don't know how to be quiet and verbose at once (" << opt << "), going to be verbose!" << std::endl;
104 quiet = false;
105 }
106 bool retryFit = options.Contains("retryfit");
107 options.ReplaceAll("retryfit", "");
108
109 // store original range
110 int firstBin = fit_hist->GetXaxis()->GetFirst();
111 int lastBin = fit_hist->GetXaxis()->GetLast();
112 // set range to fit range of peak fitter
113 fit_hist->GetXaxis()->SetRangeUser(fRangeLow, fRangeHigh);
114
115 ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2", "Combination");
116 TVirtualFitter::SetMaxIterations(100000);
117 TVirtualFitter::SetPrecision(1e-4);
118 if(fTotalFitFunction == nullptr) {
119 fTotalFitFunction = new TF1(Form("total_fit_%.0f_%.0f", fRangeLow, fRangeHigh), this, &TPeakFitter::FitFunction, fRangeLow, fRangeHigh, GetNParameters(), "TPeakFitter", "FitFunction");
120 }
121 fTotalFitFunction->SetLineColor(static_cast<Color_t>(kMagenta + fColorIndex));
123 fBGToFit->SetRange(fRangeLow, fRangeHigh);
124 // We need to initialize all of the parameters based on the peak parameters
125 if(fit_hist != nullptr && (!fInitFlag || fit_hist != fLastHistFit)) {
126 if(verbose) { std::cout << "Initializing Fit...." << std::endl; }
128 InitializeParameters(fit_hist);
129 fInitFlag = true;
130 fLastHistFit = fit_hist;
131 }
133 if(verbose) { PrintParameters(); }
134 // Do a first fit to get closer on the parameters
135 fit_hist->Fit(fTotalFitFunction, "RNQ0");
136 // Now try the fit with the options provided
137 fit_hist->Fit(fTotalFitFunction, Form("RNQO%s", options.Data()));
138 // Now do a final fit with integrated bins
139 TFitResultPtr fit_res = fit_hist->Fit(fTotalFitFunction, Form("SRI%s", options.Data()));
140
141 // check the parameter errors
142 if(!TGRSIFunctions::CheckParameterErrors(fit_res, options.Data())) {
143 if(retryFit) {
144 // fit again with all parameters released
145 if(!quiet) { std::cout << GREEN << "Re-fitting with released parameters (without any limits)" << RESET_COLOR << std::endl; }
146 for(int i = 0; i < fTotalFitFunction->GetNpar(); ++i) {
147 // as of now this only works for the first peak and is reliant on the fact that all peaks have the centroig as parameter 1
148 // might be better to check if the parameter name matches "centroid_X" where X is the peak number?
149 if(i == 1) { continue; } // skipping centroid, which should always be parameter 1
150 // only release parameters with limits, not those that have been fixed
151 double min = 0.;
152 double max = 0.;
153 fTotalFitFunction->GetParLimits(i, min, max);
154 if(min != max) { fTotalFitFunction->ReleaseParameter(i); }
155 }
156 fit_res = fit_hist->Fit(fTotalFitFunction, Form("SRI%s", options.Data()));
157 TGRSIFunctions::CheckParameterErrors(fit_res, options.Data());
158 } else {
159 // re-try using minos instead of minuit (only thing from opt that might play a role it 'q' or 'v')
160 if(!quiet) { std::cout << YELLOW << "Re-fitting with \"E\" option to get better error estimation using Minos technique." << RESET_COLOR << std::endl; }
161 fit_res = fit_hist->Fit(fTotalFitFunction, Form("SRIE%s", options.Data()));
162 }
163 }
164
165 if(fit_res.Get() != nullptr) {
166 fTotalFitFunction->SetFitResult(*fit_res); // Not sure if this is needed
167
168 // Once we do the fit, we want to update all of the Peak parameters.
169 UpdatePeakParameters(fit_res, fit_hist);
170 } else if(!quiet) {
171 std::cout << RED << "Failed to get fit result, be aware that the results and especially the error estimates might be wrong!" << RESET_COLOR << std::endl;
172 }
173 fit_hist->Draw("hist");
174 fTotalFitFunction->Draw("same");
175 fPeaksToFit.front()->DrawBackground("same");
176
177 // always print the result of the fit even if not verbose
178 if(!quiet) {
179 std::cout << "****************" << std::endl
180 << "Summary of Fit: " << std::endl;
181 Print();
182 }
183 // restore old range
184 fit_hist->GetXaxis()->SetRange(firstBin, lastBin);
185 // reset the default back to minuit in case we switched to minos in between
186 ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2", "Combination");
187
188 return fit_res;
189}
190
191void TPeakFitter::UpdatePeakParameters(const TFitResultPtr& fit_res, TH1* fit_hist)
192{
193 // We enter this function once we have performed a fit and need to tell the peaks
194 // what their new parameters should be.
195 Int_t param_counter = 0;
196
197 TF1* global_bg = new TF1("global_bg", this, &TPeakFitter::BackgroundFunction, fRangeLow, fRangeHigh, fTotalFitFunction->GetNpar(), "TPeakFitter", "BackgroundFunction");
198 global_bg->SetParameters(fTotalFitFunction->GetParameters());
199
200 // Start by looping through all of the peaks in the fitter.
201 for(auto* p_it : fPeaksToFit) {
202 TF1* peak_func = p_it->GetFitFunction();
203 // We need to make a clone of the main fit function, and set all of the non-peak parameters to 0.
204 // This will allow us to do the proper integrals. We have to do the same with the covariance matrix
205 // We get the covariance matrix that we are using to update the individual peaks
206 TF1* total_function_copy = new TF1(*fTotalFitFunction);
207 TMatrixDSym covariance_matrix = fit_res->GetCovarianceMatrix();
208 bool goodCovarianceMatrix = true;
209 if(covariance_matrix.GetNrows() < peak_func->GetNpar() || covariance_matrix.GetNcols() < peak_func->GetNpar()) {
210 goodCovarianceMatrix = false;
211 }
212 // Now we need to remove all of the parts of the covariance peak that has nothing to do with this particular peak
213 // This would be the diagonals of the background, and all of the other peaks.
214 Int_t param_to_zero_counter = 0;
215 std::vector<Int_t> param_to_zero_list;
216 for(auto* other_p_it : fPeaksToFit) {
217 if(other_p_it != p_it) {
218 for(int i = 0; i < other_p_it->GetFitFunction()->GetNpar(); ++i) {
219 param_to_zero_list.push_back(param_to_zero_counter);
220 ++param_to_zero_counter;
221 }
222 } else {
223 for(int i = 0; i < peak_func->GetNpar(); ++i) {
224 if(p_it->IsBackgroundParameter(i)) {
225 param_to_zero_list.push_back(param_to_zero_counter);
226 }
227 ++param_to_zero_counter;
228 }
229 Double_t low_range = 0.;
230 Double_t high_range = 0.;
231 fTotalFitFunction->GetRange(low_range, high_range);
232 peak_func->SetRange(low_range, high_range);
233 }
234 }
235 // We have now taken care to zero all of the non-peak parameters in the peak list, so we want to remove the total background contribution
236 for(int i = param_to_zero_counter; i < fTotalFitFunction->GetNpar(); ++i) {
237 param_to_zero_list.push_back(i);
238 }
239
240 // zero other non-peak portions of matrix
241 for(auto i : param_to_zero_list) {
242 for(auto j : param_to_zero_list) {
243 if(goodCovarianceMatrix) { covariance_matrix(i, j) = 0.0; }
244 }
245 total_function_copy->SetParameter(i, 0.0);
246 }
247 if(peak_func != nullptr) {
248 for(int i = 0; i < peak_func->GetNpar(); ++i) {
249 peak_func->SetParameter(i, fTotalFitFunction->GetParameter(param_counter));
250 peak_func->SetParError(i, fTotalFitFunction->GetParError(param_counter));
251 Double_t low = 0.;
252 Double_t high = 0.;
253 fTotalFitFunction->GetParLimits(param_counter, low, high);
254 peak_func->SetParLimits(i, low, high);
255 ++param_counter;
256 // Lets do some integrals meow.
257 }
258 p_it->SetArea(total_function_copy->Integral(p_it->Centroid() - p_it->Width() * 5., p_it->Centroid() + p_it->Width() * 5., 1e-8) / fit_hist->GetBinWidth(1));
259 if(goodCovarianceMatrix) {
260 p_it->SetAreaErr(total_function_copy->IntegralError(p_it->Centroid() - p_it->Width() * 5., p_it->Centroid() + p_it->Width() * 5., total_function_copy->GetParameters(), covariance_matrix.GetMatrixArray(), 1E-5) / fit_hist->GetBinWidth(1));
261 } else if(fVerboseLevel != EVerbosity::kQuiet) {
262 std::cout << "Not setting area error because we don't have a good covariance matrix!" << std::endl;
263 }
264 }
265 total_function_copy->Delete();
266 }
267 if(fBGToFit != nullptr) {
268 for(int i = 0; i < fBGToFit->GetNpar(); ++i) {
269 fBGToFit->SetParameter(i, fTotalFitFunction->GetParameter(param_counter));
270 fBGToFit->SetParError(i, fTotalFitFunction->GetParError(param_counter));
271 ++param_counter;
272 }
273 }
274 // We now have a copy of the global background. Copy and send to every peak
275 for(auto* p_it : fPeaksToFit) {
276 p_it->SetGlobalBackground(new TF1(*global_bg));
277 p_it->SetChi2(fTotalFitFunction->GetChisquare());
278 p_it->SetNDF(fTotalFitFunction->GetNDF());
279 }
280 global_bg->Delete();
281}
282
284{
285 for(auto* p_it : fPeaksToFit) {
286 p_it->InitializeParameters(fit_hist, fRangeLow, fRangeHigh);
287 }
288}
289
291{
292 /// This functions gets the parameters and their limits from the peak functions and
293 /// sets them for the total fit function
294 Int_t param_counter = 0;
295 Int_t peak_counter = 0;
296 for(auto* p_it : fPeaksToFit) {
297 TF1* peak_func = p_it->GetFitFunction();
298 if(peak_func != nullptr) {
299 for(int i = 0; i < peak_func->GetNpar(); ++i) {
300 fTotalFitFunction->SetParName(param_counter, Form("%s_%i", peak_func->GetParName(i), peak_counter));
301 fTotalFitFunction->SetParameter(param_counter, peak_func->GetParameter(i));
302 fTotalFitFunction->SetParError(param_counter, peak_func->GetParError(i));
303 Double_t limit_low = 0.;
304 Double_t limit_high = 0.;
305 peak_func->GetParLimits(i, limit_low, limit_high);
306 fTotalFitFunction->SetParLimits(param_counter, limit_low, limit_high);
307 ++param_counter;
308 }
309 ++peak_counter;
310 }
311 }
312 if(fBGToFit != nullptr) {
313 for(int i = 0; i < fBGToFit->GetNpar(); ++i) {
314 fTotalFitFunction->SetParName(param_counter, fBGToFit->GetParName(i));
315 fTotalFitFunction->SetParameter(param_counter, fBGToFit->GetParameter(i));
316 fTotalFitFunction->SetParError(param_counter, fBGToFit->GetParError(i));
317 Double_t limit_low = 0.;
318 Double_t limit_high = 0.;
319 fBGToFit->GetParLimits(i, limit_low, limit_high);
320 fTotalFitFunction->SetParLimits(param_counter, limit_low, limit_high);
321 ++param_counter;
322 }
323 }
324}
325
326Double_t TPeakFitter::FitFunction(Double_t* dim, Double_t* par)
327{
328 // I want to use the EvalPar command here in order to get the individual peaks
329 Double_t sum = 0;
330 Int_t params_so_far = 0;
331 if(fVerboseLevel >= EVerbosity::kSubroutines) { std::cout << __PRETTY_FUNCTION__ << ": this " << this << " has " << fPeaksToFit.size() << " peaks to be evaluated at x = " << dim[0] << std::endl; } // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay)
332 for(auto* p_it : fPeaksToFit) {
333 TF1* peakFunction = p_it->GetFitFunction();
334 if(peakFunction == nullptr) {
335 std::cerr << "Failed to get fit function for peak from " << p_it << " at " << p_it->Centroid() << std::endl;
336 return 0.;
337 }
338 if(fVerboseLevel >= EVerbosity::kLoops) { std::cout << "Evaluating fit function " << peakFunction << " using " << peakFunction->GetNpar() << " parameters starting at " << params_so_far << " (" << &par[params_so_far] << ")" << std::endl; }
339 sum += peakFunction->EvalPar(dim, &par[params_so_far]);
340 params_so_far += peakFunction->GetNpar();
341 }
342 sum += fBGToFit->EvalPar(dim, &par[params_so_far]);
343
344 return sum;
345}
346
347Double_t TPeakFitter::BackgroundFunction(Double_t* dim, Double_t* par)
348{
349 Double_t sum = 0;
350 Int_t params_so_far = 0;
351 for(auto* p_it : fPeaksToFit) {
352 TF1* peak_func = p_it->GetBackgroundFunction();
353 sum += peak_func->EvalPar(dim, &par[params_so_far]);
354 params_so_far += peak_func->GetNpar();
355 }
356 sum += fBGToFit->EvalPar(dim, &par[params_so_far]);
357
358 return sum;
359}
360
362{
363 fBGToFit->SetParName(0, "A");
364 fBGToFit->SetParName(1, "B");
365 fBGToFit->SetParName(2, "C");
366 fBGToFit->SetParName(3, "bg_offset");
367
368 Double_t lowLimit = 0.;
369 Double_t highLimit = 0.;
370 fBGToFit->GetParLimits(0, lowLimit, highLimit);
371
372 double value = fBGToFit->GetParameter(0);
373 if(value == 0 && lowLimit == 0 && highLimit == 0) {
374 fBGToFit->SetParLimits(0, 0.0, fit_hist->GetBinContent(fit_hist->FindBin(fRangeHigh)) * 100.);
375 fBGToFit->SetParameter("A", fit_hist->GetBinContent(fit_hist->FindBin(fRangeHigh)));
376 }
377 fBGToFit->GetParLimits(1, lowLimit, highLimit);
378
379 value = fBGToFit->GetParameter(1);
380 if(value == 0 && lowLimit == 0 && highLimit == 0) {
381 fBGToFit->SetParameter("B", fit_hist->GetBinWidth(1) * ((fit_hist->GetBinContent(fit_hist->FindBin(fRangeLow)) - fit_hist->GetBinContent(fit_hist->FindBin(fRangeHigh))) / (fRangeLow - fRangeHigh)));
382 }
383 fBGToFit->GetParLimits(2, lowLimit, highLimit);
384
385 value = fBGToFit->GetParameter(2);
386 if(value == 0 && lowLimit == 0 && highLimit == 0) {
387 fBGToFit->SetParameter("C", 0.0000);
388 fBGToFit->FixParameter(2, 0.00);
389 }
390 fBGToFit->GetParLimits(3, lowLimit, highLimit);
391
392 value = fBGToFit->GetParameter(3);
393 if(value == 0 && lowLimit == 0 && highLimit == 0) {
394 fBGToFit->SetParLimits(3, fRangeLow, fRangeHigh);
395 if(!fPeaksToFit.empty()) {
396 fBGToFit->SetParameter("bg_offset", fPeaksToFit.front()->Centroid());
397 } else {
398 fBGToFit->SetParameter("bg_offset", (fRangeHigh + fRangeLow) / 2.);
399 }
400 }
401}
402
403Double_t TPeakFitter::DefaultBackgroundFunction(Double_t* dim, Double_t* par) // NOLINT(readability-non-const-parameter)
404{
405 Double_t x = dim[0];
406 Double_t A = par[0];
407 Double_t B = par[1];
408 Double_t C = par[2];
409 Double_t bg_offset = par[3];
410
411 return A + B * (x - bg_offset) + C * TMath::Power(x - bg_offset, 2.);
412}
413
414void TPeakFitter::DrawPeaks(Option_t*) const
415{
416 // First we are going to draw the background
417 TF1* bg_to_draw = new TF1;
418 fTotalFitFunction->Copy(*bg_to_draw);
419 bg_to_draw->SetLineColor(static_cast<Color_t>(kRed + fColorIndex));
420
421 for(auto* p_it : fPeaksToFit) {
422 TF1* peak_func = p_it->GetFitFunction();
423 TF1* total_function_copy = new TF1;
424 fTotalFitFunction->Copy(*total_function_copy);
425 // Now we need to remove all of the parts of the covariance peak that has nothing to do with this particular peak
426 // This would be the diagonals of the background, and all of the other peaks.
427 Int_t param_to_zero_counter = 0;
428 for(auto* other_p_it : fPeaksToFit) {
429 if(other_p_it != p_it) {
430 for(int i = 0; i < other_p_it->GetFitFunction()->GetNpar(); ++i) {
431 total_function_copy->SetParameter(param_to_zero_counter, 0.0);
432 ++param_to_zero_counter;
433 }
434 } else {
435 for(int i = 0; i < peak_func->GetNpar(); ++i) {
436 if(!p_it->IsBackgroundParameter(i)) {
437 bg_to_draw->SetParameter(param_to_zero_counter, 0.0);
438 }
439 ++param_to_zero_counter;
440 }
441 }
442 }
443 total_function_copy->Draw("same");
444 }
445 bg_to_draw->Draw("same");
446}
EVerbosity
Definition Globals.h:143
@ kQuiet
Definition Globals.h:144
@ kSubroutines
Definition Globals.h:146
@ kLoops
Definition Globals.h:147
#define RED
Definition Globals.h:9
#define YELLOW
Definition Globals.h:7
#define GREEN
Definition Globals.h:8
#define RESET_COLOR
Definition Globals.h:5
void UpdatePeakParameters(const TFitResultPtr &fit_res, TH1 *fit_hist)
TF1 * fTotalFitFunction
void PrintParameters() const
void SetRange(const Double_t &low, const Double_t &high)
void Print(Option_t *opt="") const override
Int_t GetNParameters() const
void InitializeBackgroundParameters(TH1 *fit_hist)
Double_t FitFunction(Double_t *dim, Double_t *par)
Double_t fRangeLow
void InitializeParameters(TH1 *fit_hist)
int fColorIndex
this index is added to the colors kRed for the total function and kMagenta for the individual peaks
void DrawPeaks(Option_t *="") const
std::list< TSinglePeak * > fPeaksToFit
Double_t BackgroundFunction(Double_t *dim, Double_t *par)
Double_t fRangeHigh
static EVerbosity fVerboseLevel
Changes verbosity of code.
TFitResultPtr Fit(TH1 *fit_hist, Option_t *opt="")
Double_t DefaultBackgroundFunction(Double_t *dim, Double_t *par)
TH1 * fLastHistFit
void UpdateFitterParameters()
bool CheckParameterErrors(const TFitResultPtr &fitres, std::string opt="")