GRSISort "v4.1.1.0"
An extension of the ROOT analysis Framework
Loading...
Searching...
No Matches
TSourceCalibration.cxx
Go to the documentation of this file.
2
3#include <chrono>
4#include <cstdlib>
5#include <iostream>
6#include <iomanip>
7#include <string>
8#include <cstring>
9#if __cplusplus >= 201703L
10#include <filesystem>
11#endif
12#include <thread>
13
14#include "TSystem.h"
15#include "TGTableLayout.h"
16#include "TCanvas.h"
17#include "TLinearFitter.h"
18#include "TF1.h"
19#include "TSpectrum.h"
20#include "TPolyMarker.h"
21#include "TObject.h"
22#include "TFrame.h"
23#include "TVirtualX.h"
24
25#include "TChannel.h"
26#include "GRootCommands.h"
27#include "combinations.h"
28#include "Globals.h"
29
30std::map<double, std::tuple<double, double, double, double>> RoughCal(std::vector<double> peaks, std::vector<std::tuple<double, double, double, double>> sources, TSourceTab* sourceTab)
31{
32 /// This function tries to match a list of found peaks (channels) to a list of provided peaks (energies).
33 /// It does so in slightly smarter way than the brute force method `Match`, by taking the reported intensity of the source peaks into account.
34
35 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << RESET_COLOR << "\"Smart\" matching " << peaks.size() << " peaks with " << sources.size() << " source energies" << std::endl; }
37 for(size_t i = 0; i < peaks.size() || i < sources.size(); ++i) {
38 std::cout << i << ".: " << std::setw(8);
39 if(i < peaks.size()) {
40 std::cout << peaks[i];
41 } else {
42 std::cout << " ";
43 }
44 std::cout << " - " << std::setw(8);
45 if(i < sources.size()) {
46 std::cout << std::get<0>(sources[i]);
47 } else {
48 std::cout << " ";
49 }
50 std::cout << std::endl;
51 }
52 }
53
54 std::map<double, std::tuple<double, double, double, double>> result;
55 std::sort(peaks.begin(), peaks.end());
56 std::sort(sources.begin(), sources.end(), [](const std::tuple<double, double, double, double>& a, const std::tuple<double, double, double, double>& b) { return std::get<2>(a) > std::get<2>(b); });
57
58 auto maxSize = std::max(peaks.size(), sources.size());
59
60 // Peaks are the fitted points.
61 // source are the known values
62
63 TLinearFitter fitter(1, "1 ++ x");
64
65 // intermediate vectors and map
66 std::vector<double> sourceValues(sources.size());
67 std::map<double, double> tmpMap;
68
69 for(size_t num_data_points = std::min(peaks.size(), sourceValues.size()); num_data_points > 0; num_data_points--) {
70 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << num_data_points << " data points:" << std::endl; }
71 double best_chi2 = DBL_MAX;
72 int iteration = 0;
73 for(auto peak_values : combinations(peaks, num_data_points)) {
74 // Add a (0,0) point to the calibration.
75 peak_values.push_back(0);
76 // instead of going through all possible combinations of the peaks with the source energies
77 // we pick the num_data_points most intense lines and try them
78 // we don't do the same with the peaks as there might be an intense background peak in the data (511 etc.)
79 sourceValues.resize(num_data_points);
80 for(size_t i = 0; i < sourceValues.size(); ++i) { sourceValues[i] = std::get<0>(sources[i]); }
81 std::sort(sourceValues.begin(), sourceValues.end());
82 sourceValues.push_back(0);
83
85 for(size_t i = 0; i < peak_values.size(); ++i) {
86 std::cout << i << ".: " << std::setw(8) << peak_values[i] << " - " << std::setw(8) << sourceValues[i] << std::endl;
87 }
88 }
89 if(peaks.size() > 3) {
90 double pratio = peak_values.front() / peak_values.at(peak_values.size() - 2);
91 double sratio = sourceValues.front() / sourceValues.at(sourceValues.size() - 2);
92 if(TSourceCalibration::VerboseLevel() > EVerbosity::kAll) { std::cout << "ratio: " << pratio << " - " << sratio << " = " << std::abs(pratio - sratio) << std::endl; }
93 if(std::abs(pratio - sratio) > 0.02) {
94 if(TSourceCalibration::VerboseLevel() > EVerbosity::kAll) { std::cout << "skipping" << std::endl; }
95 continue;
96 }
97 }
98
99 fitter.ClearPoints();
100 fitter.AssignData(sourceValues.size(), 1, peak_values.data(), sourceValues.data());
101 fitter.Eval();
102
103 if(fitter.GetChisquare() < best_chi2) {
104 tmpMap.clear();
105 for(size_t i = 0; i < num_data_points; i++) {
106 tmpMap[peak_values[i]] = sourceValues[i];
107 }
108 best_chi2 = fitter.GetChisquare();
109 }
110 sourceTab->Status(Form("%zu/%zu - %zu - %d", num_data_points, peaks.size(), maxSize, iteration), 1);
111 ++iteration;
112 if(iteration >= TSourceCalibration::MaxIterations()) { break; }
113 }
114
115 // Remove one peak value from the best fit, make sure that we reproduce (0,0) intercept.
116 if(tmpMap.size() > 2) {
117 std::vector<double> peak_values;
118 std::vector<double> source_values;
119 for(auto& item : tmpMap) {
120 peak_values.push_back(item.first);
121 source_values.push_back(item.second);
122 }
123
124 for(size_t skipped_point = 0; skipped_point < source_values.size(); skipped_point++) {
125 std::swap(peak_values[skipped_point], peak_values.back());
126 std::swap(source_values[skipped_point], source_values.back());
127
128 fitter.ClearPoints();
129 fitter.AssignData(source_values.size() - 1, 1, peak_values.data(), source_values.data());
130 fitter.Eval();
131
132 if(std::abs(fitter.GetParameter(0)) > 10) {
134 std::cout << fitter.GetParameter(0) << " too big an offset, clearing map with " << tmpMap.size() << " points: ";
135 for(auto iter : tmpMap) { std::cout << iter.first << " - " << iter.second << "; "; }
136 std::cout << std::endl;
137 }
138 tmpMap.clear();
139 break;
140 }
141
142 std::swap(peak_values[skipped_point], peak_values.back());
143 std::swap(source_values[skipped_point], source_values.back());
144 }
145 }
146
147 // copy all values from the vectors to the result map
148 if(!tmpMap.empty()) {
149 // apparently c++14 is needed to use auto in a lambda so for now we spell it out
150 // for(auto it : tmpMap) result[*(std::find_if(peaks.begin(), peaks.end(), [&it] (auto& item) { return it.first == std::get<0>(item); }))] =
151 // *(std::find_if(sources.begin(), sources.end(), [&it] (auto& item) { return it.second == std::get<0>(item); }));
152 for(auto iter : tmpMap) {
153 result[*(std::find_if(peaks.begin(), peaks.end(), [&iter](auto& item) { return iter.first == item; }))] =
154 *(std::find_if(sources.begin(), sources.end(), [&iter](auto& item) { return iter.second == std::get<0>(item); }));
155 }
157 std::cout << "Smart matched " << num_data_points << " data points from " << peaks.size() << " peaks with " << sources.size() << " source energies" << std::endl;
158 std::cout << "Returning map with " << result.size() << " points: ";
159 for(auto iter : result) { std::cout << iter.first << " - " << std::get<0>(iter.second) << "; "; }
160 std::cout << std::endl;
161 }
162 break;
163 }
164 sourceTab->Status(Form("%zu/%zu - %zu", num_data_points, peaks.size(), maxSize), 1);
165 }
166
167 return result;
168}
169
170std::map<TGauss*, std::tuple<double, double, double, double>> Match(std::vector<TGauss*> peaks, std::vector<std::tuple<double, double, double, double>> sources, TSourceTab* sourceTab)
171{
172 /// This function tries to match a list of found peaks (channels) to a list of provided peaks (energies).
173 /// It does so in a brute force fashion where we try all combinations of channels and energies, do a linear fit through them, and keep the one with the best chi square.
174
175 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << RESET_COLOR << "Matching " << peaks.size() << " peaks with " << sources.size() << " source energies" << std::endl; }
177 std::cout << "peaks:";
178 for(auto* peak : peaks) {
179 std::cout << " " << peak->Centroid();
180 }
181 std::cout << std::endl;
182 std::cout << "energies:";
183 for(auto source : sources) {
184 std::cout << " " << std::get<0>(source);
185 }
186 std::cout << std::endl;
187 }
188
189 std::map<TGauss*, std::tuple<double, double, double, double>> result;
190 std::sort(peaks.begin(), peaks.end(), [](const TGauss* a, const TGauss* b) { return a->Centroid() < b->Centroid(); });
191 std::sort(sources.begin(), sources.end());
192
193 auto maxSize = std::max(peaks.size(), sources.size());
194
195 // Peaks are the fitted points.
196 // source are the known values
197
198 TLinearFitter fitter(1, "1 ++ x");
199
200 // intermediate vectors and map
201 std::vector<double> peakValues(peaks.size());
202 for(size_t i = 0; i < peaks.size(); ++i) { peakValues[i] = peaks[i]->Centroid(); }
203 std::vector<double> sourceValues(sources.size());
204 for(size_t i = 0; i < sources.size(); ++i) { sourceValues[i] = std::get<0>(sources[i]); }
205 std::map<double, double> tmpMap;
206
207 for(size_t num_data_points = peakValues.size(); num_data_points > 0; num_data_points--) {
208 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << num_data_points << " data points:" << std::endl; }
209 double best_chi2 = DBL_MAX;
210 int iteration = 0;
211 for(auto peak_values : combinations(peakValues, num_data_points)) {
212 // Add a (0,0) point to the calibration.
213 peak_values.push_back(0);
214 for(auto source_values : combinations(sourceValues, num_data_points)) {
215 source_values.push_back(0);
216
217 if(peakValues.size() > 3) {
218 double pratio = peak_values.front() / peak_values.at(peak_values.size() - 2);
219 double sratio = source_values.front() / source_values.at(source_values.size() - 2);
220 if(TSourceCalibration::VerboseLevel() > EVerbosity::kAll) { std::cout << "ratio: " << pratio << " - " << sratio << " = " << std::abs(pratio - sratio) << " "; }
221 if(std::abs(pratio - sratio) > 0.02) {
222 if(TSourceCalibration::VerboseLevel() > EVerbosity::kAll) { std::cout << "skipping" << std::endl; }
223 continue;
224 }
225 if(TSourceCalibration::VerboseLevel() > EVerbosity::kAll) { std::cout << std::endl; }
226 }
227
228 fitter.ClearPoints();
229 fitter.AssignData(source_values.size(), 1, peak_values.data(), source_values.data());
230 fitter.Eval();
231
232 if(fitter.GetChisquare() < best_chi2) {
233 tmpMap.clear();
234 for(size_t i = 0; i < num_data_points; i++) {
235 tmpMap[peak_values[i]] = source_values[i];
236 }
238 std::cout << "Chi^2 " << fitter.GetChisquare() << " better than " << best_chi2 << ":" << std::endl;
239 for(auto iter : tmpMap) {
240 std::cout << std::setw(10) << iter.first << " - " << iter.second << std::endl;
241 }
242 }
243 best_chi2 = fitter.GetChisquare();
244 }
245 }
246 sourceTab->Status(Form("%zu/%zu - %zu - %d", num_data_points, peakValues.size(), maxSize, iteration), 1);
247 ++iteration;
248 if(iteration >= TSourceCalibration::MaxIterations()) { break; }
249 }
250
251 // Remove one peak value from the best fit, make sure that we reproduce (0,0) intercept.
252 if(tmpMap.size() > 2) {
253 std::vector<double> peak_values;
254 std::vector<double> source_values;
255 for(auto& item : tmpMap) {
256 peak_values.push_back(item.first);
257 source_values.push_back(item.second);
258 }
259
260 for(size_t skipped_point = 0; skipped_point < source_values.size(); skipped_point++) {
261 std::swap(peak_values[skipped_point], peak_values.back());
262 std::swap(source_values[skipped_point], source_values.back());
263
264 fitter.ClearPoints();
265 fitter.AssignData(source_values.size() - 1, 1, peak_values.data(), source_values.data());
266 fitter.Eval();
267
268 if(std::abs(fitter.GetParameter(0)) > 10) {
270 std::cout << fitter.GetParameter(0) << " too big, clearing map with " << tmpMap.size() << " points: ";
271 for(auto iter : tmpMap) { std::cout << iter.first << " - " << iter.second << "; "; }
272 std::cout << std::endl;
273 }
274 tmpMap.clear();
275 break;
276 }
277
278 std::swap(peak_values[skipped_point], peak_values.back());
279 std::swap(source_values[skipped_point], source_values.back());
280 }
281 }
282
283 // copy all values from the vectors to the result map
284 if(!tmpMap.empty()) {
285 for(auto iter : tmpMap) {
286 result[*(std::find_if(peaks.begin(), peaks.end(), [&iter](auto& item) { return iter.first == item->Centroid(); }))] =
287 *(std::find_if(sources.begin(), sources.end(), [&iter](auto& item) { return iter.second == std::get<0>(item); }));
288 // apparently c++14 is needed to use auto in a lambda so for now we spell it out
289 //result[*(std::find_if(peaks.begin(), peaks.end(), [&iter](std::tuple<double, double, double, double>& item) { return iter.first == std::get<0>(item); }))] =
290 // *(std::find_if(sources.begin(), sources.end(), [&iter](std::tuple<double, double, double, double>& item) { return iter.second == std::get<0>(item); }));
291 }
293 std::cout << "Matched " << num_data_points << " data points from " << peaks.size() << " peaks with " << sources.size() << " source energies" << std::endl;
294 std::cout << "Returning map with " << result.size() << " points: ";
295 for(auto iter : result) { std::cout << iter.first->Centroid() << " - " << std::get<0>(iter.second) << "; "; }
296 std::cout << std::endl;
297 }
298 break;
299 }
300 sourceTab->Status(Form("%zu/%zu - %zu", num_data_points, peakValues.size(), maxSize), 1);
301 }
302
303 return result;
304}
305
306std::map<TGauss*, std::tuple<double, double, double, double>> SmartMatch(std::vector<TGauss*> peaks, std::vector<std::tuple<double, double, double, double>> sources, TSourceTab* sourceTab)
307{
308 /// This function tries to match a list of found peaks (channels) to a list of provided peaks (energies).
309 /// It does so in slightly smarter way than the brute force method `Match`, by taking the reported intensity of the source peaks into account.
310
311 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << RESET_COLOR << "\"Smart\" matching " << peaks.size() << " peaks with " << sources.size() << " source energies" << std::endl; }
313 for(size_t i = 0; i < peaks.size() || i < sources.size(); ++i) {
314 std::cout << i << ".: " << std::setw(8);
315 if(i < peaks.size()) {
316 std::cout << peaks[i]->Centroid();
317 } else {
318 std::cout << " ";
319 }
320 std::cout << " - " << std::setw(8);
321 if(i < sources.size()) {
322 std::cout << std::get<0>(sources[i]);
323 } else {
324 std::cout << " ";
325 }
326 std::cout << std::endl;
327 }
328 }
329
330 std::map<TGauss*, std::tuple<double, double, double, double>> result;
331 std::sort(peaks.begin(), peaks.end(), [](const TGauss* a, const TGauss* b) { return a->Centroid() < b->Centroid(); });
332 std::sort(sources.begin(), sources.end(), [](const std::tuple<double, double, double, double>& a, const std::tuple<double, double, double, double>& b) { return std::get<2>(a) > std::get<2>(b); });
333
334 auto maxSize = std::max(peaks.size(), sources.size());
335
336 // Peaks are the fitted points.
337 // source are the known values
338
339 TLinearFitter fitter(1, "1 ++ x");
340
341 // intermediate vectors and map
342 std::vector<double> peakValues(peaks.size());
343 for(size_t i = 0; i < peaks.size(); ++i) { peakValues[i] = peaks[i]->Centroid(); }
344 std::vector<double> sourceValues(sources.size());
345 std::map<double, double> tmpMap;
346
347 for(size_t num_data_points = std::min(peakValues.size(), sourceValues.size()); num_data_points > 0; num_data_points--) {
348 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << num_data_points << " data points:" << std::endl; }
349 double best_chi2 = DBL_MAX;
350 int iteration = 0;
351 for(auto peak_values : combinations(peakValues, num_data_points)) {
352 // Add a (0,0) point to the calibration.
353 peak_values.push_back(0);
354 // instead of going through all possible combinations of the peaks with the source energies
355 // we pick the num_data_points most intense lines and try them
356 // we don't do the same with the peaks as there might be an intense background peak in the data (511 etc.)
357 sourceValues.resize(num_data_points);
358 for(size_t i = 0; i < sourceValues.size(); ++i) { sourceValues[i] = std::get<0>(sources[i]); }
359 std::sort(sourceValues.begin(), sourceValues.end());
360 sourceValues.push_back(0);
361
363 for(size_t i = 0; i < peak_values.size(); ++i) {
364 std::cout << i << ".: " << std::setw(8) << peak_values[i] << " - " << std::setw(8) << sourceValues[i] << std::endl;
365 }
366 }
367 if(peakValues.size() > 3) {
368 double pratio = peak_values.front() / peak_values.at(peak_values.size() - 2);
369 double sratio = sourceValues.front() / sourceValues.at(sourceValues.size() - 2);
370 if(TSourceCalibration::VerboseLevel() > EVerbosity::kAll) { std::cout << "ratio: " << pratio << " - " << sratio << " = " << std::abs(pratio - sratio) << std::endl; }
371 if(std::abs(pratio - sratio) > 0.02) {
372 if(TSourceCalibration::VerboseLevel() > EVerbosity::kAll) { std::cout << "skipping" << std::endl; }
373 continue;
374 }
375 }
376
377 fitter.ClearPoints();
378 fitter.AssignData(sourceValues.size(), 1, peak_values.data(), sourceValues.data());
379 fitter.Eval();
380
381 if(fitter.GetChisquare() < best_chi2) {
382 tmpMap.clear();
383 for(size_t i = 0; i < num_data_points; i++) {
384 tmpMap[peak_values[i]] = sourceValues[i];
385 }
386 best_chi2 = fitter.GetChisquare();
387 }
388 sourceTab->Status(Form("%zu/%zu - %zu - %d", num_data_points, peakValues.size(), maxSize, iteration), 1);
389 ++iteration;
390 if(iteration >= TSourceCalibration::MaxIterations()) { break; }
391 }
392
393 // Remove one peak value from the best fit, make sure that we reproduce (0,0) intercept.
394 if(tmpMap.size() > 2) {
395 std::vector<double> peak_values;
396 std::vector<double> source_values;
397 for(auto& item : tmpMap) {
398 peak_values.push_back(item.first);
399 source_values.push_back(item.second);
400 }
401
402 for(size_t skipped_point = 0; skipped_point < source_values.size(); skipped_point++) {
403 std::swap(peak_values[skipped_point], peak_values.back());
404 std::swap(source_values[skipped_point], source_values.back());
405
406 fitter.ClearPoints();
407 fitter.AssignData(source_values.size() - 1, 1, peak_values.data(), source_values.data());
408 fitter.Eval();
409
410 if(std::abs(fitter.GetParameter(0)) > 10) {
412 std::cout << fitter.GetParameter(0) << " too big an offset, clearing map with " << tmpMap.size() << " points: ";
413 for(auto iter : tmpMap) { std::cout << iter.first << " - " << iter.second << "; "; }
414 std::cout << std::endl;
415 }
416 tmpMap.clear();
417 break;
418 }
419
420 std::swap(peak_values[skipped_point], peak_values.back());
421 std::swap(source_values[skipped_point], source_values.back());
422 }
423 }
424
425 // copy all values from the vectors to the result map
426 if(!tmpMap.empty()) {
427 // apparently c++14 is needed to use auto in a lambda so for now we spell it out
428 // for(auto it : tmpMap) result[*(std::find_if(peaks.begin(), peaks.end(), [&it] (auto& item) { return it.first == std::get<0>(item); }))] =
429 // *(std::find_if(sources.begin(), sources.end(), [&it] (auto& item) { return it.second == std::get<0>(item); }));
430 for(auto iter : tmpMap) {
431 result[*(std::find_if(peaks.begin(), peaks.end(), [&iter](auto& item) { return iter.first == item->Centroid(); }))] =
432 *(std::find_if(sources.begin(), sources.end(), [&iter](auto& item) { return iter.second == std::get<0>(item); }));
433 }
435 std::cout << "Smart matched " << num_data_points << " data points from " << peaks.size() << " peaks with " << sources.size() << " source energies" << std::endl;
436 std::cout << "Returning map with " << result.size() << " points: ";
437 for(auto iter : result) { std::cout << iter.first->Centroid() << " - " << std::get<0>(iter.second) << "; "; }
438 std::cout << std::endl;
439 }
440 break;
441 }
442 sourceTab->Status(Form("%zu/%zu - %zu", num_data_points, peakValues.size(), maxSize), 1);
443 }
444
445 return result;
446}
447
448double Polynomial(double* x, double* par) // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay, readability-non-const-parameter)
449{
450 double result = par[1];
451 for(int i = 1; i <= par[0]; ++i) {
452 result += par[i + 1] * TMath::Power(x[0], i);
453 }
454 return result;
455}
456
457bool FilledBin(TH2* matrix, const int& bin)
458{
459 return (matrix->Integral(bin, bin, 1, matrix->GetNbinsY()) > 1000);
460}
461
462//////////////////////////////////////// TSourceTab ////////////////////////////////////////
463TSourceTab::TSourceTab(TSourceCalibration* sourceCal, TChannelTab* channel, TGCompositeFrame* frame, GH1D* projection, const char* sourceName, std::vector<std::tuple<double, double, double, double>> sourceEnergy)
464 : fSourceCalibration(sourceCal), fChannelTab(channel), fSourceFrame(frame), fProjection(projection), fSourceName(sourceName), fSourceEnergy(std::move(sourceEnergy))
465{
466 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << DCYAN << __PRETTY_FUNCTION__ << RESET_COLOR << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
467 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << DCYAN << this << ": const. 1 fProjection = " << fProjection << " called " << fProjection->GetName() << RESET_COLOR << std::endl; }
468 // default sorting of tuples is by the first entry, which is the source energy here, so no need for a custom sort condition
469 std::sort(fSourceEnergy.begin(), fSourceEnergy.end());
471 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << DCYAN << "interface built, finding peaks ... " << RESET_COLOR << std::endl; }
472 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << DCYAN << this << ": const. 2 fProjection = " << fProjection << " called " << fProjection->GetName() << RESET_COLOR << std::endl; }
473}
474
476{
477 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << DCYAN << __PRETTY_FUNCTION__ << RESET_COLOR << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
482
483 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << DCYAN << this << ": copy const. 1 fProjection = " << fProjection << " called " << fProjection->GetName() << RESET_COLOR << std::endl; }
485 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << DCYAN << this << ": copy const. 2 fProjection = " << fProjection << " called " << fProjection->GetName() << RESET_COLOR << std::endl; }
486 fData = rhs.fData;
487 fFwhm = rhs.fFwhm;
488 // ? not sure what to do here, clearing the vectors seems unnecessary since we never filled them
489 // but we also can't just simply copy them, we would need a deep copy, which might not make sense?
490 // does is make sense in general to have assignment constructors for this class?
491 fFits.clear();
492 fBadFits.clear();
493 fPeaks.clear();
494}
495
497{
498 for(auto* fit : fFits) {
499 delete fit;
500 }
501 fFits.clear();
502 for(auto* fit : fBadFits) {
503 delete fit;
504 }
505 fBadFits.clear();
506 for(auto* peak : fPeaks) {
507 delete peak;
508 }
509 fPeaks.clear();
510 delete fSourceFrame;
511 delete fProjectionCanvas;
512 delete fSourceStatusBar;
513 //delete fProjection;
514 //delete fData;
515 //delete fFwhm;
516}
517
519{
520 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << DCYAN << __PRETTY_FUNCTION__ << RESET_COLOR << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
521 //if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << "Creating unique lock on graphics mutex to build interface " << std::endl; }
522 //std::unique_lock<std::mutex> graphicsLock(fSourceCalibration->GraphicsMutex());
523 //if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << "Got unique lock on graphics mutex!" << std::endl; }
524 // frame with canvas and status bar
525
526 fProjectionCanvas = new TRootEmbeddedCanvas(Form("ProjectionCanvas%s", fProjection->GetName()), fSourceFrame, TSourceCalibration::PanelWidth(), TSourceCalibration::PanelHeight());
527
528 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << "Created projection canvas " << fProjectionCanvas->GetName() << "!" << std::endl; }
529
530 fSourceFrame->AddFrame(fProjectionCanvas, new TGLayoutHints(kLHintsTop | kLHintsExpandX | kLHintsExpandY, 2, 2, 2, 2));
531
532 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << "Added projection canvas to source frame!" << std::endl; }
533
535 std::array<int, 3> parts = {35, 50, 15};
536 fSourceStatusBar->SetParts(parts.data(), parts.size());
537
538 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << "Created status bar!" << std::endl; }
539
540 fSourceFrame->AddFrame(fSourceStatusBar, new TGLayoutHints(kLHintsBottom | kLHintsExpandX, 2, 2, 2, 2));
541
542 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << DCYAN << this << ": buildinterface fProjection = " << fProjection << " called " << fProjection->GetName() << RESET_COLOR << std::endl; }
543 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << "Releasing unique lock on graphics mutex after building interface" << std::endl; }
544}
545
547{
548 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << DCYAN << __PRETTY_FUNCTION__ << RESET_COLOR << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
549 //std::unique_lock<std::mutex> graphicsLock(fSourceCalibration->GraphicsMutex());
550 fProjectionCanvas->GetCanvas()->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", "TSourceTab", this, "ProjectionStatus(Int_t,Int_t,Int_t,TObject*)");
551 fProjectionCanvas->Connect("ProcessedEvent(Event_t*)", "TSourceTab", this, "ProjectionStatus(Event_t*)");
552}
553
555{
556 //std::unique_lock<std::mutex> graphicsLock(fSourceCalibration->GraphicsMutex());
557 fProjectionCanvas->Disconnect("ProcessedEvent(Event_t*)", this, "ProjectionStatus(Event_t*)");
558 fProjectionCanvas->GetCanvas()->Disconnect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", this, "ProjectionStatus(Int_t,Int_t,Int_t,TObject*)");
559}
560
562{
563 // enum EGEventType {
564 // kGKeyPress = 0, kKeyRelease, = 1, kButtonPress = 2, kButtonRelease = 3,
565 // kMotionNotify = 4, kEnterNotify = 5, kLeaveNotify = 6, kFocusIn = 7, kFocusOut = 8,
566 // kExpose = 9, kConfigureNotify = 10, kMapNotify = 11, kUnmapNotify = 12, kDestroyNotify = 13,
567 // kClientMessage = 14, kSelectionClear = 15, kSelectionRequest = 16, kSelectionNotify = 17,
568 // kColormapNotify = 18, kButtonDoubleClick = 19, kOtherEvent = 20
569 // };
570 //if(TSourceCalibration::VerboseLevel() == EVerbosity::kAll) {
571 // std::cout << __PRETTY_FUNCTION__ << ": code " << event->fCode << ", count " << event->fCount << ", state " << event->fState << ", type " << event->fType << std::endl; // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay)
572 // if(event->fType == kClientMessage) {
573 // std::cout << "Client message: " << event->fUser[0] << ", " << event->fUser[1] << ", " << event->fUser[2] << std::endl;
574 // }
575 //}
576 if(event->fType == kEnterNotify) {
578 std::cout << "Entered source tab => changing gPad from " << gPad->GetName();
579 }
580 gPad = fProjectionCanvas->GetCanvas();
582 std::cout << " to " << gPad->GetName() << std::endl;
583 }
584 }
585}
586
587void TSourceTab::ProjectionStatus(Int_t event, Int_t px, Int_t py, TObject* selected)
588{
589 // looks like 51 is hovering over the object, 52 is moving the cursor over the object, and 53 is moving the cursort away from the object
590 // kButton1 = left mouse button, kButton2 = right mouse button
591 // enum EEventType {
592 // kNoEvent = 0,
593 // kButton1Down = 1, kButton2Down = 2, kButton3Down = 3, kKeyDown = 4,
594 // kWheelUp = 5, kWheelDown = 6, kButton1Shift = 7, kButton1ShiftMotion = 8,
595 // kButton1Up = 11, kButton2Up = 12, kButton3Up = 13, kKeyUp = 14,
596 // kButton1Motion = 21, kButton2Motion = 22, kButton3Motion = 23, kKeyPress = 24,
597 // kArrowKeyPress = 25, kArrowKeyRelease = 26,
598 // kButton1Locate = 41, kButton2Locate = 42, kButton3Locate = 43, kESC = 27,
599 // kMouseMotion = 51, kMouseEnter = 52, kMouseLeave = 53,
600 // kButton1Double = 61, kButton2Double = 62, kButton3Double = 63
601 //};
602 if(selected == nullptr) { return; }
603 Status(selected->GetName(), 0);
604 Status(selected->GetObjectInfo(px, py), 1);
605 //auto key = static_cast<char>(px);
606 switch(event) {
607 case kButton1Down:
608 if(selected == fProjection) {
609 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << DCYAN << "Adding new marker at " << px << ", " << py << " (pixel to user: " << fProjectionCanvas->GetCanvas()->PixeltoX(px) << ", " << fProjectionCanvas->GetCanvas()->PixeltoY(py - fProjectionCanvas->GetCanvas()->GetWh()) << ", abs pixel to user: " << fProjectionCanvas->GetCanvas()->AbsPixeltoX(px) << ", " << fProjectionCanvas->GetCanvas()->AbsPixeltoY(py) << ")" << std::endl; }
611 std::cout << fPeaks.size() << " peaks:";
612 for(auto* peak : fPeaks) {
613 std::cout << " " << peak->Centroid();
614 }
615 std::cout << std::endl;
616 }
617 //auto* polym = static_cast<TPolyMarker*>(fProjection->GetListOfFunctions()->FindObject("TPolyMarker"));
618 //if(polym == nullptr) {
619 // std::cerr << "No peaks defined yet?" << std::endl;
620 // return;
621 //}
622 //polym->SetNextPoint(fProjectionCanvas->GetCanvas()->AbsPixeltoX(px), fProjectionCanvas->GetCanvas()->AbsPixeltoX(py));
623 double range = fSourceCalibration->FitRange(); // * fProjection->GetXaxis()->GetBinWidth(1);
624 auto* pf = new TPeakFitter(fProjectionCanvas->GetCanvas()->AbsPixeltoX(px) - range, fProjectionCanvas->GetCanvas()->AbsPixeltoX(px) + range);
625 auto* peak = new TGauss(fProjectionCanvas->GetCanvas()->AbsPixeltoX(px));
626 pf->AddPeak(peak);
627 pf->Fit(fProjection, "qnretryfit");
628 if(peak->Area() > 0) {
629 fPeaks.push_back(peak);
630 fFits.push_back(pf);
631 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << "Fitted peak " << fProjectionCanvas->GetCanvas()->AbsPixeltoX(px) - range << " - " << fProjectionCanvas->GetCanvas()->AbsPixeltoX(px) + range << " -> centroid " << peak->Centroid() << std::endl; }
632 } else {
633 std::cout << "Ignoring peak at " << peak->Centroid() << " with negative area " << peak->Area() << std::endl;
634 }
635 fProjection->GetListOfFunctions()->Remove(peak);
636 fProjection->Sumw2(false); // turn errors off, makes the histogram look like a histogram when using normal plotting (hist and samefunc doesn't work for some reason)
637
638 std::sort(fPeaks.begin(), fPeaks.end(), [](const TGauss* a, const TGauss* b) { return a->Centroid() < b->Centroid(); });
639
640 auto map = Match(fPeaks, fSourceEnergy, this);
641
643
644 UpdateFits();
645
646 Add(map);
647
648 // update status
649 Status(Form("%d/%d", static_cast<int>(fData->GetN()), static_cast<int>(fPeaks.size())), 2);
650
651 // update data and re-calibrate
655 } else if(selected->IsA() == GH1D::Class() && fProjection != nullptr) {
656 std::cout << "selected object " << selected << " is a histogram called " << selected->GetName() << " and not the projection " << fProjection << " called " << fProjection->GetName() << ", going to redraw the projection" << std::endl;
657 Draw();
659 std::cout << DCYAN << "Selected object " << selected << " (" << selected->ClassName() << " called " << selected->GetName() << ") does not match projection " << fProjection << " (" << fProjection->ClassName() << " called " << fProjection->GetName() << ")" << std::endl;
660 }
661 break;
662 case kArrowKeyPress:
663 //Move1DHistogram(px, fProjection);
664 //fProjectionCanvas->GetCanvas()->Modified();
665 case kArrowKeyRelease:
666 break;
667 case kKeyDown:
668 case kKeyPress:
669 //switch(key) {
670 // case 'u':
671 // fProjection->GetXaxis()->UnZoom();
672 // fProjection->GetYaxis()->UnZoom();
673 // break;
674 // case 'U':
675 // fProjection->GetYaxis()->UnZoom();
676 // break;
677 // default:
678 // std::cout << "Key press '" << key << "' not recognized!" << std::endl;
679 // break;
680 //}
681 case kKeyUp:
682 break;
683 case kButton1Motion:
684 break;
685 case kButton1Up:
686 break;
687 case kMouseMotion:
688 //{
689 // auto* canvas = fProjectionCanvas->GetCanvas();
690 // if(canvas == nullptr) {
691 // Status("canvas nullptr", 0);
692 // break;
693 // }
694 // auto* pad = canvas->GetSelectedPad();
695 // if(pad == nullptr) {
696 // Status("pad nullptr", 0);
697 // break;
698 // }
699
700 // Status(Form("x = %f, y = %f", pad->AbsPixeltoX(px), pad->AbsPixeltoY(py)), 0);
701 //}
702 break;
703 case kMouseEnter:
704 break;
705 case kMouseLeave:
706 break;
707 default:
709 std::cout << "unprocessed event " << event << " with px, py " << px << ", " << py << " selected object is a " << selected->ClassName() << " with name \"" << selected->GetName() << "\" and title \"" << selected->GetTitle() << "\" object info is " << selected->GetObjectInfo(px, py) << std::endl;
710 }
711 break;
712 }
713}
714
716{
717 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << "UpdateRegions: clearing " << fRegions.size() << " regions" << std::endl; }
718 fRegions.clear();
720 std::cout << "projection has " << fProjection->ListOfRegions()->GetEntries() << " regions:" << std::endl;
721 fProjection->ListOfRegions()->Print();
722 }
724 for(auto* obj : *(fProjection->ListOfRegions())) {
725 if(obj->InheritsFrom(TRegion::Class())) {
726 auto* region = static_cast<TRegion*>(obj);
727 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << "UpdateRegions: found " << fRegions.size() << ". region: " << std::min(region->GetX1(), region->GetX2()) << " - " << std::max(region->GetX1(), region->GetX2()) << std::endl; }
728 if(TSourceCalibration::VerboseLevel() > EVerbosity::kLoops) { region->Print(); }
729 fRegions.emplace_back(std::min(region->GetX1(), region->GetX2()), std::max(region->GetX1(), region->GetX2()));
731 std::cout << obj->GetName() << " is not a TRegion but a " << obj->ClassName() << std::endl;
732 }
733 }
734}
735
737{
738 return TSourceCalibration::AcceptBadFits() || (peak->CentroidErr() < 0.1 * peak->Centroid() && peak->AreaErr() < peak->Area() && !std::isnan(peak->CentroidErr()) && !std::isnan(peak->AreaErr()));
739}
740
742{
743 /// This functions removes all fits named "gauss_total" from the histogram and adds instead the fit functions from all good and bad peaks to it.
744 TList* functions = fProjection->GetListOfFunctions();
746 std::cout << "Updating functions for histogram " << fProjection->GetName() << " by deleting " << functions->GetEntries() << " functions:" << std::endl;
747 for(auto* obj : *functions) {
748 std::cout << obj << ": " << obj->GetName() << "/" << obj->GetTitle() << " - " << obj->ClassName();
749 if(obj->IsA() == TF1::Class()) {
750 std::cout << " line color " << static_cast<TF1*>(obj)->GetLineColor();
751 } else {
752 std::cout << " not a TF1";
753 }
754 std::cout << std::endl;
755 }
756 }
757 // remove all old fits (all have the same name/title of "gauss_total")
758 // alternative: clone the peaks with indices in their names, that would also allow us to add a function to replace a single fit
759 TObject* fit = nullptr;
760 while((fit = functions->FindObject("gauss_total")) != nullptr) {
761 functions->Remove(fit);
762 }
763 // add all the fit functions of the good and the bad peaks
764 for(auto* obj : fFits) {
765 functions->Add(obj->GetFitFunction()->Clone("gauss_total"));
766 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << "Added clone of " << obj->GetFitFunction() << " = " << functions->Last() << std::endl; }
767 }
768 for(auto* obj : fBadFits) {
769 functions->Add(obj->GetFitFunction()->Clone("gauss_total"));
770 }
772 std::cout << "And adding " << functions->GetEntries() << " functions instead:" << std::endl;
773 for(auto* obj : *functions) {
774 std::cout << obj << ": " << obj->GetName() << "/" << obj->GetTitle() << " - " << obj->ClassName();
775 if(obj->IsA() == TF1::Class()) {
776 std::cout << " line color " << static_cast<TF1*>(obj)->GetLineColor();
777 } else {
778 std::cout << " not a TF1";
779 }
780 std::cout << std::endl;
781 }
782 }
783 if(!fFits.empty() && fFits.front()->NumberOfPeaks() > 0 && fFits.back()->NumberOfPeaks() > 0) {
784 fProjection->SetTitle(Form("%lu of %lu source peaks used, channel %.0f - %.0f", fFits.size(), fSourceEnergy.size(), fFits.front()->Peak(0)->Centroid(), fFits.back()->Peak(0)->Centroid()));
785 } else {
786 fProjection->SetTitle(Form("%lu of %lu source peaks used, channel ???? - ????", fFits.size(), fSourceEnergy.size()));
787 }
788 fProjectionCanvas->GetCanvas()->Modified();
789 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << DCYAN << this << ": updatefits fProjection = " << fProjection << " called " << fProjection->GetName() << RESET_COLOR << std::endl; }
790}
791
793{
794 fProjectionCanvas->GetCanvas()->cd();
795 fProjection->Draw(); // seems like hist + samefunc does not work. Could use hist + loop over list of functions (with same)
796 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << __PRETTY_FUNCTION__ << ": drew " << fProjection << " = " << fProjection->GetName() << " on " << fProjectionCanvas->GetCanvas()->GetName() << "/" << gPad->GetName() << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
797}
798
799void TSourceTab::InitialCalibration(const bool& force)
800{
801 /// This functions finds the peaks in the histogram, fits them, and adds the fits to the list of peaks.
802 /// This list is then used to find all peaks that lie on a straight line.
803
804 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << DCYAN << __PRETTY_FUNCTION__ << std::flush << " " << fProjection->GetName() << ": got " << fPeaks.size() << " peaks" << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
805
806 if(fPeaks.empty() || fData == nullptr || force) {
808 std::cout << __PRETTY_FUNCTION__ << ": # peaks " << fPeaks.size() << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
809 }
810 fPeaks.clear();
811 fFits.clear();
812 fBadFits.clear();
813 // Remove all associated functions from projection.
814 // These are the poly markers, fits, and the pave stat
816 TList* functions = fProjection->GetListOfFunctions();
817 std::cout << functions->GetEntries() << " functions in projection " << fProjection->GetName() << " before clearing" << std::endl;
818 for(auto* obj : *functions) {
819 std::cout << obj->GetName() << "/" << obj->GetTitle() << " - " << obj->ClassName() << std::endl;
820 }
821 }
822 fProjection->GetListOfFunctions()->Clear();
823 TSpectrum spectrum;
824 int nofPeaks = 0;
825 std::vector<double> peakPos;
827 if(fRegions.empty()) {
828 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << "No regions found, using whole spectrum" << std::endl; }
830 nofPeaks = spectrum.GetNPeaks();
831 peakPos.insert(peakPos.end(), spectrum.GetPositionX(), spectrum.GetPositionX() + nofPeaks);
832 } else {
833 for(auto& region : fRegions) {
834 fProjection->GetXaxis()->SetRangeUser(region.first, region.second);
836 int tmpNofPeaks = spectrum.GetNPeaks();
837 peakPos.insert(peakPos.end(), spectrum.GetPositionX(), spectrum.GetPositionX() + tmpNofPeaks);
838 nofPeaks += spectrum.GetNPeaks();
839 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << __PRETTY_FUNCTION__ << ": found " << spectrum.GetNPeaks() << " peaks in region " << region.first << " - " << region.second << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
840 }
841 fProjection->GetXaxis()->UnZoom();
842 }
843 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << __PRETTY_FUNCTION__ << ": found " << nofPeaks << " peaks" << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
844
845 auto map = RoughCal(peakPos, fSourceEnergy, this);
846 Add(map);
847
848 // update status
849 Status(Form("%d/%d", static_cast<int>(fData->GetN()), nofPeaks), 2);
850 }
851}
852
853void TSourceTab::FindCalibratedPeaks(const TF1* calibration)
854{
855 /// This functions finds uses the existing calibration to find all peaks from the list of source energies.
856
857 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << DCYAN << __PRETTY_FUNCTION__ << std::flush << " " << fProjection->GetName() << ": using calibration " << calibration << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
858
859 if(calibration == nullptr) {
860 std::cerr << __PRETTY_FUNCTION__ << ": provided calibration is nullptr, this function can only be called after an initial calibration has been made!" << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
861 return;
862 }
863
864 // try and make sure the right canvas is selected as all fits are drawn in the currently selected canvas
865 gPad = fProjectionCanvas->GetCanvas();
866
867 fPeaks.clear();
868 fFits.clear();
869 fBadFits.clear();
870 // Remove all associated functions from projection.
871 // These are the poly markers, fits, and the pave stat
873 std::cout << "calibration has " << calibration->GetNpar() << " parameters (";
874 for(int i = 0; i < calibration->GetNpar(); ++i) {
875 if(i != 0) {
876 std::cout << ", ";
877 }
878 std::cout << i << ": " << calibration->GetParameter(i);
879 }
880 std::cout << "), if linear offset is " << calibration->Eval(0) << " and gain is " << calibration->Eval(1) - calibration->Eval(0) << std::endl;
881 TList* functions = fProjection->GetListOfFunctions();
882 std::cout << functions->GetEntries() << " functions in projection " << fProjection->GetName() << " before clearing" << std::endl;
883 for(auto* obj : *functions) {
884 std::cout << obj->GetName() << "/" << obj->GetTitle() << " - " << obj->ClassName() << std::endl;
885 }
886 }
887 fProjection->GetListOfFunctions()->Clear();
888
889 std::map<TGauss*, std::tuple<double, double, double, double>> map;
890
891 // loop over all source energies and check if we want to use the peak or if it's too weak
892 // for each source energy we calculate the channel it would be at
893 std::vector<std::pair<std::tuple<double, double, double, double>, double>> channelPos;
894 for(auto& tuple : fSourceEnergy) {
895 if(std::get<2>(tuple) < TSourceCalibration::MinIntensity()) {
897 std::cout << "Skipping peak at " << std::get<0>(tuple) << " keV with intensity " << std::get<2>(tuple) << " below required minimum intensity " << TSourceCalibration::MinIntensity() << std::endl;
898 }
899 continue;
900 }
901 channelPos.emplace_back(tuple, calibration->GetX(std::get<0>(tuple)));
903 std::cout << "Set channel position for " << channelPos.size() - 1 << ". peak to " << channelPos.back().second << " channels = " << std::get<0>(channelPos.back().first) << " keV" << std::endl;
904 }
905 }
906 // loop over channel positions and fit the peaks
907 for(size_t i = 0; i < channelPos.size(); ++i) {
908 double lowRange = fSourceCalibration->FitRange();
909 double highRange = fSourceCalibration->FitRange();
910 // if the range for this peaks overlaps with the previous/next one, reduce the range but not more than half
911 // also check if the previous/next one is above the minimum intensity, because otherwise we would have skipped that one anyway
912 if(i > 0) {
913 if(channelPos[i].second - lowRange < channelPos[i - 1].second + highRange) {
915 std::cout << channelPos[i].second - lowRange << " < " << channelPos[i - 1].second + highRange << ": changing low range from " << lowRange << " to " << std::max(lowRange / 2., (channelPos[i].second - channelPos[i - 1].second) / 2.) << " = std::max(" << lowRange / 2. << ", " << (channelPos[i].second - channelPos[i - 1].second) / 2. << " = (" << channelPos[i].second << " - " << channelPos[i - 1].second << "/2.)" << std::endl;
916 }
917 lowRange = std::max(lowRange / 2., (channelPos[i].second - channelPos[i - 1].second) / 2.);
918 }
919 }
920 if(i + 1 < channelPos.size()) {
921 if(channelPos[i].second + highRange > channelPos[i + 1].second - lowRange) {
923 std::cout << channelPos[i].second + highRange << " > " << channelPos[i + 1].second - lowRange << ": changing high range from " << highRange << " to " << std::max(highRange / 2., (channelPos[i + 1].second - channelPos[i].second) / 2.) << " = std::max(" << highRange / 2. << ", " << (channelPos[i + 1].second - channelPos[i].second) / 2. << " = (" << channelPos[i + 1].second << " - " << channelPos[i].second << "/2.)" << std::endl;
924 }
925 highRange = std::min(highRange / 2., (channelPos[i + 1].second - channelPos[i].second) / 2.);
926 }
927 }
929 std::cout << "Trying to fit peak " << channelPos[i].second << " in range " << channelPos[i].second - lowRange << " - " << channelPos[i].second + highRange << std::endl;
930 }
931 // if we aren't already redirecting the output, we redirect it to /dev/null and do a quiet fit, otherwise we do a normal fit since the output will be redirected to a file
932 auto* pf = new TPeakFitter(channelPos[i].second - lowRange, channelPos[i].second + highRange);
933 auto* peak = new TGauss(channelPos[i].second);
934 pf->AddPeak(peak);
935 if(TSourceCalibration::LogFile().empty()) {
936 TRedirect redirect("/dev/null");
937 pf->Fit(fProjection, "qnretryfit");
938 } else {
939 pf->Fit(fProjection, "nretryfit");
940 }
941 if(peak->Area() > 0) {
942 // check if we accept bad fits (and if it was a bad fit)
943 if(Good(peak, channelPos[i].second - lowRange, channelPos[i].second + highRange)) {
944 fPeaks.push_back(peak);
945 fFits.push_back(pf);
946 map[peak] = channelPos[i].first;
948 std::cout << "Fitted peak " << channelPos[i].second - lowRange << " - " << channelPos[i].second + highRange << " -> centroid " << peak->Centroid() << ", area " << peak->Area() << std::endl;
949 }
950 } else {
951 fBadFits.push_back(pf);
953 std::cout << "We're ignoring bad fits, and this one has position " << peak->Centroid() << " +- " << peak->CentroidErr() << ", area " << peak->Area() << " +- " << peak->AreaErr() << std::endl;
954 }
955 }
957 std::cout << "Ignoring peak at " << peak->Centroid() << " with negative area " << peak->Area() << std::endl;
958 }
959 }
960
961 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << __PRETTY_FUNCTION__ << ": added " << fPeaks.size() << " peaks (" << fSourceEnergy.size() << " source energies, " << channelPos.size() << " channel positions)" << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
962 fProjection->Sumw2(false); // turn errors off, makes the histogram look like a histogram when using normal plotting (hist and samefunc doesn't work for some reason)
963
964 std::sort(fFits.begin(), fFits.end(), [](const TPeakFitter* a, const TPeakFitter* b) { return a->GetFitFunction()->GetParameter("centroid_0") < b->GetFitFunction()->GetParameter("centroid_0"); });
965 std::sort(fBadFits.begin(), fBadFits.end(), [](const TPeakFitter* a, const TPeakFitter* b) { return a->GetFitFunction()->GetParameter("centroid_0") < b->GetFitFunction()->GetParameter("centroid_0"); });
966 std::sort(fPeaks.begin(), fPeaks.end(), [](const TGauss* a, const TGauss* b) { return a->Centroid() < b->Centroid(); });
967
969
970 UpdateFits();
971
972 Add(map);
973
975 std::cout << __PRETTY_FUNCTION__ << ": found " << fData->GetN() << " peaks"; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
976 Double_t* x = fData->GetX();
977 Double_t* y = fData->GetY();
978 for(int i = 0; i < fData->GetN(); ++i) {
979 std::cout << "; " << x[i] << " - " << y[i];
980 }
981 std::cout << std::endl;
982 }
983}
984
985void TSourceTab::Add(std::map<double, std::tuple<double, double, double, double>> map)
986{
988 std::cout << DCYAN << "Adding map (double->tuple) with " << map.size() << " data points, fData = " << fData << std::endl;
989 }
990 if(fData == nullptr) {
991 fData = new TGraphErrors(map.size());
992 } else {
993 fData->Set(map.size());
994 }
995 fData->SetLineColor(2);
996 fData->SetMarkerColor(2);
997 int i = 0;
998 for(auto iter = map.begin(); iter != map.end();) {
999 // more readable variable names
1000 auto peakPos = iter->first;
1001 auto energy = std::get<0>(iter->second);
1002 fData->SetPoint(i, peakPos, energy);
1004 std::cout << "Using peak with position " << peakPos << ", energy " << energy << std::endl;
1005 }
1006 ++iter;
1007 ++i;
1008 }
1010 std::cout << "Accepted " << map.size() << " peaks" << std::endl;
1011 }
1012 fData->Sort();
1013}
1014
1015void TSourceTab::Add(std::map<TGauss*, std::tuple<double, double, double, double>> map)
1016{
1018 std::cout << DCYAN << "Adding map (TGauss*->tuple) with " << map.size() << " data points, fData = " << fData << std::endl;
1019 }
1020 if(fData == nullptr) {
1021 fData = new TGraphErrors(map.size());
1022 } else {
1023 fData->Set(map.size());
1024 }
1025 fData->SetLineColor(2);
1026 fData->SetMarkerColor(2);
1027 if(fFwhm == nullptr) {
1028 fFwhm = new TGraphErrors(map.size());
1029 } else {
1030 fFwhm->Set(map.size());
1031 }
1032 fFwhm->SetLineColor(4);
1033 fFwhm->SetMarkerColor(4);
1034 int i = 0;
1035 for(auto iter = map.begin(); iter != map.end();) {
1036 // more readable variable names
1037 auto peakPos = iter->first->Centroid();
1038 auto peakPosErr = iter->first->CentroidErr();
1039 auto peakArea = iter->first->Area();
1040 auto peakAreaErr = iter->first->AreaErr();
1041 auto fwhm = iter->first->FWHM();
1042 auto fwhmErr = iter->first->FWHMErr();
1043 auto energy = std::get<0>(iter->second);
1044 auto energyErr = std::get<1>(iter->second);
1045 // drop this peak if the uncertainties in area or position are too large
1046 if(!Good(iter->first)) {
1048 std::cout << "Dropping peak with position " << peakPos << " +- " << peakPosErr << ", area " << peakArea << " +- " << peakAreaErr << ", energy " << energy << std::endl;
1049 }
1050 iter = map.erase(iter);
1052 std::cout << "Erasing peak returned iterator " << std::distance(map.begin(), iter) << std::endl;
1053 }
1054 } else {
1055 fData->SetPoint(i, peakPos, energy);
1056 fData->SetPointError(i, peakPosErr, energyErr);
1057 fFwhm->SetPoint(i, peakPos, fwhm);
1058 fFwhm->SetPointError(i, peakPosErr, fwhmErr);
1060 std::cout << "Using peak with position " << peakPos << " +- " << peakPosErr << ", area " << peakArea << " +- " << peakAreaErr << ", energy " << energy << std::endl;
1061 }
1062 ++iter;
1063 ++i;
1064 }
1065 }
1067 std::cout << "Accepted " << map.size() << " peaks" << std::endl;
1068 }
1069 // if we dropped a peak, we need to resize the graph, if we haven't dropped any this doesn't do anything
1070 fData->Set(map.size());
1071 fFwhm->Set(map.size());
1072 fData->Sort();
1073 fFwhm->Sort();
1074 // split poly marker into those peaks that were used, and those that weren't
1075 TList* functions = fProjection->GetListOfFunctions();
1077 std::cout << functions->GetEntries() << " functions in projection " << fProjection->GetName() << std::endl;
1078 for(auto* obj : *functions) {
1079 std::cout << obj->GetName() << "/" << obj->GetTitle() << " - " << obj->ClassName() << std::endl;
1080 }
1081 }
1082 auto* polym = static_cast<TPolyMarker*>(functions->FindObject("TPolyMarker"));
1083 if(polym != nullptr) {
1084 double* oldX = polym->GetX();
1085 double* oldY = polym->GetY();
1086 int size = polym->GetN();
1087 std::vector<double> newX;
1088 std::vector<double> newY;
1089 std::vector<double> unusedX;
1090 std::vector<double> unusedY;
1091 for(i = 0; i < size; ++i) {
1092 bool used = false;
1093 for(auto point : map) {
1094 if(TMath::Abs(oldX[i] - point.first->Centroid()) < fSourceCalibration->Sigma()) {
1095 newX.push_back(oldX[i]);
1096 newY.push_back(oldY[i]);
1097 used = true;
1098 break;
1099 }
1100 }
1101 if(!used) {
1102 unusedX.push_back(oldX[i]);
1103 unusedY.push_back(oldY[i]);
1104 }
1105 }
1106 polym->SetPolyMarker(newX.size(), newX.data(), newY.data());
1107 auto* unusedMarkers = new TPolyMarker(unusedX.size(), unusedX.data(), unusedY.data());
1108 unusedMarkers->SetMarkerStyle(23); // triangle down
1109 unusedMarkers->SetMarkerColor(16); // light grey
1110 functions->Add(unusedMarkers);
1112 std::cout << fProjection->GetTitle() << ": " << size << " peaks founds total, " << polym->GetN() << " used, and " << unusedMarkers->GetN() << " unused" << std::endl;
1114 std::cout << "Used: ";
1115 for(size_t m = 0; m < newX.size(); ++m) { std::cout << newX[m] << " - " << newY[m] << ";"; }
1116 std::cout << std::endl;
1117 std::cout << "Unused: ";
1118 for(size_t m = 0; m < unusedX.size(); ++m) { std::cout << unusedX[m] << " - " << unusedY[m] << ";"; }
1119 std::cout << std::endl;
1120 }
1122 std::cout << functions->GetEntries() << " functions in projection " << fProjection->GetTitle() << std::endl;
1123 for(auto* obj : *functions) {
1124 std::cout << obj->GetName() << "/" << obj->GetTitle() << " - " << obj->ClassName() << std::endl;
1125 }
1126 }
1127 }
1128 }
1129
1130 // change color of fit functions for unused peaks
1131 TIter iter(functions);
1132 TObject* item = nullptr;
1133 while((item = iter.Next()) != nullptr) {
1134 if(item->IsA() == TF1::Class() || item->IsA() == TGauss::Class()) { // if the item is a TF1 or TGauss we see if we can find the centroid in the map of used peaks
1135 double centroid = 0.;
1136 if(item->IsA() == TF1::Class()) {
1137 centroid = static_cast<TF1*>(item)->GetParameter(1);
1138 } else {
1139 centroid = static_cast<TGauss*>(item)->Centroid();
1140 }
1141 bool found = false;
1142 for(auto point : map) {
1143 if(TMath::Abs(centroid - point.first->Centroid()) < fSourceCalibration->Sigma()) {
1144 found = true;
1145 break;
1146 }
1147 }
1148 // remove the TF1 if it wasn't found in the map
1149 if(!found) {
1150 //functions->Remove(item);
1151 if(item->IsA() == TF1::Class()) {
1152 static_cast<TF1*>(item)->SetLineColor(kGray);
1153 } else {
1154 static_cast<TGauss*>(item)->GetFitFunction()->SetLineColor(kGray);
1155 }
1156 }
1157 }
1158 }
1159}
1160
1161void TSourceTab::ReplacePeak(const size_t& index, const double& channel)
1162{
1163 /// Replace the peak at the index with one centered at channel (calculated from an initial calibration and the source energy of this peak).
1164 /// This replaces the TGauss* in fPeaks, and updates the values of this index in fData and fFwhm.
1166 std::cout << DCYAN << __PRETTY_FUNCTION__ << ": index " << index << ", channel " << channel << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
1167 for(size_t i = 0; i < fPeaks.size() && i < fFits.size(); ++i) {
1168 std::cout << i << ": " << fPeaks[i]->Centroid() << " +- " << fPeaks[i]->CentroidErr() << " - " << fFits[i]->GetFitFunction()->GetParameter("centroid_0") << std::endl;
1169 }
1170 }
1171
1172 if(index >= fPeaks.size()) {
1173 std::cerr << "Trying to replace peak at index " << index << " but there are only " << fPeaks.size() << " peaks!" << std::endl;
1174 return;
1175 }
1176
1177 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << "Replacing old peak at index " << index << ", centroid " << fPeaks.at(index)->Centroid() << " with new peak at channel " << channel << std::endl; }
1178
1179 // calculate the fitting range (low to high) and adjust it so we ignore the old peak
1180 double range = fSourceCalibration->FitRange();
1181 double low = channel - range;
1182 double high = channel + range;
1183 if(fPeaks.at(index)->Centroid() < channel) {
1184 low = (fPeaks.at(index)->Centroid() + channel) / 2.;
1185 } else {
1186 high = (fPeaks.at(index)->Centroid() + channel) / 2.;
1187 }
1188 // if we aren't already redirecting the output, we redirect it to /dev/null and do a quiet fit, otherwise we do a normal fit since the output will be redirected to a file
1189 auto* pf = new TPeakFitter(low, high);
1190 auto* peak = new TGauss(channel);
1191 pf->AddPeak(peak);
1192 if(TSourceCalibration::LogFile().empty()) {
1193 TRedirect redirect("/dev/null");
1194 pf->Fit(fProjection, "retryfit");
1195 } else {
1196 pf->Fit(fProjection, "retryfit");
1197 }
1198 if(peak->Area() > 0) {
1199 if(Good(peak, low, high)) {
1201 std::cout << index << ". peak/fit out of " << fPeaks.size() << "/" << fFits.size() << " peaks/fits: " << std::flush << fPeaks.at(index)->Centroid() << "/" << fFits.at(index)->GetFitFunction()->GetParameter("centroid_0") << std::endl;
1202 std::cout << "Deleting " << index << ". peak " << fPeaks.at(index) << std::flush;
1203 }
1204 peak->SetLineColor(fPeaks.at(index)->GetLineColor());
1205 delete fPeaks.at(index);
1206 fPeaks[index] = peak;
1208 std::cout << " and setting it to " << fPeaks.at(index) << std::endl;
1209 }
1211 std::cout << "Deleting " << index << ". fit " << fFits.at(index) << std::flush;
1212 }
1213 pf->GetFitFunction()->SetLineColor(fFits.at(index)->GetFitFunction()->GetLineColor());
1214 delete fFits.at(index);
1215 fFits[index] = pf;
1217 std::cout << " and setting it to " << fFits.at(index) << std::endl;
1218 }
1219 UpdateFits();
1220 fData->SetPoint(index, peak->Centroid(), fData->GetPointY(index));
1221 fData->SetPointError(index, peak->CentroidErr(), fData->GetErrorY(index));
1222 fFwhm->SetPoint(index, peak->Centroid(), peak->FWHM());
1223 fFwhm->SetPointError(index, peak->CentroidErr(), peak->FWHMErr());
1225 std::cout << "Fitted peak " << low << " - " << high << " -> centroid " << peak->Centroid() << ", area " << peak->Area() << std::endl;
1226 }
1227 } else {
1229 std::cout << "We're ignoring bad fits, and this one has position " << peak->Centroid() << " +- " << peak->CentroidErr() << ", area " << peak->Area() << " +- " << peak->AreaErr() << std::endl;
1230 }
1231 }
1233 std::cout << "Ignoring peak at " << peak->Centroid() << " with negative area " << peak->Area() << std::endl;
1234 }
1235}
1236
1237void TSourceTab::RemovePoint(Int_t oldPoint)
1238{
1240 std::cout << DCYAN << __PRETTY_FUNCTION__ << ": oldPoint " << oldPoint << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
1241 }
1242 // first thing to do is check if the graph the point was removed from matches the id of this tab
1244 //auto erasedPeak = fPeaks.at(oldPoint);
1245 fPeaks.erase(fPeaks.begin() + oldPoint);
1246 auto* erasedFit = fFits.at(oldPoint);
1247 fFits.erase(fFits.begin() + oldPoint);
1248 fBadFits.push_back(erasedFit);
1249 fData->RemovePoint(oldPoint);
1250
1251 SetLineColors();
1252
1253 UpdateFits();
1254
1256}
1257
1259{
1260 /// This function sets the line colors of good and bad fits to alternating colors.
1261 /// Red and blue for the good fits, and grey and dark grey for the bad fits.
1262
1263 // set colors of good fits alternating to red or blue (peaks themselves are set to darker versions of those colors)
1264 for(size_t i = 0; i < fPeaks.size() && i < fFits.size(); ++i) {
1265 if(i % 2 == 0) {
1266 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << i << ": setting line colors to (dark) red for peak " << fPeaks[i] << ", fit function " << fFits[i]->GetFitFunction() << std::endl; }
1267 fPeaks[i]->SetLineColor(kRed + 2);
1268 fFits[i]->GetFitFunction()->SetLineColor(kRed);
1269 } else {
1270 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << i << ": setting line colors to (dark) blue for peak " << fPeaks[i] << ", fit function " << fFits[i]->GetFitFunction() << std::endl; }
1271 fPeaks[i]->SetLineColor(kBlue + 2);
1272 fFits[i]->GetFitFunction()->SetLineColor(kBlue);
1273 }
1274 }
1275
1276 // set colors of bad fits alterating to grey and dark grey
1277 for(size_t i = 0; i < fBadFits.size(); ++i) {
1278 if(i % 2 == 0) {
1279 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << i << ": setting line colors to grey for fit function " << fBadFits[i]->GetFitFunction() << std::endl; }
1280 fBadFits[i]->GetFitFunction()->SetLineColor(kGray);
1281 } else {
1282 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << i << ": setting line colors to dark grey for fit function " << fBadFits[i]->GetFitFunction() << std::endl; }
1283 fBadFits[i]->GetFitFunction()->SetLineColor(kGray + 2);
1284 }
1285 }
1286}
1287
1288void TSourceTab::Status(const char* status, int position)
1289{
1290 fSourceStatusBar->SetText(status, position);
1291 fSourceStatusBar->MapWindow();
1292 fSourceFrame->MapSubwindows();
1293 gVirtualX->Update();
1294}
1295
1297{
1298 std::cout << "Projection " << fProjection;
1299 if(fProjection != nullptr) {
1300 std::cout << " = " << fProjection->GetName() << "/" << fProjection->GetTitle();
1301 } else {
1302 std::cout << " is null";
1303 }
1304 std::cout << std::endl;
1305
1306 std::cout << "Data " << fData;
1307 if(fData != nullptr) {
1308 std::cout << " = " << fData->GetName() << "/" << fData->GetTitle();
1309 } else {
1310 std::cout << " is null";
1311 }
1312 std::cout << std::endl;
1313
1314 std::cout << "Fwhm " << fFwhm;
1315 if(fFwhm != nullptr) {
1316 std::cout << " = " << fFwhm->GetName() << "/" << fFwhm->GetTitle();
1317 } else {
1318 std::cout << " is null";
1319 }
1320 std::cout << std::endl;
1321
1322 std::cout << fFits.size() << " good fits:";
1323 for(size_t i = 0; i < fFits.size(); ++i) {
1324 std::cout << " " << std::setw(2) << i << " - " << std::setw(8) << fFits.at(i)->GetFitFunction()->GetParameter("centroid_0");
1325 }
1326 std::cout << std::endl;
1327
1328 std::cout << fBadFits.size() << " bad fits:";
1329 for(size_t i = 0; i < fBadFits.size(); ++i) {
1330 std::cout << " " << std::setw(2) << i << " - " << std::setw(8) << fBadFits.at(i)->GetFitFunction()->GetParameter("centroid_0");
1331 }
1332 std::cout << std::endl;
1333
1334 std::cout << fPeaks.size() << " peaks:";
1335 for(size_t i = 0; i < fPeaks.size(); ++i) {
1336 std::cout << " " << std::setw(2) << i << " - " << std::setw(8) << fPeaks.at(i)->Centroid();
1337 }
1338 std::cout << std::endl;
1339}
1340
1342{
1343 std::cout << this << ": TSourceTab " << fSourceName << " projection canvas:" << std::endl;
1344 fProjectionCanvas->GetCanvas()->Print();
1345 for(auto* obj : *(fProjectionCanvas->GetCanvas()->GetListOfPrimitives())) {
1346 std::cout << obj << ": " << obj->ClassName() << " called " << obj->GetName() << std::endl;
1347 obj->Print();
1348 }
1349 std::cout << "--------------------" << std::endl;
1350}
1351
1353{
1354 std::cout << "TSourceTab frame:" << std::endl;
1355 fSourceFrame->Print();
1356 std::cout << "TSourceTab canvas:" << std::endl;
1357 fProjectionCanvas->Print();
1358 std::cout << "TSourceTab status bar:" << std::endl;
1359 fSourceStatusBar->Print();
1360}
1361
1362//////////////////////////////////////// TChannelTab ////////////////////////////////////////
1363TChannelTab::TChannelTab(TSourceCalibration* sourceCal, std::vector<TNucleus*> nuclei, std::vector<GH1D*> projections, std::string name, TGCompositeFrame* frame, std::vector<std::vector<std::tuple<double, double, double, double>>> sourceEnergies, TGHProgressBar* progressBar)
1364 : fChannelFrame(frame), fSourceTab(new TGTab(frame, TSourceCalibration::PanelWidth(), TSourceCalibration::PanelHeight() + 2 * TSourceCalibration::StatusbarHeight())), fCanvasTab(new TGTab(frame, TSourceCalibration::PanelWidth(), TSourceCalibration::PanelHeight() + 2 * TSourceCalibration::StatusbarHeight())), fProgressBar(progressBar), fSourceCalibration(sourceCal), fNuclei(std::move(nuclei)), fProjections(std::move(projections)), fName(std::move(name)), fSourceEnergies(std::move(sourceEnergies))
1365{
1367 std::cout << DYELLOW << "========================================" << std::endl;
1368 std::cout << __PRETTY_FUNCTION__ << std::endl // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
1369 << " name = " << fName << ", fData = " << fData << std::endl;
1370 std::cout << "========================================" << std::endl;
1371 }
1372
1374
1376 std::cout << DYELLOW << "----------------------------------------" << std::endl;
1377 std::cout << __PRETTY_FUNCTION__ << std::endl // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
1378 << " channel " << fName << " done" << std::endl;
1379 std::cout << "----------------------------------------" << std::endl;
1380 }
1381}
1382
1384{
1385 delete fChannelFrame;
1386 delete fSourceTab;
1387 for(auto* channel : fSources) {
1388 delete channel;
1389 }
1390 delete fCanvasTab;
1391 delete fCalibrationFrame;
1392 delete fCalibrationCanvas;
1393 delete fResidualPad;
1394 delete fCalibrationPad;
1395 delete fChannelStatusBar;
1396 delete fLegend;
1397 delete fChi2Label;
1398 delete fFwhmFrame;
1399 delete fFwhmCanvas;
1400 // delete all fProjections and all fNuclei?
1401 delete fData;
1402 delete fFwhm;
1403}
1404
1406{
1407 {
1408 //std::unique_lock<std::mutex> graphicsLock(fSourceCalibration->GraphicsMutex());
1409 fChannelFrame->SetLayoutManager(new TGHorizontalLayout(fChannelFrame));
1410 }
1411
1412 fSources.resize(fNuclei.size(), nullptr);
1413 //if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << __PRETTY_FUNCTION__ << ": creating channels for bin 1 to " << fMatrix->GetNbinsX() << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
1414 for(size_t source = 0; source < fNuclei.size(); ++source) {
1415 CreateSourceTab(source);
1416 }
1417
1418 for(auto iter = fSources.begin(); iter != fSources.end(); ++iter) {
1419 if(*iter == nullptr) {
1420 fSources.erase(iter--); // erase iterator and then go to the element before this one (and then the loop goes to the next one)
1421 }
1422 }
1423
1424 {
1425 //std::unique_lock<std::mutex> graphicsLock(fSourceCalibration->GraphicsMutex());
1426 fChannelFrame->AddFrame(fSourceTab, new TGLayoutHints(kLHintsLeft | kLHintsTop | kLHintsExpandX | kLHintsExpandY, 2, 2, 2, 2));
1427
1428 fCalibrationFrame = fCanvasTab->AddTab("Calibration");
1429 fCalibrationFrame->SetLayoutManager(new TGVerticalLayout(fCalibrationFrame));
1430 fCalibrationCanvas = new TRootEmbeddedCanvas(Form("CalibrationCanvas%s", fName.c_str()), fCalibrationFrame, TSourceCalibration::PanelWidth(), TSourceCalibration::PanelHeight());
1431 fCalibrationFrame->AddFrame(fCalibrationCanvas, new TGLayoutHints(kLHintsTop | kLHintsExpandX | kLHintsExpandY, 2, 2, 2, 2));
1433 std::array<int, 3> parts = {25, 35, 40};
1434 fChannelStatusBar->SetParts(parts.data(), parts.size());
1435 fCalibrationFrame->AddFrame(fChannelStatusBar, new TGLayoutHints(kLHintsBottom | kLHintsExpandX, 2, 2, 2, 2));
1436 fFwhmFrame = fCanvasTab->AddTab("FWHM");
1437 fFwhmFrame->SetLayoutManager(new TGVerticalLayout(fFwhmFrame));
1438 fFwhmCanvas = new TRootEmbeddedCanvas(Form("FwhmCanvas%s", fName.c_str()), fFwhmFrame, TSourceCalibration::PanelWidth(), TSourceCalibration::PanelHeight());
1439 fFwhmFrame->AddFrame(fFwhmCanvas, new TGLayoutHints(kLHintsTop | kLHintsExpandX | kLHintsExpandY, 2, 2, 2, 2));
1440 //fChannelFrame->AddFrame(fCalibrationFrame, new TGLayoutHints(kLHintsRight | kLHintsBottom | kLHintsExpandX | kLHintsExpandY, 2, 2, 2, 2));
1441
1442 fInitFrame = fCanvasTab->AddTab("Initial");
1443 fInitFrame->SetLayoutManager(new TGVerticalLayout(fInitFrame));
1444 fInitCanvas = new TRootEmbeddedCanvas(Form("InitCanvas%s", fName.c_str()), fInitFrame, TSourceCalibration::PanelWidth(), TSourceCalibration::PanelHeight());
1445 fInitFrame->AddFrame(fInitCanvas, new TGLayoutHints(kLHintsTop | kLHintsExpandX | kLHintsExpandY, 2, 2, 2, 2));
1446
1447 fChannelFrame->AddFrame(fCanvasTab, new TGLayoutHints(kLHintsRight | kLHintsTop | kLHintsExpandX | kLHintsExpandY, 2, 2, 2, 2));
1448 }
1449
1450 UpdateData();
1451
1452 {
1453 //std::unique_lock<std::mutex> graphicsLock(fSourceCalibration->GraphicsMutex());
1454 fCalibrationCanvas->GetCanvas()->cd();
1455 fCalibrationPad = new TPad(Form("cal_%s", fName.c_str()), Form("calibration for %s", fName.c_str()), 0.2, 0., 1., 1.);
1456 fCalibrationPad->SetNumber(1);
1457 fCalibrationPad->Draw();
1458 fCalibrationPad->AddExec("zoom", "TChannelTab::ZoomY()");
1459
1460 fLegend = new TLegend(0.8, 0.3, 0.95, 0.3 + static_cast<double>(fNuclei.size()) * 0.05); // x1, y1, x2, y2
1461
1462 fCalibrationCanvas->GetCanvas()->cd();
1463 fResidualPad = new TPad(Form("res_%s", fName.c_str()), Form("residual for %s", fName.c_str()), 0.0, 0., 0.2, 1.);
1464 fResidualPad->SetNumber(2);
1465 fResidualPad->Draw();
1466 fResidualPad->AddExec("zoom", "TChannelTab::ZoomY()");
1467 }
1468 //Calibrate(); // also creates the residual and chi^2 label
1469
1470 // check whether we want to use this initial calibration to find more peaks
1471 //if(TSourceCalibration::UseCalibratedPeaks()) {
1472 // FindCalibratedPeaks();
1473 //}
1474
1475 // get the fwhm graphs and plot them
1476 UpdateFwhm();
1477 //{
1478 // std::unique_lock<std::mutex> graphicsLock(fSourceCalibration->GraphicsMutex());
1479 // fFwhmCanvas->GetCanvas()->cd();
1480 // fFwhm->DrawCalibration("*");
1481 // fLegend->Draw();
1482 //}
1483}
1484
1486{
1487 if(fProjections[source]->GetEntries() > 1000) {
1489 std::cout << DYELLOW << "Creating source tab " << source << ", using fSourceCalibration " << fSourceCalibration << std::endl;
1490 }
1491 TGCompositeFrame* tmpTab = nullptr;
1492 {
1493 //if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << "Creating unique lock on graphics mutex to add tab" << std::endl; }
1494 //std::unique_lock<std::mutex> graphicsLock(fSourceCalibration->GraphicsMutex());
1495 tmpTab = fSourceTab->AddTab(Form("%s_%s", fNuclei[source]->GetName(), fName.c_str()));
1496 //if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << "Releasing unique lock on graphics mutex after adding tab" << std::endl; }
1497 }
1499 std::cout << DYELLOW << "Creating source tab " << source << " = " << fProjections[source]->GetName() << ", tab name " << tmpTab->GetName() << std::endl;
1500 }
1501 fSources[source] = new TSourceTab(fSourceCalibration, this, tmpTab, fProjections[source], fNuclei[source]->GetName(), fSourceEnergies[source]);
1502 {
1503 //if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << "Creating unique lock on graphics mutex to increment status bar" << std::endl; }
1504 //std::unique_lock<std::mutex> graphicsLock(fSourceCalibration->GraphicsMutex());
1505 fProgressBar->Increment(1);
1506 //if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << "Releasing unique lock on graphics mutex after incrementing status bar" << std::endl; }
1507 }
1508 } else {
1509 fSources[source] = nullptr;
1511 std::cout << DYELLOW << "Skipping projection of source " << source << " = " << fProjections[source]->GetName() << ", only " << fProjections[source]->GetEntries() << " entries" << std::endl;
1512 }
1513 }
1515 std::cout << __PRETTY_FUNCTION__ << " source " << source << " done" << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
1516 }
1517}
1518
1520{
1521 //std::unique_lock<std::mutex> graphicsLock(fSourceCalibration->GraphicsMutex());
1522 fSourceTab->Connect("Selected(Int_t)", "TChannelTab", this, "SelectedTab(Int_t)");
1523 for(auto* source : fSources) {
1524 source->MakeConnections();
1525 }
1526 fCalibrationCanvas->GetCanvas()->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", "TChannelTab", this, "CalibrationStatus(Int_t,Int_t,Int_t,TObject*)");
1527 fCalibrationCanvas->Connect("ProcessedEvent(Event_t*)", "TChannelTab", this, "SelectCanvas(Event_t*)");
1528 if(fData != nullptr) {
1530 std::cout << DYELLOW << __PRETTY_FUNCTION__ << ": connecting fData " << fData << ":" << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
1531 fData->Print();
1532 }
1533 fData->Connect("RemovePoint(Int_t, Int_t)", "TChannelTab", this, "RemovePoint(Int_t, Int_t)");
1534 } else {
1535 std::cerr << DRED << "Failed to create connections for the data graph fData since it hasn't been created yet!" << std::endl;
1536 }
1537}
1538
1540{
1541 //std::unique_lock<std::mutex> graphicsLock(fSourceCalibration->GraphicsMutex());
1542 fSourceTab->Disconnect("Selected(Int_t)", this, "SelectedTab(Int_t)");
1543 for(auto* source : fSources) {
1544 source->Disconnect();
1545 }
1546 fCalibrationCanvas->GetCanvas()->Disconnect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", this, "CalibrationStatus(Int_t,Int_t,Int_t,TObject*)");
1547 fData->Disconnect("RemovePoint(Int_t, Int_t)", this, "RemovePoint(Int_t, Int_t)");
1548}
1549
1551{
1552 /// Simple function that enables and disables the previous and next buttons depending on which tab was selected.
1553 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << DYELLOW << __PRETTY_FUNCTION__ << ": id " << id << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
1554 fActiveSourceTab = id;
1555 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << "active source tab " << fActiveSourceTab << RESET_COLOR << std::endl; }
1556 if(fSources[fActiveSourceTab] == nullptr) {
1557 std::cout << "Failed to get fSources[" << fActiveSourceTab << "]" << RESET_COLOR << std::endl;
1558 return;
1559 }
1560 if(fSources[fActiveSourceTab]->ProjectionCanvas() == nullptr) {
1561 std::cout << "Failed to get fSources[" << fActiveSourceTab << "]->ProjectionCanvas()" << RESET_COLOR << std::endl;
1562 return;
1563 }
1564 if(fSources[fActiveSourceTab]->ProjectionCanvas()->GetCanvas() == nullptr) {
1565 std::cout << "Failed to get fSources[" << fActiveSourceTab << "]->ProjectionCanvas()->GetCanvas()" << RESET_COLOR << std::endl;
1566 return;
1567 }
1568 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << DGREEN << "Changing gpad from \"" << gPad->GetName() << "\" to \""; }
1569 gPad = fSources[fActiveSourceTab]->ProjectionCanvas()->GetCanvas();
1570 fSources[fActiveSourceTab]->ProjectionCanvas()->GetCanvas()->cd();
1571 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << gPad->GetName() << "\"" << RESET_COLOR << std::endl; }
1572}
1573
1574void TChannelTab::SelectCanvas(Event_t* event)
1575{
1576 if(event->fType == kEnterNotify) {
1578 std::cout << "Entered channel tab => changing gPad from " << gPad->GetName();
1579 }
1580 gPad = fCalibrationCanvas->GetCanvas();
1582 std::cout << " to " << gPad->GetName() << std::endl;
1583 }
1584 }
1585}
1586
1587void TChannelTab::RemovePoint(Int_t oldGraph, Int_t oldPoint)
1588{
1590 std::cout << DCYAN << __PRETTY_FUNCTION__ << ": oldGraph " << oldGraph << ", oldPoint " << oldPoint << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
1591 }
1592 if(0 <= oldGraph && oldGraph < static_cast<int>(fSources.size())) {
1593 if(oldPoint >= 0) {
1594 fSources[oldGraph]->RemovePoint(oldPoint);
1595 return;
1596 }
1597 std::cout << "Can't remove negative point " << oldPoint << " from graph " << oldGraph << std::endl;
1598 }
1599 std::cout << "Graph the point was removed from was " << oldGraph << ", but we only have " << fSources.size() << " source tabs?" << std::endl;
1600}
1601
1602void TChannelTab::CalibrationStatus(Int_t event, Int_t px, Int_t py, TObject* selected)
1603{
1604 // looks like 51 is hovering over the object, 52 is moving the cursor over the object, and 53 is moving the cursort away from the object
1605 // kButton1 = left mouse button, kButton2 = right mouse button
1606 // enum EEventType {
1607 // kNoEvent = 0,
1608 // kButton1Down = 1, kButton2Down = 2, kButton3Down = 3, kKeyDown = 4,
1609 // kWheelUp = 5, kWheelDown = 6, kButton1Shift = 7, kButton1ShiftMotion = 8,
1610 // kButton1Up = 11, kButton2Up = 12, kButton3Up = 13, kKeyUp = 14,
1611 // kButton1Motion = 21, kButton2Motion = 22, kButton3Motion = 23, kKeyPress = 24,
1612 // kArrowKeyPress = 25, kArrowKeyRelease = 26,
1613 // kButton1Locate = 41, kButton2Locate = 42, kButton3Locate = 43, kESC = 27,
1614 // kMouseMotion = 51, kMouseEnter = 52, kMouseLeave = 53,
1615 // kButton1Double = 61, kButton2Double = 62, kButton3Double = 63
1616 //};
1617 fChannelStatusBar->SetText(selected->GetName(), 0);
1618 if(event == kKeyPress) {
1619 fChannelStatusBar->SetText(Form("%c", static_cast<char>(px)), 1);
1620 } else {
1621 auto* canvas = fCalibrationCanvas->GetCanvas();
1622 if(canvas == nullptr) {
1623 fChannelStatusBar->SetText("canvas nullptr");
1624 return;
1625 }
1626 auto* pad = canvas->GetSelectedPad();
1627 if(pad == nullptr) {
1628 fChannelStatusBar->SetText("pad nullptr");
1629 return;
1630 }
1631 //if(TSourceCalibration::VerboseLevel() > EVerbosity::kLoops) {
1632 // std::cout << "px " << px << ", py " << py << "; gPad: px " << gPad->AbsPixeltoX(px) << ", py " << gPad->AbsPixeltoY(py) << "; fCalibrationCanvas: px " << fCalibrationCanvas->GetCanvas()->AbsPixeltoX(px) << ", py " << fCalibrationCanvas->GetCanvas()->AbsPixeltoY(py) << "; fResidualPad: px " << fResidualPad->AbsPixeltoX(px) << ", py " << fResidualPad->AbsPixeltoY(py) << "; fCalibrationPad: px " << fCalibrationPad->AbsPixeltoX(px) << ", py " << fCalibrationPad->AbsPixeltoY(py) << "; pad: px " << pad->AbsPixeltoX(px) << ", py " << pad->AbsPixeltoY(py) << std::endl;
1633 //}
1634
1635 fChannelStatusBar->SetText(Form("x = %f, y = %f", pad->AbsPixeltoX(px), pad->AbsPixeltoY(py)), 1);
1636 }
1637}
1638
1640{
1641 /// Copy data from all sources into one graph (which we use for the calib && source < fSources.size()ration).
1642 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << DYELLOW << __PRETTY_FUNCTION__ << " fData = " << fData << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
1643 if(fData == nullptr) {
1644 fData = new TCalibrationGraphSet("ADC channel", "Energy [keV]");
1645 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << __PRETTY_FUNCTION__ << " created fData = " << fData << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
1646 }
1647 fData->Clear();
1648
1649 fData->SetName(Form("data%s", fName.c_str()));
1651 std::cout << "set name of fData using " << fName.c_str() << ": " << fData->GetName() << std::endl;
1652 std::cout << "fData " << fData << ": " << (fData != nullptr ? fData->GetN() : -2) << " data points after creation" << std::endl;
1653 std::cout << "Looping over " << fNuclei.size() << " sources, fSources.size() = " << fSources.size() << std::endl;
1654 }
1655 for(size_t source = 0; source < fNuclei.size() && source < fSources.size(); ++source) {
1656 if(fSources[source]->Data() != nullptr && fSources[source]->Data()->GetN() > 0) {
1657 int index = fData->Add(fSources[source]->Data(), fNuclei[source]->GetName());
1658 if(index >= 0) {
1659 fData->SetLineColor(index, static_cast<Color_t>(source + 1)); //+1 for the color so that we start with 1 = black instead of 0 = white
1660 fData->SetMarkerColor(index, static_cast<Color_t>(source + 1));
1661 }
1662 }
1663 }
1665 std::cout << "fData " << fData << ": " << (fData != nullptr ? fData->GetN() : -1) << " data points after adding " << fNuclei.size() << " graphs" << std::endl;
1666 fData->Print();
1667 }
1668}
1669
1671{
1672 /// Copy data from all sources into one graph (which we use for the calib && source < fSources.size()ration).
1673 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << DYELLOW << __PRETTY_FUNCTION__ << " fFwhm = " << fFwhm << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
1674 if(fFwhm == nullptr) {
1675 fFwhm = new TCalibrationGraphSet("ADC channel", "FWHM [ADC channel]");
1676 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << __PRETTY_FUNCTION__ << " created fFwhm = " << fFwhm << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
1677 }
1678 fFwhm->Clear();
1679
1680 fFwhm->SetName(Form("fwhm%s", fName.c_str()));
1682 std::cout << "set name of fFwhm using " << fName.c_str() << ": " << fFwhm->GetName() << std::endl;
1683 std::cout << "fFwhm " << fFwhm << ": " << (fFwhm != nullptr ? fFwhm->GetN() : -2) << " data points after creation" << std::endl;
1684 std::cout << "Looping over " << fNuclei.size() << " sources, fSources.size() = " << fSources.size() << std::endl;
1685 }
1686 for(size_t source = 0; source < fNuclei.size() && source < fSources.size(); ++source) {
1687 if(fSources[source]->Fwhm() != nullptr && fSources[source]->Fwhm()->GetN() > 0) {
1688 int index = fFwhm->Add(fSources[source]->Fwhm(), fNuclei[source]->GetName());
1689 if(index >= 0) {
1690 fFwhm->SetLineColor(index, static_cast<Color_t>(source + 1)); //+1 for the color so that we start with 1 = black instead of 0 = white
1691 fFwhm->SetMarkerColor(index, static_cast<Color_t>(source + 1));
1692 }
1693 }
1694 }
1696 std::cout << "fFwhm " << fFwhm << ": " << (fFwhm != nullptr ? fFwhm->GetN() : -1) << " data points after adding " << fNuclei.size() << " graphs" << std::endl;
1697 fFwhm->Print();
1698 }
1699}
1700
1702{
1703 /// Uses the statis TSourceCalibration::Degree for the degree of the polynomial used to calibrate the data.
1705}
1706
1708{
1709 /// This function fit's the final data of the given channel.
1710 /// It requires all other elements to have been created already.
1711 if(TSourceCalibration::VerboseLevel() >= EVerbosity::kSubroutines) { std::cout << __PRETTY_FUNCTION__ << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
1712 TF1* calibration = new TF1("fitfunction", ::Polynomial, 0., 10000., degree + 2);
1713 calibration->FixParameter(0, degree);
1715 std::cout << DYELLOW << "fData " << fData << ": ";
1716 if(fData != nullptr) {
1717 std::cout << fData->GetN() << " data points being fit, ";
1718 }
1719 double min = 0.;
1720 double max = 0.;
1721 calibration->GetRange(min, max);
1722 std::cout << "range " << min << " - " << max << std::endl;
1723 }
1724 fData->Fit(calibration, "Q");
1725 TString text = Form("%.6f + %.6f*x", calibration->GetParameter(1), calibration->GetParameter(2));
1726 for(int i = 2; i <= degree; ++i) {
1727 text.Append(Form(" + %.6f*x^%d", calibration->GetParameter(i + 1), i));
1728 }
1729 fChannelStatusBar->SetText(text.Data(), 2);
1730 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << "re-calculating residuals and clearing fields" << std::endl; }
1731 // re-calculate the residuals
1732 fData->SetResidual(true);
1733
1734 fLegend->Clear();
1735 fCalibrationPad->cd();
1737 fLegend->Draw();
1738 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << "set chi2 label" << std::endl; }
1739 // calculate the corners of the chi^2 label from the minimum and maximum x/y-values of the graph
1740 // we position it in the top left corner about 50% of the width and 10% of the height of the graph
1741 if(fChi2Label == nullptr) {
1742 double left = fData->GetMinimumX();
1743 double right = left + (fData->GetMaximumX() - left) * 0.5;
1744 double top = fData->GetMaximumY();
1745 double bottom = top - (top - fData->GetMinimumY()) * 0.1;
1746
1747 fChi2Label = new TPaveText(left, bottom, right, top);
1749 fChi2Label->ConvertNDCtoPad();
1750 std::cout << "fChi2Label created " << fChi2Label << " (" << left << " - " << right << ", " << bottom << " - " << top << ", from " << fData->GetMinimumX() << "-" << fData->GetMaximumX() << ", " << fData->GetMinimumY() << "-" << fData->GetMaximumY() << ") on gPad " << gPad->GetName() << " coordinates: " << fChi2Label->GetX1() << "/" << fChi2Label->GetX1NDC() << " - " << fChi2Label->GetX2() << "/" << fChi2Label->GetX2NDC() << ", " << fChi2Label->GetY1() << "/" << fChi2Label->GetY1NDC() << " - " << fChi2Label->GetY2() << "/" << fChi2Label->GetY2NDC() << std::endl;
1751 fData->Print("e");
1752 }
1753 } else {
1754 fChi2Label->Clear();
1755 }
1756 fChi2Label->AddText(Form("#chi^{2}/NDF = %f", calibration->GetChisquare() / calibration->GetNDF()));
1757 fChi2Label->SetFillColor(kWhite);
1758 fChi2Label->Draw();
1760 std::cout << "fChi2Label set to:" << std::endl;
1761 fChi2Label->GetListOfLines()->Print();
1762 std::cout << "Text size " << fChi2Label->GetTextSize() << std::endl;
1763 }
1764 if(fInfoLabel == nullptr) {
1765 fInfoLabel = new TPaveText(0, 0.95, 1., 1., "NDCNB");
1767 fInfoLabel->ConvertNDCtoPad();
1768 std::cout << "fInfoLabel created " << fInfoLabel << " (" << fInfoLabel->GetX1() << "/" << fInfoLabel->GetX1NDC() << " - " << fInfoLabel->GetX2() << "/" << fInfoLabel->GetX2NDC() << ", " << fInfoLabel->GetY1() << "/" << fInfoLabel->GetY1NDC() << " - " << fInfoLabel->GetY2() << "/" << fInfoLabel->GetY2NDC() << ", from 0, 0.95, 1., 1.) on gPad " << gPad->GetName() << std::endl;
1769 }
1770 } else {
1771 fInfoLabel->Clear();
1773 fInfoLabel->ConvertNDCtoPad();
1774 std::cout << "fInfoLabel using existing " << fInfoLabel << " (" << fInfoLabel->GetX1() << "/" << fInfoLabel->GetX1NDC() << " - " << fInfoLabel->GetX2() << "/" << fInfoLabel->GetX2NDC() << ", " << fInfoLabel->GetY1() << "/" << fInfoLabel->GetY1NDC() << " - " << fInfoLabel->GetY2() << "/" << fInfoLabel->GetY2NDC() << ") on gPad " << gPad->GetName() << std::endl;
1775 }
1776 }
1777 fInfoLabel->AddText(Form("Using %d peaks from %.0f to %.0f keV", fData->TotalGraph()->GetN(), fData->GetMinimumY(), fData->GetMaximumY()));
1778 fInfoLabel->SetFillColor(kWhite);
1779 fInfoLabel->Draw();
1780
1781 fResidualPad->cd();
1782 fData->DrawResidual("*r0");
1783
1784 fCalibrationCanvas->GetCanvas()->Modified();
1785
1786 fData->FitFunction()->SetRange(0., 10000.);
1788 double min = 0.;
1789 double max = 0.;
1790 calibration->GetRange(min, max);
1791 std::cout << DYELLOW << "calibration has range " << min << " - " << max;
1792 if(fData->FitFunction() != nullptr) {
1793 fData->FitFunction()->GetRange(min, max);
1794 std::cout << ", fit function has range " << min << " - " << max;
1795 }
1796 std::cout << std::endl;
1797 }
1798 delete calibration;
1799
1800 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << __PRETTY_FUNCTION__ << " done" << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
1801}
1802
1803void TChannelTab::RecursiveRemove(const double& maxResidual)
1804{
1805 /// This function goes through the residuals of the calibration and recursively removes each point that is above the maxResidual parameter.
1806 /// The function requires an initial calibration to have residuals available.
1807 /// After each point that is removed, it re-calibrates the data, producing new residuals.
1808
1809 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << DYELLOW << __PRETTY_FUNCTION__ << ", max. residual " << maxResidual << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
1810
1811 if(fData == nullptr) {
1812 std::cout << "No longer have a data for channel " << Name() << "?" << std::endl;
1813 return;
1814 }
1815
1816 if(fData->TotalResidualGraph() == nullptr) {
1817 std::cout << "No longer have a total residual graph for channel " << Name() << ". Maybe removed all points? Total graph has " << fData->TotalGraph()->GetN() << " points." << std::endl;
1818 return;
1819 }
1820
1821 if(fData->TotalResidualGraph()->GetN() == 0) {
1822 std::cout << "Recursively removed all data points for channel " << Name() << std::endl;
1823 return;
1824 }
1825
1826 // get the maximum residual
1827 auto* maxElement = std::max_element(fData->TotalResidualGraph()->GetX(), fData->TotalResidualGraph()->GetX() + fData->TotalResidualGraph()->GetN(), [](double a, double b) { return std::abs(a) < std::abs(b); });
1828 auto index = std::distance(fData->TotalResidualGraph()->GetX(), maxElement);
1829
1830 if(std::abs(*maxElement) < maxResidual) {
1831 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << DYELLOW << "largest residual is " << *maxElement << " at position " << index << " and is smaller than max. residual " << maxResidual << ": stopping recursive remove" << std::endl; }
1832 return;
1833 }
1834
1835 // remove the index with the maximum residual
1836 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << DYELLOW << "std::abs(" << *maxElement << ") > " << maxResidual << ": removing index " << index << "/" << fData->TotalResidualGraph()->GetN() << std::endl; }
1837 fData->RemovePoint(index);
1838
1839 // update the data, and re-calibrate
1840 UpdateData();
1841 UpdateFwhm();
1842 Calibrate();
1843
1844 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << "RecursiveRemove(" << maxResidual << ") being called again" << RESET_COLOR << std::endl; }
1845 RecursiveRemove(maxResidual);
1846}
1847
1849{
1850 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << DYELLOW << __PRETTY_FUNCTION__ << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
1851 // write the actual parameters of the fit
1852 std::stringstream str;
1853 str << std::hex << fName;
1854 int address = 0;
1855 str >> address;
1856 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << "Got address " << hex(address, 4) << " from name " << fName << std::endl; }
1857 TF1* calibration = fData->FitFunction();
1858 if(calibration == nullptr) {
1859 std::cout << "Failed to find calibration in fData" << std::endl;
1860 fData->TotalGraph()->GetListOfFunctions()->Print();
1861 return;
1862 }
1863 std::vector<Float_t> parameters;
1864 for(int i = 0; i <= calibration->GetParameter(0); ++i) {
1865 parameters.push_back(static_cast<Float_t>(calibration->GetParameter(i + 1)));
1866 }
1867 TChannel* channel = TChannel::GetChannel(address, false);
1868 if(channel == nullptr) {
1869 std::cerr << "Failed to get channel for address " << hex(address, 4) << std::endl;
1870 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << __PRETTY_FUNCTION__ << ": done" << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
1871 return;
1872 }
1873 channel->SetENGCoefficients(parameters);
1874 channel->DestroyEnergyNonlinearity();
1876 double* x = fData->GetX();
1877 double* y = fData->GetY();
1878 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << "Going to add " << fData->GetN() << " points to nonlinearity graph" << std::endl; }
1879 for(int i = 0; i < fData->GetN(); ++i) {
1880 // nonlinearity expects y to be the true source energy minus the calibrated energy of the peak
1881 // the source energy is y, the peak is x
1882 channel->AddEnergyNonlinearityPoint(y[i], y[i] - calibration->Eval(x[i]));
1883 if(TSourceCalibration::VerboseLevel() > EVerbosity::kLoops) { std::cout << "Added " << channel->GetEnergyNonlinearity().GetN() << ". point " << y[i] << ", " << y[i] - calibration->Eval(x[i]) << " = " << y[i] << " - " << calibration->Eval(x[i]) << std::endl; }
1884 }
1885 }
1886}
1887
1888void TChannelTab::Initialize(const bool& force)
1889{
1890 // loop through sources and get rough calibration for each
1891 for(auto* source : fSources) {
1893 std::cout << DYELLOW << "Rough calibration in source tab:" << std::endl;
1894 source->Print();
1895 }
1896 source->InitialCalibration(force);
1897 fProgressBar->Increment(1);
1898 }
1899 // get this rough calibration data and calibrate with linear calibration
1900 UpdateData();
1901 Calibrate(1);
1902 // copy the data to the initial data and draw it
1903 fInit = static_cast<TCalibrationGraphSet*>(fData->Clone(Form("init%s", fName.c_str())));
1905 std::cout << "Cloned data " << fData->GetName() << " to " << fInit->GetName() << ":" << std::endl;
1906 fInit->Print();
1907 }
1908 auto* oldPad = gPad;
1909 fInitCanvas->GetCanvas()->cd();
1910 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << "Switched from " << oldPad->GetName() << " to " << gPad->GetName() << std::endl; }
1912 // calculate the corners of the label from the minimum and maximum x/y-values of the graph
1913 // we position it in the top left corner about 50% of the width and 10% of the height of the graph
1914 if(fCalLabel == nullptr) {
1915 double left = fInit->GetMinimumX();
1916 double right = left + (fInit->GetMaximumX() - left) * 0.5;
1917 double top = fInit->GetMaximumY();
1918 double bottom = top - (top - fInit->GetMinimumY()) * 0.1;
1919
1920 fCalLabel = new TPaveText(left, bottom, right, top);
1922 std::cout << "fCalLabel created " << fCalLabel << " (" << left << " - " << right << ", " << bottom << " - " << top << ", from " << fInit->GetMinimumX() << "-" << fInit->GetMaximumX() << ", " << fInit->GetMinimumY() << "-" << fInit->GetMaximumY() << ") on gPad " << gPad->GetName() << std::endl;
1923 }
1924 } else {
1925 fCalLabel->Clear();
1926 }
1927 std::stringstream str;
1928 if(fInit->FitFunction()->GetNpar() == 3) {
1929 str << "Fit function: " << fInit->FitFunction()->GetParameter(1) << " + " << fInit->FitFunction()->GetParameter(2) << " x";
1930 } else {
1931 str << "Unknown fit function with " << fInit->FitFunction()->GetNpar() << " parameters";
1932 }
1934 std::cout << "Created label text \"" << str.str() << "\"" << std::endl;
1935 }
1936
1937 fCalLabel->AddText(str.str().c_str());
1938 fCalLabel->SetFillColor(kWhite);
1939 fCalLabel->Draw();
1940 oldPad->cd();
1941 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << "Switched back to " << gPad->GetName() << std::endl; }
1942}
1943
1945{
1947 std::cout << DYELLOW << "Finding calibrated peaks in source tab " << fActiveSourceTab << ". Old: " << fSourceTab->GetCurrent() << " = " << fSources[fSourceTab->GetCurrent()] << ", current tab element " << fSourceTab->GetCurrentTab() << " is enabled = " << fSourceTab->GetCurrentTab()->IsEnabled() << std::endl;
1948 fSourceTab->Print();
1949 for(int tab = 0; tab < fSourceTab->GetNumberOfTabs(); ++tab) {
1950 std::cout << "Tab " << tab << " = " << fSourceTab->GetTabTab(tab) << " = " << fSourceTab->GetTabTab(tab)->GetText()->GetString() << (fSourceTab->GetTabTab(tab)->IsActive() ? " is active" : " is inactive") << (fSourceTab->GetTabTab(tab)->IsEnabled() ? " is enabled" : " is not enabled") << std::endl;
1951 }
1952 }
1953 fSources[fActiveSourceTab]->FindCalibratedPeaks(fData->FitFunction());
1954 UpdateData();
1955 UpdateFwhm();
1957 std::cout << DYELLOW << "Done finding calibrated peaks in source tab " << fActiveSourceTab << ". Old: " << fSourceTab->GetCurrent() << " = " << fSources[fSourceTab->GetCurrent()] << ", current tab element " << fSourceTab->GetCurrentTab() << " is enabled = " << fSourceTab->GetCurrentTab()->IsEnabled() << std::endl;
1958 }
1959}
1960
1962{
1964 std::cout << DYELLOW << "Finding calibrated peaks in all source tabs" << std::endl;
1965 }
1966 for(auto* source : fSources) {
1968 std::cout << DYELLOW << "Finding calibrated peaks in source tab " << source->SourceName() << std::endl;
1969 }
1970 source->FindCalibratedPeaks(fData->FitFunction());
1971 }
1972 UpdateData();
1973 UpdateFwhm();
1975 std::cout << DYELLOW << "Done finding calibrated peaks in all source tabs" << std::endl;
1976 }
1977}
1978
1980{
1981 fFwhmCanvas->GetCanvas()->cd();
1982 fFwhm->DrawCalibration("*");
1983 fLegend->Draw();
1984 for(auto* source : fSources) {
1985 source->Draw();
1986 }
1987}
1988
1989void TChannelTab::Write(TFile* output)
1990{
1991 /// Write graphs to output.
1992 output->cd();
1993 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << DYELLOW << "writing " << fName << std::endl; }
1994 fData->Write(Form("calGraph_%s", fName.c_str()), TObject::kOverwrite);
1995 fFwhm->Write(Form("fwhm_%s", fName.c_str()), TObject::kOverwrite);
1996 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << DYELLOW << "wrote data " << fName << std::endl; }
1997}
1998
2000{
2001 /// finds corresponding graph (res_xxx for cal_xxx and vice versa) and sets it's x-axis range to be the same as the selected graphs
2002 // find the histogram in the current pad
2003 TH1* hist1 = nullptr;
2004 for(const auto&& obj : *(gPad->GetListOfPrimitives())) {
2005 if(obj->InheritsFrom(TGraph::Class())) {
2006 hist1 = static_cast<TGraph*>(obj)->GetHistogram();
2007 break;
2008 }
2009 }
2010 if(hist1 == nullptr) {
2011 std::cout << "ZoomX - Failed to find histogram for pad " << gPad->GetName() << std::endl;
2012 if(TSourceCalibration::VerboseLevel() > EVerbosity::kLoops) { gPad->GetListOfPrimitives()->Print(); }
2013 return;
2014 }
2015
2016 // extract base name and channel from pad name
2017 std::string padName = gPad->GetName();
2018 std::string baseName = padName.substr(0, 3);
2019 std::string channelName = padName.substr(4);
2020
2021 // find the corresponding partner pad to the selected one
2022 std::string newName;
2023 if(baseName == "cal") {
2024 newName = "res_" + channelName;
2025 } else if(baseName == "res") {
2026 newName = "cal_" + channelName;
2027 } else {
2028 std::cout << "Unknown combination of pad " << gPad->GetName() << " and histogram " << hist1->GetName() << std::endl;
2029 return;
2030 }
2031 auto* newObj = gPad->GetCanvas()->FindObject(newName.c_str());
2032 if(newObj == nullptr) {
2033 std::cout << "Failed to find " << newName << std::endl;
2034 return;
2035 }
2036
2037 if(newObj->IsA() != TPad::Class()) {
2038 std::cout << newObj << " = " << newObj->GetName() << ", " << newObj->GetTitle() << " is a " << newObj->ClassName() << " not a TPad" << std::endl;
2039 return;
2040 }
2041 auto* pad2 = static_cast<TPad*>(newObj);
2042 if(pad2 == nullptr) {
2043 std::cout << "Failed to find partner pad " << newName << std::endl;
2044 return;
2045 }
2046 // find the histogram in the partner pad
2047 TH1* hist2 = nullptr;
2048 for(const auto&& obj : *(pad2->GetListOfPrimitives())) {
2049 if(obj->InheritsFrom(TGraph::Class())) {
2050 hist2 = static_cast<TGraph*>(obj)->GetHistogram();
2051 break;
2052 }
2053 }
2054 if(hist2 == nullptr) {
2055 std::cout << "ZoomX - Failed to find histogram for partner pad " << newName << std::endl;
2056 if(TSourceCalibration::VerboseLevel() > EVerbosity::kLoops) { pad2->GetListOfPrimitives()->Print(); }
2057 return;
2058 }
2059
2060 TAxis* axis1 = hist1->GetXaxis();
2061 Int_t binmin = axis1->GetFirst();
2062 Int_t binmax = axis1->GetLast();
2063 Double_t xmin = axis1->GetBinLowEdge(binmin);
2064 Double_t xmax = axis1->GetBinLowEdge(binmax);
2065 TAxis* axis2 = hist2->GetXaxis();
2066 Int_t newmin = axis2->FindBin(xmin);
2067 Int_t newmax = axis2->FindBin(xmax);
2068 axis2->SetRange(newmin, newmax);
2069 pad2->Modified();
2070 pad2->Update();
2071}
2072
2074{
2075 /// finds corresponding graph (res_xxx for cal_xxx and vice versa) and sets it's y-axis range to be the same as the selected graphs
2076 // find the histogram in the current pad
2077 TH1* hist1 = nullptr;
2078 for(const auto&& obj : *(gPad->GetListOfPrimitives())) {
2079 if(obj->InheritsFrom(TGraph::Class())) {
2080 hist1 = static_cast<TGraph*>(obj)->GetHistogram();
2081 break;
2082 }
2083 }
2084 if(hist1 == nullptr) {
2085 std::cout << "ZoomY - Failed to find histogram for pad " << gPad->GetName() << std::endl;
2086 if(TSourceCalibration::VerboseLevel() > EVerbosity::kLoops) { gPad->GetListOfPrimitives()->Print(); }
2087 return;
2088 }
2089
2090 // extract base name and channel from pad name
2091 std::string padName = gPad->GetName();
2092 std::string baseName = padName.substr(0, 3);
2093 std::string channelName = padName.substr(4);
2094
2095 // find the corresponding partner pad to the selected one
2096 std::string newName;
2097 if(baseName == "cal") {
2098 newName = "res_" + channelName;
2099 } else if(baseName == "res") {
2100 newName = "cal_" + channelName;
2101 } else {
2102 std::cout << "Unknown combination of pad " << gPad->GetName() << " and histogram " << hist1->GetName() << std::endl;
2103 return;
2104 }
2105 auto* newObj = gPad->GetCanvas()->FindObject(newName.c_str());
2106 if(newObj == nullptr) {
2107 std::cout << "Failed to find " << newName << std::endl;
2108 return;
2109 }
2110
2111 if(newObj->IsA() != TPad::Class()) {
2112 std::cout << newObj << " = " << newObj->GetName() << ", " << newObj->GetTitle() << " is a " << newObj->ClassName() << " not a TPad" << std::endl;
2113 return;
2114 }
2115 auto* pad2 = static_cast<TPad*>(newObj);
2116 if(pad2 == nullptr) {
2117 std::cout << "Failed to find partner pad " << newName << std::endl;
2118 return;
2119 }
2120 // find the histogram in the partner pad
2121 TH1* hist2 = nullptr;
2122 for(const auto&& obj : *(pad2->GetListOfPrimitives())) {
2123 if(obj->InheritsFrom(TGraph::Class())) {
2124 hist2 = static_cast<TGraph*>(obj)->GetHistogram();
2125 break;
2126 }
2127 }
2128 if(hist2 == nullptr) {
2129 std::cout << "ZoomY - Failed to find histogram for partner pad " << newName << std::endl;
2130 if(TSourceCalibration::VerboseLevel() > EVerbosity::kLoops) { pad2->GetListOfPrimitives()->Print(); }
2131 return;
2132 }
2133
2134 hist2->SetMinimum(hist1->GetMinimum());
2135 hist2->SetMaximum(hist1->GetMaximum());
2136
2137 // update residual line
2138 TLine* resLine = nullptr;
2139 for(const auto&& obj : *(gPad->GetListOfPrimitives())) {
2140 if(obj->IsA() == TLine::Class()) {
2141 if(resLine != nullptr) { std::cout << "ZoomY found residual line, old was " << resLine << std::endl; }
2142 resLine = static_cast<TLine*>(obj);
2143 }
2144 }
2145 for(const auto&& obj : *(pad2->GetListOfPrimitives())) {
2146 if(obj->IsA() == TLine::Class()) {
2147 if(resLine != nullptr) { std::cout << "ZoomY found residual line, old was " << resLine << std::endl; }
2148 resLine = static_cast<TLine*>(obj);
2149 }
2150 }
2151
2152 // update the y-limits of the line at zero residual
2153 if(resLine != nullptr) {
2154 resLine->SetY1(hist1->GetMinimum());
2155 resLine->SetY2(hist1->GetMaximum());
2156 }
2157
2158 pad2->Modified();
2159 pad2->Update();
2160}
2161
2163{
2164 std::cout << "----------------------------------------" << std::endl;
2165 std::cout << this << ": TChannelTab " << Name() << " calibration canvas:" << std::endl;
2166 fCalibrationCanvas->GetCanvas()->Print();
2167 for(auto* obj : *(fCalibrationCanvas->GetCanvas()->GetListOfPrimitives())) {
2168 obj->Print();
2169 }
2170 std::cout << "TChannelTab fwhm canvas:" << std::endl;
2171 fFwhmCanvas->GetCanvas()->Print();
2172 for(auto* obj : *(fFwhmCanvas->GetCanvas()->GetListOfPrimitives())) {
2173 obj->Print();
2174 }
2175 for(auto* sourceTab : fSources) {
2176 sourceTab->PrintCanvases();
2177 }
2178 std::cout << "----------------------------------------" << std::endl;
2179}
2180
2182{
2183 std::cout << "TChannelTab frame:" << std::endl;
2184 fChannelFrame->Print();
2185 std::cout << "TChannelTab canvas:" << std::endl;
2186 fCalibrationCanvas->Print();
2187 std::cout << "TChannelTab status bar:" << std::endl;
2188 //fChannelStatusBar->Print();
2189 for(auto* sourceTab : fSources) {
2190 sourceTab->PrintLayout();
2191 }
2192}
2193
2194//////////////////////////////////////// TSourceCalibration ////////////////////////////////////////
2206bool TSourceCalibration::fFast = false;
2210std::vector<int> TSourceCalibration::fBadBins;
2211
2212TSourceCalibration::TSourceCalibration(double sigma, double threshold, int degree, int count...)
2213 : TGMainFrame(nullptr, 2 * fPanelWidth, fPanelHeight + 2 * fStatusbarHeight), fDefaultSigma(sigma), fThreshold(threshold), fDefaultDegree(degree)
2214{
2215 TH1::AddDirectory(false); // turns off warnings about multiple histograms with the same name because ROOT doesn't manage them anymore
2216
2217 va_list args;
2218 va_start(args, count); // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2219 for(int i = 0; i < count; ++i) {
2220 fMatrices.push_back(va_arg(args, TH2*)); // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2221 if(fMatrices.back() == nullptr) {
2222 std::cout << "Error, got nullptr as matrix input?" << std::endl;
2223 fMatrices.pop_back();
2224 }
2225 }
2226 va_end(args); // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2228 std::cout << DGREEN << __PRETTY_FUNCTION__ << ": verbose level " << static_cast<std::underlying_type_t<EVerbosity>>(fVerboseLevel) << std::endl // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2229 << "using " << count << "/" << fMatrices.size() << " matrices:" << std::endl;
2230 for(auto* mat : fMatrices) {
2231 std::cout << mat << std::flush << " = " << mat->GetName() << std::endl;
2232 }
2233 }
2234
2235 fOldErrorLevel = gErrorIgnoreLevel;
2236 gErrorIgnoreLevel = kError;
2237 gPrintViaErrorHandler = true; // redirects all printf's to root's normal message system
2238
2239 if(!fLogFile.empty()) {
2240 fRedirect = new TRedirect(fLogFile.c_str());
2241 std::cout << "Starting redirect to " << fLogFile << std::endl;
2242 }
2243
2244 // check matrices (# of filled bins and bin labels) and resize some vectors for later use
2245 // use the first matrix to get a reference for everything
2246 fNofBins = 0;
2247 std::map<int, int> channelToIndex; // used to be a member, but only used here
2248 for(int bin = 1; bin <= fMatrices[0]->GetNbinsX(); ++bin) {
2249 if(FilledBin(fMatrices[0], bin) && std::find(fBadBins.begin(), fBadBins.end(), bin) == fBadBins.end()) {
2250 fActualChannelId.push_back(fNofBins); // at this point fNofBins is the index at which this projection will end up
2251 fActiveBins.push_back(bin);
2252 channelToIndex[bin] = fNofBins; // this map simply stores which bin ends up at which index
2253 fChannelLabel.push_back(fMatrices[0]->GetXaxis()->GetBinLabel(bin));
2254 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << bin << ". bin: fNofBins " << fNofBins << ", channelToIndex[" << bin << "] " << channelToIndex[bin] << ", fActualChannelId.back() " << fActualChannelId.back() << std::endl; }
2255 ++fNofBins;
2256 } else {
2257 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << "skipping bin " << bin << std::endl; }
2258 }
2259 }
2260 // now loop over all other matrices and check vs. the first one
2261 for(size_t i = 1; i < fMatrices.size(); ++i) {
2262 int tmpBins = 0;
2263 for(int bin = 1; bin <= fMatrices[i]->GetNbinsX(); ++bin) {
2264 if(FilledBin(fMatrices[i], bin) && std::find(fBadBins.begin(), fBadBins.end(), bin) == fBadBins.end()) { // good bin in the current matrix
2265 // current index is tmpBins, so we check if the bin agrees with what we have
2266 if(channelToIndex.find(bin) == channelToIndex.end() || tmpBins != channelToIndex[bin]) { // found a full bin, but the corresponding bin in the first matrix isn't full or a different bin
2267 std::ostringstream error;
2268 error << "Mismatch in " << i << ". matrix (" << fMatrices[i]->GetName() << "), bin " << bin << " is " << tmpBins << ". filled bin, but should be " << (channelToIndex.find(bin) == channelToIndex.end() ? "not filled" : Form("%d", channelToIndex[bin])) << std::endl;
2269 throw std::invalid_argument(error.str());
2270 }
2271 if(strcmp(fMatrices[0]->GetXaxis()->GetBinLabel(bin), fMatrices[i]->GetXaxis()->GetBinLabel(bin)) != 0) { // bin is full and matches the bin in the first matrix so we check the labels
2272 std::ostringstream error;
2273 error << i << ". matrix, " << bin << ". bin: label (" << fMatrices[i]->GetXaxis()->GetBinLabel(bin) << ") doesn't match bin label of the first matrix (" << fMatrices[0]->GetXaxis()->GetBinLabel(bin) << ")" << std::endl;
2274 throw std::invalid_argument(error.str());
2275 }
2276 ++tmpBins;
2277 } else {
2278 if(channelToIndex.find(bin) != channelToIndex.end()) { //found an empty bin that was full in the first matrix
2279 std::ostringstream error;
2280 error << "Mismatch in " << i << ". matrix (" << fMatrices[i]->GetName() << "), bin " << bin << " is empty, but should be " << channelToIndex[bin] << ". filled bin" << std::endl;
2281 throw std::invalid_argument(error.str());
2282 }
2283 }
2284 }
2285 }
2286
2287 if(fVerboseLevel > EVerbosity::kQuiet) { std::cout << fMatrices.size() << " matrices with " << fMatrices[0]->GetNbinsX() << " x-bins, fNofBins " << fNofBins << ", fActualChannelId.size() " << fActualChannelId.size() << std::endl; }
2288
2289 fOutput = new TFile("TSourceCalibration.root", "recreate");
2290 if(!fOutput->IsOpen()) {
2291 throw std::runtime_error("Unable to open output file \"TSourceCalibration.root\"!");
2292 }
2293
2294 // build the first screen
2297
2298 // Set a name to the main frame
2299 SetWindowName("SourceCalibration");
2300
2301 // Map all subwindows of main frame
2302 MapSubwindows();
2303
2304 // Initialize the layout algorithm
2305 Resize(GetDefaultSize());
2306
2307 // Map main frame
2308 MapWindow();
2309}
2310
2312{
2313 fOutput->Close();
2314 delete fOutput;
2315
2316 gErrorIgnoreLevel = fOldErrorLevel;
2317
2318 delete fRedirect;
2319}
2320
2322{
2323 /// Build initial interface with histogram <-> source assignment
2324
2325 auto* layoutManager = new TGTableLayout(this, fMatrices.size() + 1, 2, true, 5);
2326 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << "created table layout manager with 2 columns, " << fMatrices.size() + 1 << " rows" << std::endl; }
2327 SetLayoutManager(layoutManager);
2328
2329 // find the path for the .sou files
2330 std::string path = TNucleus::SourceDirectory();
2331 // The matrices and source selection boxes
2332 size_t i = 0;
2333 for(i = 0; i < fMatrices.size(); ++i) {
2334 fMatrixNames.push_back(new TGLabel(this, fMatrices[i]->GetName()));
2335 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "Text height " << fMatrixNames.back()->GetFont()->TextHeight() << std::endl; }
2336 fSource.push_back(nullptr);
2337 fSourceEnergy.emplace_back();
2338 fSourceBox.push_back(new TGComboBox(this, "Select source", kSourceBox + fSourceBox.size()));
2339
2340 int index = 0;
2341 // convert name to lower-case
2342 std::string name = fMatrices[i]->GetName();
2343 std::transform(name.begin(), name.end(), name.begin(), ::tolower);
2344 fMatrices[i]->SetName(name.c_str());
2345#if __cplusplus >= 201703L
2346 for(const auto& file : std::filesystem::directory_iterator(path)) {
2347 if(file.is_regular_file() && file.path().extension().compare(".sou") == 0) {
2348 fSourceBox.back()->AddEntry(file.path().stem().c_str(), index);
2349 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << i << "/" << index << ": Comparing matrix name " << fMatrices[i]->GetName() << " to file name " << file.path().stem().c_str() << std::endl; }
2350 if(std::strstr(fMatrices[i]->GetName(), file.path().stem().c_str()) != nullptr) {
2351 fSourceBox.back()->Select(index);
2352 SetSource(kSourceBox + fSourceBox.size() - 1, index);
2353 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << i << ": selected source " << index << std::endl; }
2355 std::cout << "matrix name " << fMatrices[i]->GetName() << " not matching " << file.path().stem().c_str() << std::endl;
2356 }
2357 ++index;
2358 }
2359 }
2360 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << i << ": created list with " << index << " sources" << std::endl; }
2361#else
2362 fSourceBox.back()->AddEntry("22Na", index);
2363 if(std::strstr(fMatrices[i]->GetName(), "22na") != nullptr) {
2364 fSourceBox.back()->Select(index);
2365 SetSource(kSourceBox + fSourceBox.size() - 1, index);
2366 }
2367 ++index;
2368 fSourceBox.back()->AddEntry("56Co", index);
2369 if(std::strstr(fMatrices[i]->GetName(), "56co") != nullptr) {
2370 fSourceBox.back()->Select(index);
2371 SetSource(kSourceBox + fSourceBox.size() - 1, index);
2372 }
2373 ++index;
2374 fSourceBox.back()->AddEntry("60Co", index);
2375 if(std::strstr(fMatrices[i]->GetName(), "60co") != nullptr) {
2376 fSourceBox.back()->Select(index);
2377 SetSource(kSourceBox + fSourceBox.size() - 1, index);
2378 }
2379 ++index;
2380 fSourceBox.back()->AddEntry("133Ba", index);
2381 if(std::strstr(fMatrices[i]->GetName(), "133ba") != nullptr) {
2382 fSourceBox.back()->Select(index);
2383 SetSource(kSourceBox + fSourceBox.size() - 1, index);
2384 }
2385 ++index;
2386 fSourceBox.back()->AddEntry("152Eu", index);
2387 if(std::strstr(fMatrices[i]->GetName(), "152eu") != nullptr) {
2388 fSourceBox.back()->Select(index);
2389 SetSource(kSourceBox + fSourceBox.size() - 1, index);
2390 }
2391 ++index;
2392 fSourceBox.back()->AddEntry("241Am", index);
2393 if(std::strstr(fMatrices[i]->GetName(), "241am") != nullptr) {
2394 fSourceBox.back()->Select(index);
2395 SetSource(kSourceBox + fSourceBox.size() - 1, index);
2396 }
2397#endif
2398
2399 fSourceBox.back()->SetMinHeight(fPanelHeight / 2);
2400
2401 fSourceBox.back()->Resize(fSourceboxWidth, fLineHeight);
2402 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "Attaching " << i << ". label to 0, 1, " << i << ", " << i + 1 << ", and box to 1, 2, " << i << ", " << i + 1 << std::endl; }
2403 AddFrame(fMatrixNames.back(), new TGTableLayoutHints(0, 1, i, i + 1, kLHintsRight | kLHintsCenterY));
2404 AddFrame(fSourceBox.back(), new TGTableLayoutHints(1, 2, i, i + 1, kLHintsLeft | kLHintsCenterY));
2405 }
2406
2407 // The buttons
2408 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << "Attaching start button to 0, 2, " << i << ", " << i + 1 << std::endl; }
2409 fStartButton = new TGTextButton(this, "Accept && Continue", kStartButton);
2410 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << "start button " << fStartButton << std::endl; }
2411 AddFrame(fStartButton, new TGTableLayoutHints(0, 2, i, i + 1, kLHintsCenterX | kLHintsCenterY));
2412 Layout();
2413}
2414
2416{
2417 /// Create connections for the histogram <-> source assignment interface
2418
2419 // Connect the selection of the source
2420 for(auto* box : fSourceBox) {
2421 box->Connect("Selected(Int_t, Int_t)", "TSourceCalibration", this, "SetSource(Int_t, Int_t)");
2422 }
2423
2424 //Connect the clicking of buttons
2425 fStartButton->Connect("Clicked()", "TSourceCalibration", this, "Start()");
2426}
2427
2429{
2430 /// Disconnect all signals from histogram <-> source assignment interface
2431 for(auto* box : fSourceBox) {
2432 box->Disconnect("Selected(Int_t, Int_t)", this, "SetSource(Int_t, Int_t)");
2433 }
2434
2435 fStartButton->Disconnect("Clicked()", this, "Start()");
2436}
2437
2439{
2440 HideFrame(element);
2441 RemoveFrame(element);
2442 delete element;
2443 element = nullptr;
2444}
2445
2447{
2448 for(auto* name : fMatrixNames) {
2449 DeleteElement(name);
2450 }
2451 fMatrixNames.clear();
2452 for(auto* box : fSourceBox) {
2453 DeleteElement(box);
2454 }
2455 fSourceBox.clear();
2457 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "Deleted start button " << fStartButton << std::endl; }
2458}
2459
2460void TSourceCalibration::SetSource(Int_t windowId, Int_t entryId)
2461{
2462 int index = windowId - kSourceBox;
2463 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << __PRETTY_FUNCTION__ << ": windowId " << windowId << ", entryId " << entryId << " => " << index << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2464 delete fSource[index];
2465 auto* nucleus = new TNucleus(fSourceBox[index]->GetListBox()->GetEntry(entryId)->GetTitle());
2466 TIter iter(nucleus->GetTransitionList());
2467 fSourceEnergy[index].clear();
2468 while(auto* transition = static_cast<TTransition*>(iter.Next())) {
2469 fSourceEnergy[index].emplace_back(transition->GetEnergy(), transition->GetEnergyUncertainty(), transition->GetIntensity(), transition->GetIntensityUncertainty());
2470 }
2471 fSource[index] = nucleus;
2472}
2473
2475{
2476 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << __PRETTY_FUNCTION__ << ": fEmitter " << fEmitter << ", fStartButton " << fStartButton << ", fAcceptAllButton " << fAcceptAllButton << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2477 if(fEmitter == fStartButton) {
2478 SecondWindow();
2479 } else if(fEmitter == fAcceptAllButton) {
2480 for(auto& channelTab : fChannelTab) {
2481 if(channelTab == nullptr) { continue; }
2482 channelTab->Write(fOutput);
2483 }
2485 std::cout << "Closing window" << std::endl;
2486 CloseWindow();
2487 std::cout << "all done" << std::endl;
2488 exit(0);
2489 }
2490}
2491
2493{
2494 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << __PRETTY_FUNCTION__ << ": fEmitter " << fEmitter << ", fStartButton " << fStartButton << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2495 if(fEmitter == nullptr) { // we only want to do this once at the beginning (after fEmitter was initialized to nullptr)
2497 TTimer::SingleShot(fWaitMs, "TSourceCalibration", this, "HandleTimer()");
2498 }
2499}
2500
2502{
2503 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << __PRETTY_FUNCTION__ << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2504 // check that all sources have been set
2505 for(size_t i = 0; i < fSource.size(); ++i) {
2506 if(fSource[i] == nullptr) {
2507 std::cerr << "Source " << i << " not set!" << std::endl;
2508 return;
2509 }
2510 }
2511 // now check that we don't have the same source twice (which wouldn't make sense)
2512 for(size_t i = 0; i < fSource.size(); ++i) {
2513 for(size_t j = i + 1; j < fSource.size(); ++j) {
2514 if(*(fSource[i]) == *(fSource[j])) {
2515 std::cerr << "Duplicate sources: " << i << " - " << fSource[i]->GetName() << " and " << j << " - " << fSource[j]->GetName() << std::endl;
2516 return;
2517 }
2518 }
2519 }
2520
2521 // disconnect signals of first screen and remove all elements
2523 RemoveAll();
2524 DeleteFirst();
2525
2526 SetLayoutManager(new TGVerticalLayout(this));
2527
2528 // create intermediate progress bar
2529 fProgressBar = new TGHProgressBar(this, TGProgressBar::kFancy, fPanelWidth);
2530 fProgressBar->SetBarColor("lightblue");
2531 fProgressBar->SetRange(0., static_cast<Float_t>(fMatrices.size() * fActiveBins.size()));
2532 fProgressBar->ShowPosition(true, false, "Creating source tabs: %.0f");
2533 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "Set range of progress bar to 0. - " << fProgressBar->GetMax() << " = " << fMatrices.size() * fActiveBins.size() << " = " << fMatrices.size() << "*" << fActiveBins.size() << std::endl; }
2534 AddFrame(fProgressBar, new TGLayoutHints(kLHintsCenterX | kLHintsCenterY | kLHintsExpandX | kLHintsExpandY, 0, 0, 0, 0));
2535
2536 // Map all subwindows of main frame
2537 MapSubwindows();
2538
2539 // Initialize the layout algorithm
2540 Resize(TGDimension(fPanelWidth, fLineHeight));
2541
2542 // Map main frame
2543 MapWindow();
2544
2545 // create second screen and its connections
2548
2549 // remove progress bar
2551
2552 // Map all subwindows of main frame
2553 MapSubwindows();
2554
2555 // Initialize the layout algorithm
2556 Resize(TGDimension(2 * fPanelWidth, fPanelHeight + 2 * fStatusbarHeight));
2557
2558 // Map main frame
2559 MapWindow();
2560 // re-draw to try and fix issue with last source histogram being on the wrong tab
2561 for(auto* channelTab : fChannelTab) {
2562 if(channelTab != nullptr) {
2563 channelTab->Draw();
2564 }
2565 }
2566 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << DGREEN << __PRETTY_FUNCTION__ << " done" << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2567}
2568
2570{
2571 SetLayoutManager(new TGVerticalLayout(this));
2572 fTab = new TGTab(this, 2 * fPanelWidth, fPanelHeight + 2 * fStatusbarHeight);
2573
2574 fChannelTab.resize(fMatrices[0]->GetNbinsX());
2575 for(auto& bin : fActiveBins) {
2576 // create vector with projections of this channel for each source
2577 std::vector<GH1D*> projections(fMatrices.size());
2578 size_t index = 0;
2579 for(auto* matrix : fMatrices) {
2580 projections[index] = new GH1D(matrix->ProjectionY(Form("%s_%s", fSource[index]->GetName(), matrix->GetXaxis()->GetBinLabel(bin)), bin, bin));
2581 ++index;
2582 }
2583 auto* frame = fTab->AddTab(fMatrices[0]->GetXaxis()->GetBinLabel(bin));
2584 // create a new tab in a new thread
2585 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << DGREEN << "Creating channel for bin " << bin << std::endl; }
2586 fChannelTab[bin - 1] = new TChannelTab(this, fSource, projections, fMatrices[0]->GetXaxis()->GetBinLabel(bin), frame, fSourceEnergy, fProgressBar);
2587 //// right now the lambda uses copies of all variables, not sure if that is neccesary, or if we only need that for projections?
2588 //fFutures.emplace(std::make_pair(bin - 1, std::async(std::launch::async, [=] {
2589 // std::cout << "creating channel tab with these arguments: " << this << ", fSource vector of size " << fSource.size() << ", projection vector of size " << projections.size() << ", " << fMatrices[0]->GetXaxis()->GetBinLabel(bin) << ", " << frame << ", " << fDefaultSigma << ", " << fThreshold << ", " << fDefaultDegree << ", " << fSourceEnergy.size() << " source energy vectors, " << fProgressBar << "."<< std::endl;
2590 // auto* newTab = new TChannelTab(this, fSource, projections, fMatrices[0]->GetXaxis()->GetBinLabel(bin), frame, fDefaultSigma, fThreshold, fDefaultDegree, fSourceEnergy, fProgressBar);
2591 // std::cout << "created new tab " << newTab << std::endl;
2592 // return newTab;
2593 // })));
2594 //// if we have too many threads running, get the result of the first one (this means waiting on that thread to finish) and assign the results
2595 //if(fFutures.size() >= fNumberOfThreads) {
2596 // if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << fFutures.size() << " threads running, waiting for thread #" << fFutures.front().first << " to finish before starting next, max. # threads " << fNumberOfThreads << std::endl; }
2597 // auto tmp = fFutures.front().second.get();
2598 // fChannelTab[fFutures.front().first] = tmp;
2599 // fFutures.pop();
2600 //}
2601 }
2602 // wait for all remaining threads to finish
2603 //while(!fFutures.empty()) {
2604 // if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << fFutures.size() << " threads left, waiting for thread #" << fFutures.front().first << " to finish before starting next, max. # threads " << fNumberOfThreads << std::endl; }
2605 // auto tmp = fFutures.front().second.get();
2606 // fChannelTab[fFutures.front().first] = tmp;
2607 // fFutures.pop();
2608 //}
2609
2610 //if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << fMatrices.size() << " matrices with " << fMatrices[0]->GetNbinsX() << " x-bins, fNofBins " << fNofBins << ", fActualChannelId.size() " << fActualChannelId.size() << ", fChannelTab[0]->SourceTab()->GetNumberOfTabs() " << fChannelTab[0]->SourceTab()->GetNumberOfTabs() << std::endl; }
2611
2612 //for(int i = 0; i < fChannelTab[0]->SourceTab()->GetNumberOfTabs(); ++i) {
2613 //if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << i << ": " << fChannelTab[0]->SourceTab()->GetTabTab(i)->GetText()->GetString() << std::endl; }
2614 //}
2615
2616 AddFrame(fTab, new TGLayoutHints(kLHintsTop | kLHintsExpandX | kLHintsExpandY, 0, 0, 0, 0));
2617
2618 // bottom frame with navigation button group, text entries, etc.
2619 fBottomFrame = new TGHorizontalFrame(this, 2 * fPanelWidth, TSourceCalibration::StatusbarHeight());
2620
2621 fLeftFrame = new TGVerticalFrame(fBottomFrame, fPanelWidth, fParameterHeight);
2622 fNavigationGroup = new TGHButtonGroup(fLeftFrame, "Navigation");
2623 fPreviousButton = new TGTextButton(fNavigationGroup, "&Previous");
2624 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << DGREEN << "prev button " << fPreviousButton << std::endl; }
2625 fPreviousButton->SetEnabled(false);
2626 fDiscardButton = new TGTextButton(fNavigationGroup, "&Discard");
2627 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "discard button " << fDiscardButton << std::endl; }
2628 fAcceptButton = new TGTextButton(fNavigationGroup, "&Accept");
2629 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "accept button " << fAcceptButton << std::endl; }
2630 fAcceptAllButton = new TGTextButton(fNavigationGroup, "Accept All && Finish");
2631 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "accept all button " << fAcceptAllButton << std::endl; }
2632 fWriteButton = new TGTextButton(fNavigationGroup, "&Write");
2633 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "write button " << fWriteButton << std::endl; }
2634 fNextButton = new TGTextButton(fNavigationGroup, "&Next");
2635 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "next button " << fNextButton << std::endl; }
2636
2637 fLeftFrame->AddFrame(fNavigationGroup, new TGLayoutHints(kLHintsTop | kLHintsExpandX, 2, 2, 2, 2));
2638
2639 fFittingGroup = new TGHButtonGroup(fLeftFrame, "Fitting/Calibrating");
2640 fFindSourcePeaksButton = new TGTextButton(fFittingGroup, "&Find Source Peaks");
2641 fFindChannelPeaksButton = new TGTextButton(fFittingGroup, "Find Channel Peaks");
2642 fFindAllPeaksButton = new TGTextButton(fFittingGroup, "Find All Peaks");
2643 fCalibrateButton = new TGTextButton(fFittingGroup, "&Calibrate");
2644 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "cal button " << fCalibrateButton << std::endl; }
2645 fRecursiveRemoveButton = new TGTextButton(fFittingGroup, "Recursive Remove");
2646 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "remove button " << fRecursiveRemoveButton << std::endl; }
2647
2648 fLeftFrame->AddFrame(fFittingGroup, new TGLayoutHints(kLHintsBottom | kLHintsExpandX, 2, 2, 2, 2));
2649
2650 fBottomFrame->AddFrame(fLeftFrame, new TGLayoutHints(kLHintsLeft, 2, 2, 2, 2));
2651
2652 fRightFrame = new TGVerticalFrame(fBottomFrame, fPanelWidth, fParameterHeight);
2653 fParameterFrame = new TGGroupFrame(fRightFrame, "Parameters", kHorizontalFrame);
2654 fParameterFrame->SetLayoutManager(new TGMatrixLayout(fParameterFrame, 0, 4, 2));
2655 fSigmaLabel = new TGLabel(fParameterFrame, "Sigma (in channels)");
2656 fSigmaEntry = new TGNumberEntry(fParameterFrame, fDefaultSigma, fDigitWidth, kSigmaEntry, TGNumberFormat::EStyle::kNESRealTwo, TGNumberFormat::EAttribute::kNEAPositive);
2657 fFitRangeLabel = new TGLabel(fParameterFrame, "Fit range (in channels)");
2658 fFitRangeEntry = new TGNumberEntry(fParameterFrame, fFitRange, fDigitWidth, kFitRangeEntry, TGNumberFormat::EStyle::kNESInteger, TGNumberFormat::EAttribute::kNEAPositive);
2659 fDegreeLabel = new TGLabel(fParameterFrame, "Degree of polynomial");
2660 fDegreeEntry = new TGNumberEntry(fParameterFrame, fDefaultDegree, 2, kDegreeEntry, TGNumberFormat::EStyle::kNESInteger, TGNumberFormat::EAttribute::kNEAPositive);
2661 fMaxResidualLabel = new TGLabel(fParameterFrame, "Max. residual (in keV)");
2662 fMaxResidualEntry = new TGNumberEntry(fParameterFrame, fDefaultMaxResidual, fDigitWidth, kMaxResidualEntry, TGNumberFormat::EStyle::kNESRealTwo, TGNumberFormat::EAttribute::kNEAPositive);
2663
2664 fParameterFrame->AddFrame(fSigmaLabel);
2665 fParameterFrame->AddFrame(fSigmaEntry);
2668 fParameterFrame->AddFrame(fDegreeLabel);
2669 fParameterFrame->AddFrame(fDegreeEntry);
2672
2673 fRightFrame->AddFrame(fParameterFrame, new TGLayoutHints(kLHintsTop | kLHintsRight | kLHintsExpandX, 2, 2, 2, 2));
2674
2675 fSettingsFrame = new TGGroupFrame(fRightFrame, "Settings", kHorizontalFrame);
2676 fSettingsFrame->SetLayoutManager(new TGMatrixLayout(fSettingsFrame, 0, 2, 2));
2677 fApplyToAll = new TGCheckButton(fSettingsFrame, "Apply to all channels", kApplyToAll);
2678 fApplyToAll->SetState(kButtonUp);
2679 fSettingsFrame->AddFrame(fApplyToAll);
2680 fWriteNonlinearities = new TGCheckButton(fSettingsFrame, "Write nonlinearities", kWriteNonlinearities);
2681 fWriteNonlinearities->SetState(kButtonUp);
2683
2684 fRightFrame->AddFrame(fSettingsFrame, new TGLayoutHints(kLHintsBottom | kLHintsRight | kLHintsExpandX, 2, 2, 2, 2));
2685
2686 fBottomFrame->AddFrame(fRightFrame, new TGLayoutHints(kLHintsBottom | kLHintsRight, 2, 2, 2, 2));
2687
2688 AddFrame(fBottomFrame, new TGLayoutHints(kLHintsBottom | kLHintsExpandX, 2, 2, 2, 2));
2689
2690 SelectedTab(0);
2691
2692 fProgressBar->Reset();
2693 fProgressBar->ShowPosition(true, false, "Finding peaks for tab #%.0f");
2694 for(auto* channelTab : fChannelTab) {
2695 if(channelTab != nullptr) {
2697 for(int i = 0; i < 20; ++i) { std::cout << std::endl; }
2698 for(int i = 0; i < 5; ++i) {
2699 std::cout << "================================================================================" << std::endl;
2700 }
2701 for(int i = 0; i < 20; ++i) { std::cout << std::endl; }
2702 }
2703 channelTab->Initialize(true);
2704 channelTab->FindAllCalibratedPeaks();
2705 channelTab->Calibrate();
2706 channelTab->Draw();
2708 channelTab->PrintCanvases();
2709 }
2710 }
2711 }
2712
2713 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << DGREEN << "Second interface done" << std::endl; }
2714}
2715
2717{
2718 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << DGREEN << __PRETTY_FUNCTION__ << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2719 fNavigationGroup->Connect("Clicked(Int_t)", "TSourceCalibration", this, "Navigate(Int_t)");
2720 fFittingGroup->Connect("Clicked(Int_t)", "TSourceCalibration", this, "Fitting(Int_t)");
2721 fTab->Connect("Selected(Int_t)", "TSourceCalibration", this, "SelectedTab(Int_t)");
2722 // we don't need to connect the sigma, threshold, and degree number entries, those are automatically read when we start the calibration
2723 for(auto* channelTab : fChannelTab) {
2724 if(channelTab != nullptr) {
2725 channelTab->MakeConnections();
2726 }
2727 }
2728}
2729
2731{
2732 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << DGREEN << __PRETTY_FUNCTION__ << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2733 fNavigationGroup->Disconnect("Clicked(Int_t)", this, "Navigate(Int_t)");
2734 fFittingGroup->Disconnect("Clicked(Int_t)", this, "Fitting(Int_t)");
2735 fTab->Disconnect("Selected(Int_t)", this, "SelectedTab(Int_t)");
2736 for(auto* channelTab : fChannelTab) {
2737 if(channelTab != nullptr) {
2738 channelTab->Disconnect();
2739 }
2740 }
2741}
2742
2744{
2745 // we get the current source tab id and use it to get the channel tab from the right source tab
2746 // since the current id only refers to the position within the open tabs we need to keep track of the actual id it relates to
2747 // for this we created a vector with all ids at the beginning and now remove the position correspoding to the current id when we remove the tab
2748 int currentChannelId = fTab->GetCurrent();
2749 int actualChannelId = fActualChannelId[currentChannelId];
2750 int nofTabs = fTab->GetNumberOfTabs();
2751 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << DGREEN << __PRETTY_FUNCTION__ << ": id " << id << ", channel tab id " << currentChannelId << ", actual channel tab id " << actualChannelId << ", # of tabs " << nofTabs << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2752 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << "Before: active source tab of tab " << fChannelTab[fActiveBins[fTab->GetCurrent()] - 1]->Name() << ", #" << fTab->GetCurrent() << " = bin " << fActiveBins[fTab->GetCurrent()] << ": " << fChannelTab[fActiveBins[fTab->GetCurrent()] - 1]->ActiveSourceTab() << std::endl; }
2753 switch(id) {
2754 case ENavigate::kPrevious: // previous
2755 fTab->SetTab(currentChannelId - 1);
2756 SelectedTab(currentChannelId - 1);
2757 break;
2758 case ENavigate::kDiscard: // discard
2759 // select the next (or if we are on the last tab, the previous) tab
2760 if(currentChannelId < nofTabs - 1) {
2761 fTab->SetTab(currentChannelId + 1);
2762 } else {
2763 fTab->SetTab(currentChannelId - 1);
2764 }
2765 // remove the original active tab
2766 fTab->RemoveTab(currentChannelId);
2767 break;
2768 case ENavigate::kAccept: // accept
2769 AcceptChannel(currentChannelId);
2770 break;
2771 case ENavigate::kAcceptAll: // accept all (no argument = -1 = all)
2772 AcceptChannel();
2773 return;
2774 case ENavigate::kWrite: // write
2776 break;
2777 case ENavigate::kNext: // next
2778 fTab->SetTab(currentChannelId + 1);
2779 SelectedTab(currentChannelId + 1);
2780 break;
2781 default:
2782 break;
2783 }
2784 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << DGREEN << "After: active source tab of tab " << fChannelTab[fActiveBins[fTab->GetCurrent()] - 1]->Name() << ", #" << fTab->GetCurrent() << " = bin " << fActiveBins[fTab->GetCurrent()] << ": " << fChannelTab[fActiveBins[fTab->GetCurrent()] - 1]->ActiveSourceTab() << std::endl; }
2785}
2786
2788{
2789 // we get the current source tab id and use it to get the channel tab from the right source tab
2790 // since the current id only refers to the position within the open tabs we need to keep track of the actual id it relates to
2791 // for this we created a vector with all ids at the beginning and now remove the position correspoding to the current id when we remove the tab
2792 int currentChannelId = fTab->GetCurrent();
2793 int actualChannelId = fActualChannelId[currentChannelId];
2794 int nofTabs = fTab->GetNumberOfTabs();
2795 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << DGREEN << __PRETTY_FUNCTION__ << ": id " << id << ", channel tab id " << currentChannelId << ", actual channel tab id " << actualChannelId << ", # of tabs " << nofTabs << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2796 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << "Before: active source tab of tab " << fChannelTab[fActiveBins[fTab->GetCurrent()] - 1]->Name() << ", #" << fTab->GetCurrent() << " = bin " << fActiveBins[fTab->GetCurrent()] << ": " << fChannelTab[fActiveBins[fTab->GetCurrent()] - 1]->ActiveSourceTab() << std::endl; }
2797 switch(id) {
2798 case EFitting::kFindSourcePeaks: // find peaks for selected source
2800 break;
2801 case EFitting::kFindChannelPeaks: // find peaks for selected channel
2803 break;
2804 case EFitting::kFindAllPeaks: // find all peaks
2805 FindAllPeaks();
2806 break;
2807 case EFitting::kCalibrate: // calibrate
2808 Calibrate();
2809 break;
2810 case EFitting::kRecursiveRemove: // recursively remove outliers
2812 break;
2813 default:
2814 break;
2815 }
2816 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << DGREEN << "After: active source tab of tab " << fChannelTab[fActiveBins[fTab->GetCurrent()] - 1]->Name() << ", #" << fTab->GetCurrent() << " = bin " << fActiveBins[fTab->GetCurrent()] << ": " << fChannelTab[fActiveBins[fTab->GetCurrent()] - 1]->ActiveSourceTab() << std::endl; }
2817}
2818
2820{
2821 /// Simple function that enables and disables the previous and next buttons depending on which tab was selected.
2822 /// Also sets gPad to the projection of the source selected in this tab.
2823 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << DGREEN << __PRETTY_FUNCTION__ << ": id " << id << ", nofTabs " << fTab->GetNumberOfTabs() << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2824 if(id < 0) { return; }
2825 if(id == 0) {
2826 fPreviousButton->SetEnabled(false);
2827 } else {
2828 fPreviousButton->SetEnabled(true);
2829 }
2830
2831 if(id >= fTab->GetNumberOfTabs() - 1) {
2832 fNextButton->SetEnabled(false);
2833 } else {
2834 fNextButton->SetEnabled(true);
2835 }
2836
2837 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << DGREEN << "active source tab of tab " << fChannelTab[fActiveBins[fTab->GetCurrent()] - 1]->Name() << ", #" << fTab->GetCurrent() << " = bin " << fActiveBins[fTab->GetCurrent()] << ": " << fChannelTab[fActiveBins[fTab->GetCurrent()] - 1]->ActiveSourceTab() << RESET_COLOR << std::endl; }
2838
2839 if(fChannelTab[fActiveBins[id] - 1] == nullptr) {
2840 std::cout << "Failed to get fChannelTab[" << fActiveBins[id] - 1 << " = fActiveBins[" << id << "]-1]" << RESET_COLOR << std::endl;
2841 return;
2842 }
2843 if(fChannelTab[fActiveBins[id] - 1]->SelectedSourceTab() == nullptr) {
2844 std::cout << "Failed to get fChannelTab[" << fActiveBins[id] - 1 << " = fActiveBins[" << id << "]-1]->SelectedSourceTab()" << RESET_COLOR << std::endl;
2845 return;
2846 }
2847 if(fChannelTab[fActiveBins[id] - 1]->SelectedSourceTab()->ProjectionCanvas() == nullptr) {
2848 std::cout << "Failed to get fChannelTab[" << fActiveBins[id] - 1 << " = fActiveBins[" << id << "]-1]->SelectedSourceTab()->ProjectionCanvas()" << RESET_COLOR << std::endl;
2849 return;
2850 }
2851 if(fChannelTab[fActiveBins[id] - 1]->SelectedSourceTab()->ProjectionCanvas()->GetCanvas() == nullptr) {
2852 std::cout << "Failed to get fChannelTab[" << fActiveBins[id] - 1 << " = fActiveBins[" << id << "]-1]->SelectedSourceTab()->ProjectionCanvas()->GetCanvas()" << RESET_COLOR << std::endl;
2853 return;
2854 }
2855 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << DGREEN << "Changing gpad from \"" << gPad->GetName() << "\" to \""; }
2856 gPad = fChannelTab[fActiveBins[id] - 1]->SelectedSourceTab()->ProjectionCanvas()->GetCanvas();
2857 fChannelTab[fActiveBins[id] - 1]->SelectedSourceTab()->ProjectionCanvas()->GetCanvas()->cd();
2858 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << gPad->GetName() << "\"" << RESET_COLOR << std::endl; }
2859}
2860
2861void TSourceCalibration::AcceptChannel(const int& channelId)
2862{
2863 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << DGREEN << __PRETTY_FUNCTION__ << ": channelId " << channelId << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2864
2865 // select the next (or if we are on the last tab, the previous) tab
2866 int nofTabs = fTab->GetNumberOfTabs();
2867 int minChannel = 0;
2868 int maxChannel = nofTabs - 1;
2869 if(channelId >= 0) {
2870 minChannel = channelId;
2871 maxChannel = channelId;
2872 }
2873
2874 // we need to loop backward, because removing the first channel would make the second one the new first and so on
2875 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << __PRETTY_FUNCTION__ << ": accepting channels " << maxChannel << " to " << minChannel << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2876 for(int currentChannelId = maxChannel; currentChannelId >= minChannel; --currentChannelId) {
2877 int actualChannelId = fActualChannelId[currentChannelId];
2878 if(fVerboseLevel > EVerbosity::kLoops) { std::cout << __PRETTY_FUNCTION__ << ": currentChannelId " << currentChannelId << ", actualChannelId " << actualChannelId << ", fChannelTab.size() " << fChannelTab.size() << ", fActualChannelId.size() " << fActualChannelId.size() << ", nofTabs " << nofTabs << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2879
2880 if(minChannel == maxChannel) { // we don't need to select the tab if we close all
2881 if(currentChannelId < nofTabs) {
2882 fTab->SetTab(currentChannelId + 1);
2883 if(fVerboseLevel > EVerbosity::kLoops) { std::cout << DGREEN << __PRETTY_FUNCTION__ << ": set tab to " << currentChannelId + 1 << " = " << fTab->GetCurrent() << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2884 } else {
2885 fTab->SetTab(currentChannelId - 1);
2886 if(fVerboseLevel > EVerbosity::kLoops) { std::cout << DGREEN << __PRETTY_FUNCTION__ << ": set tab to " << currentChannelId - 1 << " = " << fTab->GetCurrent() << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2887 }
2888 }
2889
2890 // remove the original active tab
2891 fTab->RemoveTab(currentChannelId); // this sets the selected tab to zero if the deleted tab is the currently selected one, which is why we do it now (after we've selected a new tab)
2892 fActualChannelId.erase(fActualChannelId.begin() + currentChannelId);
2893 fActiveBins.erase(fActiveBins.begin() + currentChannelId);
2894 }
2895 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << DGREEN << __PRETTY_FUNCTION__ << ": # of channel tabs " << fActualChannelId.size() << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2896 // check if this changes which buttons are enabled or not (not using the variables from the beginning of this block because they might have changed by now!)
2897 SelectedTab(fTab->GetCurrent());
2898 // if this was also the last source vector we initiate the last screen
2899 if(fActualChannelId.empty() || fActiveBins.empty()) {
2900 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "last channel tab done - writing files now" << std::endl; }
2902 TTimer::SingleShot(fWaitMs, "TSourceCalibration", this, "HandleTimer()");
2903 }
2904 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << RESET_COLOR << std::flush; }
2905}
2906
2908{
2909 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << DGREEN << __PRETTY_FUNCTION__ << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2910 if(fApplyToAll->IsDown()) {
2911 auto* window = CreateProgressbar("Finding current source peaks in channel #%.0f");
2912
2913 for(auto& bin : fActiveBins) {
2914 if(fChannelTab[bin - 1] != nullptr) {
2915 fChannelTab[bin - 1]->FindCalibratedPeaks();
2916 fChannelTab[bin - 1]->Calibrate();
2917 }
2918 fProgressBar->Increment(1);
2919 }
2920 window->CloseWindow();
2921 delete fProgressBar;
2922 } else if(fChannelTab[fActiveBins[fTab->GetCurrent()] - 1] != nullptr) {
2923 fChannelTab[fActiveBins[fTab->GetCurrent()] - 1]->FindCalibratedPeaks();
2924 fChannelTab[fActiveBins[fTab->GetCurrent()] - 1]->Calibrate();
2925 } else {
2926 std::cout << __PRETTY_FUNCTION__ << ": fChannelTab[" << fActiveBins[fTab->GetCurrent()] << " = fActiveBins[" << fTab->GetCurrent() << "]] is a nullptr (fApplyToAll is " << (fApplyToAll->IsDown() ? "down" : "up") << ")!" << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2927 }
2928 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << RESET_COLOR << std::flush; }
2929}
2930
2932{
2933 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << DGREEN << __PRETTY_FUNCTION__ << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2934 if(fChannelTab[fActiveBins[fTab->GetCurrent()] - 1] != nullptr) {
2935 fChannelTab[fActiveBins[fTab->GetCurrent()] - 1]->FindAllCalibratedPeaks();
2936 fChannelTab[fActiveBins[fTab->GetCurrent()] - 1]->Calibrate();
2937 } else {
2938 std::cout << __PRETTY_FUNCTION__ << ": fChannelTab[" << fActiveBins[fTab->GetCurrent()] << " = fActiveBins[" << fTab->GetCurrent() << "]] is a nullptr!" << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2939 }
2940 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << RESET_COLOR << std::flush; }
2941}
2942
2944{
2945 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << DGREEN << __PRETTY_FUNCTION__ << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2946 auto* window = CreateProgressbar("Finding all source peaks for channel #%.0f");
2947
2948 for(auto& bin : fActiveBins) {
2949 if(fChannelTab[bin - 1] != nullptr) {
2950 fChannelTab[bin - 1]->FindAllCalibratedPeaks();
2951 fChannelTab[bin - 1]->Calibrate();
2952 } else {
2953 std::cout << __PRETTY_FUNCTION__ << ": fChannelTab[" << bin << "] is a nullptr!" << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2954 }
2955 fProgressBar->Increment(1);
2956 }
2957 window->CloseWindow();
2958 delete fProgressBar;
2959 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << RESET_COLOR << std::flush; }
2960}
2961
2963{
2964 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << DGREEN << __PRETTY_FUNCTION__ << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2965 if(fApplyToAll->IsDown()) {
2966 auto* window = CreateProgressbar("Calibrating channel #%.0f");
2967
2968 for(auto& bin : fActiveBins) {
2969 if(fChannelTab[bin - 1] != nullptr) {
2970 fChannelTab[bin - 1]->Calibrate();
2971 }
2972 fProgressBar->Increment(1);
2973 }
2974 window->CloseWindow();
2975 delete fProgressBar;
2976 } else if(fChannelTab[fActiveBins[fTab->GetCurrent()] - 1] != nullptr) {
2977 fChannelTab[fActiveBins[fTab->GetCurrent()] - 1]->Calibrate();
2978 } else {
2979 std::cout << __PRETTY_FUNCTION__ << ": fChannelTab[" << fActiveBins[fTab->GetCurrent()] << " = fActiveBins[" << fTab->GetCurrent() << "]] is a nullptr (fApplyToAll is " << (fApplyToAll->IsDown() ? "down" : "up") << ")!" << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2980 }
2981 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << RESET_COLOR << std::flush; }
2982}
2983
2985{
2986 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << DGREEN << __PRETTY_FUNCTION__ << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2987 if(fApplyToAll->IsDown()) {
2988 auto* window = CreateProgressbar("Recursively removing peaks for channel #%.0f");
2989
2990 for(auto& bin : fActiveBins) {
2991 if(fChannelTab[bin - 1] != nullptr) {
2992 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << "Starting recursive remove for bin " << bin << ", fChannelTab[bin - 1] " << fChannelTab[bin - 1] << std::endl; }
2993 fChannelTab[bin - 1]->RecursiveRemove(MaxResidual());
2994 }
2995 fProgressBar->Increment(1);
2996 }
2997 window->CloseWindow();
2998 delete fProgressBar;
2999 } else if(fChannelTab[fActiveBins[fTab->GetCurrent()] - 1] != nullptr) {
3000 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << "Starting recursive remove for bin " << fActiveBins[fTab->GetCurrent() - 1] << ", fChannelTab[fActiveBins[fTab->GetCurrent()] - 1] " << fChannelTab[fActiveBins[fTab->GetCurrent()] - 1] << std::endl; }
3001 fChannelTab[fActiveBins[fTab->GetCurrent()] - 1]->RecursiveRemove(MaxResidual());
3002 } else {
3003 std::cout << __PRETTY_FUNCTION__ << ": fChannelTab[" << fActiveBins[fTab->GetCurrent()] << " = fActiveBins[" << fTab->GetCurrent() << "]] is a nullptr (fApplyToAll is " << (fApplyToAll->IsDown() ? "down" : "up") << ")!" << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
3004 }
3005 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << RESET_COLOR << std::flush; }
3006}
3007
3009{
3010 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << DGREEN << __PRETTY_FUNCTION__ << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
3011 for(auto& channelTab : fChannelTab) {
3012 if(channelTab == nullptr) { continue; }
3013 channelTab->UpdateChannel();
3014 }
3015 std::ostringstream fileName;
3016 for(auto* source : fSource) {
3017 fileName << source->GetName();
3018 }
3019 fileName << ".cal";
3020 TChannel::WriteCalFile(fileName.str());
3021 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << DGREEN << __PRETTY_FUNCTION__ << ": wrote " << fileName.str() << " with " << TChannel::GetNumberOfChannels() << " channels" << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
3023}
3024
3025TGTransientFrame* TSourceCalibration::CreateProgressbar(const char* format)
3026{
3027 auto* window = new TGTransientFrame(gClient->GetRoot(), this, fPanelWidth, fStatusbarHeight);
3028 fProgressBar = new TGHProgressBar(window, TGProgressBar::kFancy, fPanelWidth);
3029 fProgressBar->SetBarColor("lightblue");
3030 fProgressBar->SetRange(0., static_cast<Float_t>(fActiveBins.size()));
3031 fProgressBar->ShowPosition(true, false, format);
3032 window->AddFrame(fProgressBar, new TGLayoutHints(kLHintsCenterX | kLHintsCenterY | kLHintsExpandX | kLHintsExpandY, 0, 0, 0, 0));
3033 window->CenterOnParent();
3034 window->MapSubwindows();
3035 window->MapWindow();
3036 window->Resize(window->GetDefaultSize());
3037
3038 return window;
3039}
3040
3042{
3043 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << DGREEN << __PRETTY_FUNCTION__ << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
3044
3045 Print();
3046 if(fBottomFrame != nullptr) { fBottomFrame->Print(); }
3047 if(fLeftFrame != nullptr) { fLeftFrame->Print(); }
3048 if(fRightFrame != nullptr) { fRightFrame->Print(); }
3049 if(fTab != nullptr) { fTab->Print(); }
3050 std::cout << "TSourceCalibration all " << fChannelTab.size() << " channel tabs:" << std::endl;
3051 for(auto* channelTab : fChannelTab) {
3052 if(channelTab != nullptr) { channelTab->PrintLayout(); }
3053 }
3054}
3055
3057{
3058 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << DGREEN << __PRETTY_FUNCTION__ << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
3059
3060 for(auto* channelTab : fChannelTab) {
3061 if(channelTab != nullptr) { channelTab->PrintCanvases(); }
3062 }
3063}
#define kArrowKeyRelease
Definition GCanvas.cxx:49
#define kArrowKeyPress
Definition GCanvas.cxx:48
#define DYELLOW
Definition Globals.h:16
#define DRED
Definition Globals.h:18
#define DGREEN
Definition Globals.h:17
#define DCYAN
Definition Globals.h:21
std::string hex(T val, int width=-1)
Definition Globals.h:129
EVerbosity
Definition Globals.h:143
#define RESET_COLOR
Definition Globals.h:5
TH2D * mat
Definition UserFillObj.h:12
Definition GH1D.h:19
TList * ListOfRegions()
Definition GH1D.h:70
void Print(Option_t *opt="") const override
Definition GH1D.cxx:95
void Clear(Option_t *opt="") override
Definition GH1D.cxx:89
void Draw(Option_t *opt="") override
Definition GH1D.cxx:108
void DrawResidual(Option_t *opt="", TLegend *legend=nullptr)
int GetN()
Returns GetN(), i.e. number of points of the total graph.
void SetMarkerColor(int index, Color_t color)
bool SetResidual(const bool &force=false)
TGraphErrors * TotalResidualGraph()
double GetMaximumY() const
Return maximum y-value.
void SetLineColor(int index, Color_t color)
double * GetY()
Returns an array of y-values of the total graph.
void Print(Option_t *opt="") const override
void Fit(TF1 *function, Option_t *opt="")
Fits the provided function to the total graph.
double GetMinimumY() const
Return minimum y-value.
double GetMinimumX() const
Return minimum x-value.
TGraphErrors * TotalGraph()
double GetMaximumX() const
Return maximum x-value.
int Add(TGraphErrors *, const std::string &label)
Add new graph to set, using the label when creating legends during plotting.
void DrawCalibration(Option_t *opt="", TLegend *legend=nullptr)
void Clear(Option_t *option="") override
Int_t RemovePoint(const Int_t &px, const Int_t &py)
double * GetX()
Returns an array of x-values of the total graph.
TF1 * FitFunction()
Gets the calibration from the total graph (might be nullptr!).
static void WriteCalFile(const std::string &outfilename="")
TGraph GetEnergyNonlinearity() const
Definition TChannel.h:195
static TChannel * GetChannel(unsigned int temp_address, bool warn=false)
Definition TChannel.cxx:476
void SetENGCoefficients(const std::vector< Float_t > &tmp, size_t range=0)
Definition TChannel.h:225
void DestroyEnergyNonlinearity()
Definition TChannel.cxx:610
void AddEnergyNonlinearityPoint(double x, double y)
Definition TChannel.h:214
static size_t GetNumberOfChannels()
Definition TChannel.h:63
TSourceCalibration * fSourceCalibration
the parent of this tab
std::vector< GH1D * > fProjections
vector of all projections
void PrintLayout() const
TGCompositeFrame * fChannelFrame
main frame of this tab
std::vector< TSourceTab * > fSources
tabs for all sources
TCalibrationGraphSet * fFwhm
combined fwhm from all sources
TGraphErrors * Data(int channelId) const
TPaveText * fChi2Label
TGTab * fSourceTab
tab for sources
void SelectCanvas(Event_t *event)
void Write(TFile *output)
TRootEmbeddedCanvas * fInitCanvas
TChannelTab(TSourceCalibration *sourceCal, std::vector< TNucleus * > nuclei, std::vector< GH1D * > projections, std::string name, TGCompositeFrame *frame, std::vector< std::vector< std::tuple< double, double, double, double > > > sourceEnergies, TGHProgressBar *progressBar)
void CalibrationStatus(Int_t event, Int_t px, Int_t py, TObject *selected)
void Initialize(const bool &force)
static void ZoomY()
void RemovePoint(Int_t oldGraph, Int_t oldPoint)
TGCompositeFrame * fCalibrationFrame
frame of tab with calibration
std::vector< std::vector< std::tuple< double, double, double, double > > > fSourceEnergies
vector with source energies and uncertainties
std::string Name() const
void CreateSourceTab(size_t source)
TGCompositeFrame * fFwhmFrame
frame of tab with fwhm
TCalibrationGraphSet * fInit
combined initial data from all sources
TRootEmbeddedCanvas * fFwhmCanvas
int fActiveSourceTab
id of the currently active source tab
std::string fName
name of this tab
TPaveText * fInfoLabel
static void ZoomX()
std::vector< TNucleus * > fNuclei
the source nucleus
void SelectedTab(Int_t id)
TGStatusBar * fChannelStatusBar
void RecursiveRemove(const double &maxResidual)
TGTab * fCanvasTab
tab for canvases (calibration with residuals, and FWHM)
TPaveText * fCalLabel
TGCompositeFrame * fInitFrame
frame of tab with initial calibration
void PrintCanvases() const
TGHProgressBar * fProgressBar
TRootEmbeddedCanvas * fCalibrationCanvas
TCalibrationGraphSet * fData
combined data from all sources
Double_t CentroidErr() const override
Definition TGauss.cxx:65
void Centroid(const Double_t &centroid) override
Definition TGauss.cxx:13
static std::string & SourceDirectory()
Definition TNucleus.cxx:19
void SetLineColor(Color_t color)
Definition TSinglePeak.h:91
Double_t AreaErr() const
Definition TSinglePeak.h:52
Double_t Area() const
Definition TSinglePeak.h:51
unsigned int fLineHeight
Height of text boxes and progress bar.
int fWaitMs
How many milliseconds we wait before we process the navigation input (to avoid double triggers?...
TGTextButton * fPreviousButton
TRedirect * fRedirect
redirect, created in constructor and destroyed in destructor if a log file name has been provided
static std::string LogFile()
static int fStatusbarHeight
Height of status bar (also extra height needed for tabs)
TGVerticalFrame * fRightFrame
TGTransientFrame * CreateProgressbar(const char *format)
static std::string fLogFile
name of log file, if empty no log file is written
static void FitRange(int val)
TSourceCalibration(double sigma, double threshold, int degree, int count...)
TGCheckButton * fWriteNonlinearities
TGNumberEntry * fMaxResidualEntry
int fDefaultDegree
The default degree of the polynomial used for calibrating, can be changed later.
static EVerbosity VerboseLevel()
TGGroupFrame * fParameterFrame
std::vector< TGLabel * > fMatrixNames
TGTextButton * fFindSourcePeaksButton
TGTextButton * fCalibrateButton
static double MinIntensity()
TGTextButton * fAcceptAllButton
void AcceptChannel(const int &channelId=-1)
TGNumberEntry * fSigmaEntry
double fDefaultSigma
The default sigma used for the peak finding algorithm, can be changed later.
static int fMaxIterations
Maximum iterations over combinations in Match and SmartMatch.
TGTextButton * fFindAllPeaksButton
TGHProgressBar * fProgressBar
TGHButtonGroup * fNavigationGroup
void SetSource(Int_t windowId, Int_t entryId)
std::vector< int > fActiveBins
TGTextButton * fDiscardButton
std::vector< TNucleus * > fSource
TGNumberEntry * fFitRangeEntry
static bool fFast
Do we use the fast peak finding method on startup or not.
std::vector< const char * > fChannelLabel
std::vector< TChannelTab * > fChannelTab
static bool fAcceptBadFits
Do we accept peaks where the fit was bad (position or area uncertainties too large or not numbers)
TGTextButton * fAcceptButton
static int fPanelWidth
Width of one panel.
static size_t fNumberOfThreads
Maximum number of threads to run while creating the channel tabs.
double fDefaultMaxResidual
The default maximum residual accepted when trying to iteratively find peaks.
static int fPanelHeight
Height of one panel.
TGTextButton * fRecursiveRemoveButton
std::vector< TH2 * > fMatrices
static bool fUseCalibratedPeaks
Do we use the initial calibration to find more peaks in the source spectra?
int fOldErrorLevel
Used to store old value of gErrorIgnoreLevel (set to kError for the scope of the class)
static int fDigitWidth
Number of digits used for parameter entries (if they are floating point)
TGHButtonGroup * fFittingGroup
TGTextButton * fWriteButton
std::vector< TGComboBox * > fSourceBox
std::vector< int > fActualChannelId
int fNofBins
Number of filled bins in first matrix.
TGTextButton * fFindChannelPeaksButton
TGNumberEntry * fDegreeEntry
static double fMinIntensity
Minimum intensity of source peaks to be considered when trying to find them in the spectrum.
void DeleteElement(TGFrame *element)
static int fParameterHeight
Height of the frame for the parameters.
std::vector< std::vector< std::tuple< double, double, double, double > > > fSourceEnergy
vector to hold source energy, energy uncertainty, intensity, and intensity uncertainty
static bool AcceptBadFits()
static std::vector< int > fBadBins
Bins of the 2D matrix to be ignored even if there are enough counts in them.
static int fSourceboxWidth
Width of the box to select which source each histogram is from.
TGTextButton * fStartButton
TGCheckButton * fApplyToAll
TGTextButton * fNextButton
TGVerticalFrame * fLeftFrame
static EVerbosity fVerboseLevel
Changes verbosity from kQuiet to kAll.
TGGroupFrame * fSettingsFrame
TGHorizontalFrame * fBottomFrame
static int fFitRange
range in sigma used for fitting peaks (peak position - range to peas position + range)
void InitialCalibration(const bool &force)
void Status(const char *status, int position)
void ReplacePeak(const size_t &index, const double &channel)
const char * fSourceName
TGraphErrors * fFwhm
bool Good(TGauss *peak)
TRootEmbeddedCanvas * fProjectionCanvas
std::vector< std::pair< double, double > > fRegions
void RemovePoint(Int_t oldPoint)
TChannelTab * fChannelTab
void ProjectionStatus(Event_t *event)
TSourceTab(TSourceCalibration *sourceCal, TChannelTab *channel, TGCompositeFrame *frame, GH1D *projection, const char *sourceName, std::vector< std::tuple< double, double, double, double > > sourceEnergy)
TGCompositeFrame * fSourceFrame
TGraphErrors * fData
std::vector< TPeakFitter * > fBadFits
all bad fits (centroid err > 10%, area err > 100%, or either of them not a number at all)
std::vector< std::tuple< double, double, double, double > > fSourceEnergy
gamma rays from the source, with their energies, uncertainties in the energies, intensities,...
std::vector< TPeakFitter * > fFits
all good fits
void Add(std::map< double, std::tuple< double, double, double, double > > map)
TSourceCalibration * fSourceCalibration
void PrintLayout() const
std::vector< TGauss * > fPeaks
all peaks that have been fitted and are good
void FindCalibratedPeaks(const TF1 *calibration)
TGStatusBar * fSourceStatusBar
void PrintCanvases() const
std::map< TGauss *, std::tuple< double, double, double, double > > SmartMatch(std::vector< TGauss * > peaks, std::vector< std::tuple< double, double, double, double > > sources, TSourceTab *sourceTab)
std::map< TGauss *, std::tuple< double, double, double, double > > Match(std::vector< TGauss * > peaks, std::vector< std::tuple< double, double, double, double > > sources, TSourceTab *sourceTab)
bool FilledBin(TH2 *matrix, const int &bin)
std::map< double, std::tuple< double, double, double, double > > RoughCal(std::vector< double > peaks, std::vector< std::tuple< double, double, double, double > > sources, TSourceTab *sourceTab)
double Polynomial(double *x, double *par)