GRSISort "v4.0.0.5"
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 // default sorting of tuples is by the first entry, which is the source energy here, so no need for a custom sort condition
468 std::sort(fSourceEnergy.begin(), fSourceEnergy.end());
470 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << DCYAN << "interface built, finding peaks ... " << RESET_COLOR << std::endl; }
471}
472
474{
475 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)
480
482 fData = rhs.fData;
483 fFwhm = rhs.fFwhm;
484 // ? not sure what to do here, clearing the vectors seems unnecessary since we never filled them
485 // but we also can't just simply copy them, we would need a deep copy, which might not make sense?
486 // does is make sense in general to have assignment constructors for this class?
487 fFits.clear();
488 fBadFits.clear();
489 fPeaks.clear();
490}
491
493{
494 for(auto* fit : fFits) {
495 delete fit;
496 }
497 fFits.clear();
498 for(auto* fit : fBadFits) {
499 delete fit;
500 }
501 fBadFits.clear();
502 for(auto* peak : fPeaks) {
503 delete peak;
504 }
505 fPeaks.clear();
506 delete fSourceFrame;
507 delete fProjectionCanvas;
508 delete fSourceStatusBar;
509 //delete fProjection;
510 //delete fData;
511 //delete fFwhm;
512}
513
515{
516 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)
517 //if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << "Creating unique lock on graphics mutex to build interface " << std::endl; }
518 //std::unique_lock<std::mutex> graphicsLock(fSourceCalibration->GraphicsMutex());
519 //if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << "Got unique lock on graphics mutex!" << std::endl; }
520 // frame with canvas and status bar
521 fProjectionCanvas = new TRootEmbeddedCanvas(Form("ProjectionCanvas%s", fProjection->GetName()), fSourceFrame, TSourceCalibration::PanelWidth(), TSourceCalibration::PanelHeight());
522
523 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << "Created projection canvas!" << std::endl; }
524
525 fSourceFrame->AddFrame(fProjectionCanvas, new TGLayoutHints(kLHintsTop | kLHintsExpandX | kLHintsExpandY, 2, 2, 2, 2));
526
527 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << "Added projection canvas to source frame!" << std::endl; }
528
530 std::array<int, 3> parts = {35, 50, 15};
531 fSourceStatusBar->SetParts(parts.data(), parts.size());
532
533 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << "Created status bar!" << std::endl; }
534
535 fSourceFrame->AddFrame(fSourceStatusBar, new TGLayoutHints(kLHintsBottom | kLHintsExpandX, 2, 2, 2, 2));
536
537 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << "Releasing unique lock on graphics mutex after building interface" << std::endl; }
538}
539
541{
542 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)
543 //std::unique_lock<std::mutex> graphicsLock(fSourceCalibration->GraphicsMutex());
544 fProjectionCanvas->GetCanvas()->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", "TSourceTab", this, "ProjectionStatus(Int_t,Int_t,Int_t,TObject*)");
545 fProjectionCanvas->Connect("ProcessedEvent(Event_t*)", "TSourceTab", this, "ProjectionStatus(Event_t*)");
546}
547
549{
550 //std::unique_lock<std::mutex> graphicsLock(fSourceCalibration->GraphicsMutex());
551 fProjectionCanvas->Disconnect("ProcessedEvent(Event_t*)", this, "ProjectionStatus(Event_t*)");
552 fProjectionCanvas->GetCanvas()->Disconnect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", this, "ProjectionStatus(Int_t,Int_t,Int_t,TObject*)");
553}
554
556{
557 // enum EGEventType {
558 // kGKeyPress = 0, kKeyRelease, = 1, kButtonPress = 2, kButtonRelease = 3,
559 // kMotionNotify = 4, kEnterNotify = 5, kLeaveNotify = 6, kFocusIn = 7, kFocusOut = 8,
560 // kExpose = 9, kConfigureNotify = 10, kMapNotify = 11, kUnmapNotify = 12, kDestroyNotify = 13,
561 // kClientMessage = 14, kSelectionClear = 15, kSelectionRequest = 16, kSelectionNotify = 17,
562 // kColormapNotify = 18, kButtonDoubleClick = 19, kOtherEvent = 20
563 // };
564 //if(TSourceCalibration::VerboseLevel() == EVerbosity::kAll) {
565 // 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)
566 // if(event->fType == kClientMessage) {
567 // std::cout << "Client message: " << event->fUser[0] << ", " << event->fUser[1] << ", " << event->fUser[2] << std::endl;
568 // }
569 //}
570 if(event->fType == kEnterNotify) {
572 std::cout << "Entered source tab => changing gPad from " << gPad->GetName();
573 }
574 gPad = fProjectionCanvas->GetCanvas();
576 std::cout << " to " << gPad->GetName() << std::endl;
577 }
578 }
579}
580
581void TSourceTab::ProjectionStatus(Int_t event, Int_t px, Int_t py, TObject* selected)
582{
583 // 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
584 // kButton1 = left mouse button, kButton2 = right mouse button
585 // enum EEventType {
586 // kNoEvent = 0,
587 // kButton1Down = 1, kButton2Down = 2, kButton3Down = 3, kKeyDown = 4,
588 // kWheelUp = 5, kWheelDown = 6, kButton1Shift = 7, kButton1ShiftMotion = 8,
589 // kButton1Up = 11, kButton2Up = 12, kButton3Up = 13, kKeyUp = 14,
590 // kButton1Motion = 21, kButton2Motion = 22, kButton3Motion = 23, kKeyPress = 24,
591 // kArrowKeyPress = 25, kArrowKeyRelease = 26,
592 // kButton1Locate = 41, kButton2Locate = 42, kButton3Locate = 43, kESC = 27,
593 // kMouseMotion = 51, kMouseEnter = 52, kMouseLeave = 53,
594 // kButton1Double = 61, kButton2Double = 62, kButton3Double = 63
595 //};
596 Status(selected->GetName(), 0);
597 Status(selected->GetObjectInfo(px, py), 1);
598 //auto key = static_cast<char>(px);
599 switch(event) {
600 case kButton1Down:
601 if(selected == fProjection) {
602 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; }
604 std::cout << fPeaks.size() << " peaks:";
605 for(auto* peak : fPeaks) {
606 std::cout << " " << peak->Centroid();
607 }
608 std::cout << std::endl;
609 }
610 auto* polym = static_cast<TPolyMarker*>(fProjection->GetListOfFunctions()->FindObject("TPolyMarker"));
611 if(polym == nullptr) {
612 std::cerr << "No peaks defined yet?" << std::endl;
613 return;
614 }
615 polym->SetNextPoint(fProjectionCanvas->GetCanvas()->AbsPixeltoX(px), fProjectionCanvas->GetCanvas()->AbsPixeltoX(py));
616 double range = 4 * fSourceCalibration->Sigma(); // * fProjection->GetXaxis()->GetBinWidth(1);
617 auto* pf = new TPeakFitter(fProjectionCanvas->GetCanvas()->AbsPixeltoX(px) - range, fProjectionCanvas->GetCanvas()->AbsPixeltoX(px) + range);
618 auto* peak = new TGauss(fProjectionCanvas->GetCanvas()->AbsPixeltoX(px));
619 pf->AddPeak(peak);
620 pf->Fit(fProjection, "qnretryfit");
621 if(peak->Area() > 0) {
622 fPeaks.push_back(peak);
623 fFits.push_back(pf);
624 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << "Fitted peak " << fProjectionCanvas->GetCanvas()->AbsPixeltoX(px) - range << " - " << fProjectionCanvas->GetCanvas()->AbsPixeltoX(px) + range << " -> centroid " << peak->Centroid() << std::endl; }
625 } else {
626 std::cout << "Ignoring peak at " << peak->Centroid() << " with negative area " << peak->Area() << std::endl;
627 }
628 fProjection->GetListOfFunctions()->Remove(peak);
629 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)
630
631 std::sort(fPeaks.begin(), fPeaks.end(), [](const TGauss* a, const TGauss* b) { return a->Centroid() < b->Centroid(); });
632
633 auto map = Match(fPeaks, fSourceEnergy, this);
634 Add(map);
635
636 // update status
637 Status(Form("%d/%d", static_cast<int>(fData->GetN()), static_cast<int>(fPeaks.size())), 2);
638
639 // update data and re-calibrate
643 }
644 break;
645 case kArrowKeyPress:
646 //Move1DHistogram(px, fProjection);
647 //fProjectionCanvas->GetCanvas()->Modified();
648 case kArrowKeyRelease:
649 break;
650 case kKeyDown:
651 case kKeyPress:
652 //switch(key) {
653 // case 'u':
654 // fProjection->GetXaxis()->UnZoom();
655 // fProjection->GetYaxis()->UnZoom();
656 // break;
657 // case 'U':
658 // fProjection->GetYaxis()->UnZoom();
659 // break;
660 // default:
661 // std::cout << "Key press '" << key << "' not recognized!" << std::endl;
662 // break;
663 //}
664 case kKeyUp:
665 break;
666 case kButton1Motion:
667 break;
668 case kButton1Up:
669 break;
670 case kMouseMotion:
671 //{
672 // auto* canvas = fProjectionCanvas->GetCanvas();
673 // if(canvas == nullptr) {
674 // Status("canvas nullptr", 0);
675 // break;
676 // }
677 // auto* pad = canvas->GetSelectedPad();
678 // if(pad == nullptr) {
679 // Status("pad nullptr", 0);
680 // break;
681 // }
682
683 // Status(Form("x = %f, y = %f", pad->AbsPixeltoX(px), pad->AbsPixeltoY(py)), 0);
684 //}
685 break;
686 case kMouseEnter:
687 break;
688 case kMouseLeave:
689 break;
690 default:
692 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;
693 }
694 break;
695 }
696}
697
699{
700 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << "UpdateRegions: clearing " << fRegions.size() << " regions" << std::endl; }
701 fRegions.clear();
703 std::cout << "projection has " << fProjection->ListOfRegions()->GetEntries() << " regions:" << std::endl;
704 fProjection->ListOfRegions()->Print();
705 }
707 for(auto* obj : *(fProjection->ListOfRegions())) {
708 if(obj->InheritsFrom(TRegion::Class())) {
709 auto* region = static_cast<TRegion*>(obj);
710 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; }
711 if(TSourceCalibration::VerboseLevel() > EVerbosity::kLoops) { region->Print(); }
712 fRegions.emplace_back(std::min(region->GetX1(), region->GetX2()), std::max(region->GetX1(), region->GetX2()));
714 std::cout << obj->GetName() << " is not a TRegion but a " << obj->ClassName() << std::endl;
715 }
716 }
717}
718
720{
721 return TSourceCalibration::AcceptBadFits() || (peak->CentroidErr() < 0.1 * peak->Centroid() && peak->AreaErr() < peak->Area() && !std::isnan(peak->CentroidErr()) && !std::isnan(peak->AreaErr()));
722}
723
725{
726 /// 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.
727 TList* functions = fProjection->GetListOfFunctions();
729 std::cout << "Updating functions for histogram " << fProjection->GetName() << " by deleting " << functions->GetEntries() << " functions:" << std::endl;
730 for(auto* obj : *functions) {
731 std::cout << obj << ": " << obj->GetName() << "/" << obj->GetTitle() << " - " << obj->ClassName();
732 if(obj->IsA() == TF1::Class()) {
733 std::cout << " line color " << static_cast<TF1*>(obj)->GetLineColor();
734 } else {
735 std::cout << " not a TF1";
736 }
737 std::cout << std::endl;
738 }
739 }
740 // remove all old fits (all have the same name/title of "gauss_total")
741 // alternative: clone the peaks with indices in their names, that would also allow us to add a function to replace a single fit
742 TObject* fit = nullptr;
743 while((fit = functions->FindObject("gauss_total")) != nullptr) {
744 functions->Remove(fit);
745 }
746 // add all the fit functions of the good and the bad peaks
747 for(auto* obj : fFits) {
748 functions->Add(obj->GetFitFunction()->Clone("gauss_total"));
749 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << "Added clone of " << obj->GetFitFunction() << " = " << functions->Last() << std::endl; }
750 }
751 for(auto* obj : fBadFits) {
752 functions->Add(obj->GetFitFunction()->Clone("gauss_total"));
753 }
755 std::cout << "And adding " << functions->GetEntries() << " functions instead:" << std::endl;
756 for(auto* obj : *functions) {
757 std::cout << obj << ": " << obj->GetName() << "/" << obj->GetTitle() << " - " << obj->ClassName();
758 if(obj->IsA() == TF1::Class()) {
759 std::cout << " line color " << static_cast<TF1*>(obj)->GetLineColor();
760 } else {
761 std::cout << " not a TF1";
762 }
763 std::cout << std::endl;
764 }
765 }
766 if(!fFits.empty() && fFits.front()->NumberOfPeaks() > 0 && fFits.back()->NumberOfPeaks() > 0) {
767 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()));
768 } else {
769 fProjection->SetTitle(Form("%lu of %lu source peaks used, channel ???? - ????", fFits.size(), fSourceEnergy.size()));
770 }
771 fProjectionCanvas->GetCanvas()->Modified();
772}
773
775{
776 fProjectionCanvas->GetCanvas()->cd();
777 fProjection->Draw(); // seems like hist + samefunc does not work. Could use hist + loop over list of functions (with same)
778 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << __PRETTY_FUNCTION__ << ": drew " << fProjection->GetName() << " on " << fProjectionCanvas->GetCanvas()->GetName() << "/" << gPad->GetName() << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
779}
780
781void TSourceTab::InitialCalibration(const bool& force)
782{
783 /// This functions finds the peaks in the histogram, fits them, and adds the fits to the list of peaks.
784 /// This list is then used to find all peaks that lie on a straight line.
785
786 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)
787
788 if(fPeaks.empty() || fData == nullptr || force) {
790 std::cout << __PRETTY_FUNCTION__ << ": # peaks " << fPeaks.size() << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
791 }
792 fPeaks.clear();
793 fFits.clear();
794 fBadFits.clear();
795 // Remove all associated functions from projection.
796 // These are the poly markers, fits, and the pave stat
798 TList* functions = fProjection->GetListOfFunctions();
799 std::cout << functions->GetEntries() << " functions in projection " << fProjection->GetName() << " before clearing" << std::endl;
800 for(auto* obj : *functions) {
801 std::cout << obj->GetName() << "/" << obj->GetTitle() << " - " << obj->ClassName() << std::endl;
802 }
803 }
804 fProjection->GetListOfFunctions()->Clear();
805 TSpectrum spectrum;
806 int nofPeaks = 0;
807 std::vector<double> peakPos;
809 if(fRegions.empty()) {
810 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << "No regions found, using whole spectrum" << std::endl; }
812 nofPeaks = spectrum.GetNPeaks();
813 peakPos.insert(peakPos.end(), spectrum.GetPositionX(), spectrum.GetPositionX() + nofPeaks);
814 } else {
815 for(auto& region : fRegions) {
816 fProjection->GetXaxis()->SetRangeUser(region.first, region.second);
818 int tmpNofPeaks = spectrum.GetNPeaks();
819 peakPos.insert(peakPos.end(), spectrum.GetPositionX(), spectrum.GetPositionX() + tmpNofPeaks);
820 nofPeaks += spectrum.GetNPeaks();
821 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)
822 }
823 fProjection->GetXaxis()->UnZoom();
824 }
825 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)
826
827 auto map = RoughCal(peakPos, fSourceEnergy, this);
828 Add(map);
829
830 // update status
831 Status(Form("%d/%d", static_cast<int>(fData->GetN()), nofPeaks), 2);
832 }
833}
834
835void TSourceTab::FindCalibratedPeaks(const TF1* calibration)
836{
837 /// This functions finds uses the existing calibration to find all peaks from the list of source energies.
838
839 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)
840
841 if(calibration == nullptr) {
842 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)
843 return;
844 }
845
846 // try and make sure the right canvas is selected as all fits are drawn in the currently selected canvas
847 gPad = fProjectionCanvas->GetCanvas();
848
849 fPeaks.clear();
850 fFits.clear();
851 fBadFits.clear();
852 // Remove all associated functions from projection.
853 // These are the poly markers, fits, and the pave stat
855 std::cout << "calibration has " << calibration->GetNpar() << " parameters (";
856 for(int i = 0; i < calibration->GetNpar(); ++i) {
857 if(i != 0) {
858 std::cout << ", ";
859 }
860 std::cout << i << ": " << calibration->GetParameter(i);
861 }
862 std::cout << "), if linear offset is " << calibration->Eval(0) << " and gain is " << calibration->Eval(1) - calibration->Eval(0) << std::endl;
863 TList* functions = fProjection->GetListOfFunctions();
864 std::cout << functions->GetEntries() << " functions in projection " << fProjection->GetName() << " before clearing" << std::endl;
865 for(auto* obj : *functions) {
866 std::cout << obj->GetName() << "/" << obj->GetTitle() << " - " << obj->ClassName() << std::endl;
867 }
868 }
869 fProjection->GetListOfFunctions()->Clear();
870
871 std::map<TGauss*, std::tuple<double, double, double, double>> map;
872
873 // loop over all source energies and check if we want to use the peak or if it's too weak
874 // for each source energy we calculate the channel it would be at
875 std::vector<std::pair<std::tuple<double, double, double, double>, double>> channelPos;
876 for(auto& tuple : fSourceEnergy) {
877 if(std::get<2>(tuple) < TSourceCalibration::MinIntensity()) {
879 std::cout << "Skipping peak at " << std::get<0>(tuple) << " keV with intensity " << std::get<2>(tuple) << " below required minimum intensity " << TSourceCalibration::MinIntensity() << std::endl;
880 }
881 continue;
882 }
883 channelPos.emplace_back(tuple, calibration->GetX(std::get<0>(tuple)));
885 std::cout << "Set channel position for " << channelPos.size() - 1 << ". peak to " << channelPos.back().second << " channels = " << std::get<0>(channelPos.back().first) << " keV" << std::endl;
886 }
887 }
888 // loop over channel positions and fit the peaks
889 for(size_t i = 0; i < channelPos.size(); ++i) {
891 double highRange = TSourceCalibration::FitRange() * fSourceCalibration->Sigma();
892 // if the range for this peaks overlaps with the previous/next one, reduce the range but not more than half
893 // also check if the previous/next one is above the minimum intensity, because otherwise we would have skipped that one anyway
894 if(i > 0) {
895 if(channelPos[i].second - lowRange < channelPos[i - 1].second + highRange) {
897 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;
898 }
899 lowRange = std::max(lowRange / 2., (channelPos[i].second - channelPos[i - 1].second) / 2.);
900 }
901 }
902 if(i + 1 < channelPos.size()) {
903 if(channelPos[i].second + highRange > channelPos[i + 1].second - lowRange) {
905 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;
906 }
907 highRange = std::min(highRange / 2., (channelPos[i + 1].second - channelPos[i].second) / 2.);
908 }
909 }
911 std::cout << "Trying to fit peak " << channelPos[i].second << " in range " << channelPos[i].second - lowRange << " - " << channelPos[i].second + highRange << std::endl;
912 }
913 // 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
914 auto* pf = new TPeakFitter(channelPos[i].second - lowRange, channelPos[i].second + highRange);
915 auto* peak = new TGauss(channelPos[i].second);
916 pf->AddPeak(peak);
917 if(TSourceCalibration::LogFile().empty()) {
918 TRedirect redirect("/dev/null");
919 pf->Fit(fProjection, "qnretryfit");
920 } else {
921 pf->Fit(fProjection, "nretryfit");
922 }
923 if(peak->Area() > 0) {
924 // check if we accept bad fits (and if it was a bad fit)
925 if(Good(peak, channelPos[i].second - lowRange, channelPos[i].second + highRange)) {
926 fPeaks.push_back(peak);
927 fFits.push_back(pf);
928 map[peak] = channelPos[i].first;
930 std::cout << "Fitted peak " << channelPos[i].second - lowRange << " - " << channelPos[i].second + highRange << " -> centroid " << peak->Centroid() << ", area " << peak->Area() << std::endl;
931 }
932 } else {
933 fBadFits.push_back(pf);
935 std::cout << "We're ignoring bad fits, and this one has position " << peak->Centroid() << " +- " << peak->CentroidErr() << ", area " << peak->Area() << " +- " << peak->AreaErr() << std::endl;
936 }
937 }
939 std::cout << "Ignoring peak at " << peak->Centroid() << " with negative area " << peak->Area() << std::endl;
940 }
941 }
942
943 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)
944 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)
945
946 std::sort(fFits.begin(), fFits.end(), [](const TPeakFitter* a, const TPeakFitter* b) { return a->GetFitFunction()->GetParameter("centroid_0") < b->GetFitFunction()->GetParameter("centroid_0"); });
947 std::sort(fBadFits.begin(), fBadFits.end(), [](const TPeakFitter* a, const TPeakFitter* b) { return a->GetFitFunction()->GetParameter("centroid_0") < b->GetFitFunction()->GetParameter("centroid_0"); });
948 std::sort(fPeaks.begin(), fPeaks.end(), [](const TGauss* a, const TGauss* b) { return a->Centroid() < b->Centroid(); });
949
951
952 UpdateFits();
953
954 Add(map);
955
957 std::cout << __PRETTY_FUNCTION__ << ": found " << fData->GetN() << " peaks"; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
958 Double_t* x = fData->GetX();
959 Double_t* y = fData->GetY();
960 for(int i = 0; i < fData->GetN(); ++i) {
961 std::cout << "; " << x[i] << " - " << y[i];
962 }
963 std::cout << std::endl;
964 }
965}
966
967void TSourceTab::Add(std::map<double, std::tuple<double, double, double, double>> map)
968{
970 std::cout << DCYAN << "Adding map with " << map.size() << " data points, fData = " << fData << std::endl;
971 }
972 if(fData == nullptr) {
973 fData = new TGraphErrors(map.size());
974 } else {
975 fData->Set(map.size());
976 }
977 fData->SetLineColor(2);
978 fData->SetMarkerColor(2);
979 int i = 0;
980 for(auto iter = map.begin(); iter != map.end();) {
981 // more readable variable names
982 auto peakPos = iter->first;
983 auto energy = std::get<0>(iter->second);
984 fData->SetPoint(i, peakPos, energy);
986 std::cout << "Using peak with position " << peakPos << ", energy " << energy << std::endl;
987 }
988 ++iter;
989 ++i;
990 }
992 std::cout << "Accepted " << map.size() << " peaks" << std::endl;
993 }
994 fData->Sort();
995}
996
997void TSourceTab::Add(std::map<TGauss*, std::tuple<double, double, double, double>> map)
998{
1000 std::cout << DCYAN << "Adding map with " << map.size() << " data points, fData = " << fData << std::endl;
1001 }
1002 if(fData == nullptr) {
1003 fData = new TGraphErrors(map.size());
1004 } else {
1005 fData->Set(map.size());
1006 }
1007 fData->SetLineColor(2);
1008 fData->SetMarkerColor(2);
1009 if(fFwhm == nullptr) {
1010 fFwhm = new TGraphErrors(map.size());
1011 } else {
1012 fFwhm->Set(map.size());
1013 }
1014 fFwhm->SetLineColor(4);
1015 fFwhm->SetMarkerColor(4);
1016 int i = 0;
1017 for(auto iter = map.begin(); iter != map.end();) {
1018 // more readable variable names
1019 auto peakPos = iter->first->Centroid();
1020 auto peakPosErr = iter->first->CentroidErr();
1021 auto peakArea = iter->first->Area();
1022 auto peakAreaErr = iter->first->AreaErr();
1023 auto fwhm = iter->first->FWHM();
1024 auto fwhmErr = iter->first->FWHMErr();
1025 auto energy = std::get<0>(iter->second);
1026 auto energyErr = std::get<1>(iter->second);
1027 // drop this peak if the uncertainties in area or position are too large
1028 if(!Good(iter->first)) {
1030 std::cout << "Dropping peak with position " << peakPos << " +- " << peakPosErr << ", area " << peakArea << " +- " << peakAreaErr << ", energy " << energy << std::endl;
1031 }
1032 iter = map.erase(iter);
1034 std::cout << "Erasing peak returned iterator " << std::distance(map.begin(), iter) << std::endl;
1035 }
1036 } else {
1037 fData->SetPoint(i, peakPos, energy);
1038 fData->SetPointError(i, peakPosErr, energyErr);
1039 fFwhm->SetPoint(i, peakPos, fwhm);
1040 fFwhm->SetPointError(i, peakPosErr, fwhmErr);
1042 std::cout << "Using peak with position " << peakPos << " +- " << peakPosErr << ", area " << peakArea << " +- " << peakAreaErr << ", energy " << energy << std::endl;
1043 }
1044 ++iter;
1045 ++i;
1046 }
1047 }
1049 std::cout << "Accepted " << map.size() << " peaks" << std::endl;
1050 }
1051 // if we dropped a peak, we need to resize the graph, if we haven't dropped any this doesn't do anything
1052 fData->Set(map.size());
1053 fFwhm->Set(map.size());
1054 fData->Sort();
1055 fFwhm->Sort();
1056 // split poly marker into those peaks that were used, and those that weren't
1057 TList* functions = fProjection->GetListOfFunctions();
1059 std::cout << functions->GetEntries() << " functions in projection " << fProjection->GetName() << std::endl;
1060 for(auto* obj : *functions) {
1061 std::cout << obj->GetName() << "/" << obj->GetTitle() << " - " << obj->ClassName() << std::endl;
1062 }
1063 }
1064 auto* polym = static_cast<TPolyMarker*>(functions->FindObject("TPolyMarker"));
1065 if(polym != nullptr) {
1066 double* oldX = polym->GetX();
1067 double* oldY = polym->GetY();
1068 int size = polym->GetN();
1069 std::vector<double> newX;
1070 std::vector<double> newY;
1071 std::vector<double> unusedX;
1072 std::vector<double> unusedY;
1073 for(i = 0; i < size; ++i) {
1074 bool used = false;
1075 for(auto point : map) {
1076 if(TMath::Abs(oldX[i] - point.first->Centroid()) < fSourceCalibration->Sigma()) {
1077 newX.push_back(oldX[i]);
1078 newY.push_back(oldY[i]);
1079 used = true;
1080 break;
1081 }
1082 }
1083 if(!used) {
1084 unusedX.push_back(oldX[i]);
1085 unusedY.push_back(oldY[i]);
1086 }
1087 }
1088 polym->SetPolyMarker(newX.size(), newX.data(), newY.data());
1089 auto* unusedMarkers = new TPolyMarker(unusedX.size(), unusedX.data(), unusedY.data());
1090 unusedMarkers->SetMarkerStyle(23); // triangle down
1091 unusedMarkers->SetMarkerColor(16); // light grey
1092 functions->Add(unusedMarkers);
1094 std::cout << fProjection->GetTitle() << ": " << size << " peaks founds total, " << polym->GetN() << " used, and " << unusedMarkers->GetN() << " unused" << std::endl;
1096 std::cout << "Used: ";
1097 for(size_t m = 0; m < newX.size(); ++m) { std::cout << newX[m] << " - " << newY[m] << ";"; }
1098 std::cout << std::endl;
1099 std::cout << "Unused: ";
1100 for(size_t m = 0; m < unusedX.size(); ++m) { std::cout << unusedX[m] << " - " << unusedY[m] << ";"; }
1101 std::cout << std::endl;
1102 }
1104 std::cout << functions->GetEntries() << " functions in projection " << fProjection->GetTitle() << std::endl;
1105 for(auto* obj : *functions) {
1106 std::cout << obj->GetName() << "/" << obj->GetTitle() << " - " << obj->ClassName() << std::endl;
1107 }
1108 }
1109 }
1110 }
1111
1112 // change color of fit functions for unused peaks
1113 TIter iter(functions);
1114 TObject* item = nullptr;
1115 while((item = iter.Next()) != nullptr) {
1116 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
1117 double centroid = 0.;
1118 if(item->IsA() == TF1::Class()) {
1119 centroid = static_cast<TF1*>(item)->GetParameter(1);
1120 } else {
1121 centroid = static_cast<TGauss*>(item)->Centroid();
1122 }
1123 bool found = false;
1124 for(auto point : map) {
1125 if(TMath::Abs(centroid - point.first->Centroid()) < fSourceCalibration->Sigma()) {
1126 found = true;
1127 break;
1128 }
1129 }
1130 // remove the TF1 if it wasn't found in the map
1131 if(!found) {
1132 //functions->Remove(item);
1133 if(item->IsA() == TF1::Class()) {
1134 static_cast<TF1*>(item)->SetLineColor(kGray);
1135 } else {
1136 static_cast<TGauss*>(item)->GetFitFunction()->SetLineColor(kGray);
1137 }
1138 }
1139 }
1140 }
1141}
1142
1143void TSourceTab::ReplacePeak(const size_t& index, const double& channel)
1144{
1145 /// Replace the peak at the index with one centered at channel (calculated from an initial calibration and the source energy of this peak).
1146 /// This replaces the TGauss* in fPeaks, and updates the values of this index in fData and fFwhm.
1148 std::cout << DCYAN << __PRETTY_FUNCTION__ << ": index " << index << ", channel " << channel << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
1149 for(size_t i = 0; i < fPeaks.size() && i < fFits.size(); ++i) {
1150 std::cout << i << ": " << fPeaks[i]->Centroid() << " +- " << fPeaks[i]->CentroidErr() << " - " << fFits[i]->GetFitFunction()->GetParameter("centroid_0") << std::endl;
1151 }
1152 }
1153
1154 if(index >= fPeaks.size()) {
1155 std::cerr << "Trying to replace peak at index " << index << " but there are only " << fPeaks.size() << " peaks!" << std::endl;
1156 return;
1157 }
1158
1159 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; }
1160
1161 // calculate the fitting range (low to high) and adjust it so we ignore the old peak
1163 double low = channel - range;
1164 double high = channel + range;
1165 if(fPeaks.at(index)->Centroid() < channel) {
1166 low = (fPeaks.at(index)->Centroid() + channel) / 2.;
1167 } else {
1168 high = (fPeaks.at(index)->Centroid() + channel) / 2.;
1169 }
1170 // 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
1171 auto* pf = new TPeakFitter(low, high);
1172 auto* peak = new TGauss(channel);
1173 pf->AddPeak(peak);
1174 if(TSourceCalibration::LogFile().empty()) {
1175 TRedirect redirect("/dev/null");
1176 pf->Fit(fProjection, "retryfit");
1177 } else {
1178 pf->Fit(fProjection, "retryfit");
1179 }
1180 if(peak->Area() > 0) {
1181 if(Good(peak, low, high)) {
1183 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;
1184 std::cout << "Deleting " << index << ". peak " << fPeaks.at(index) << std::flush;
1185 }
1186 peak->SetLineColor(fPeaks.at(index)->GetLineColor());
1187 delete fPeaks.at(index);
1188 fPeaks[index] = peak;
1190 std::cout << " and setting it to " << fPeaks.at(index) << std::endl;
1191 }
1193 std::cout << "Deleting " << index << ". fit " << fFits.at(index) << std::flush;
1194 }
1195 pf->GetFitFunction()->SetLineColor(fFits.at(index)->GetFitFunction()->GetLineColor());
1196 delete fFits.at(index);
1197 fFits[index] = pf;
1199 std::cout << " and setting it to " << fFits.at(index) << std::endl;
1200 }
1201 UpdateFits();
1202 fData->SetPoint(index, peak->Centroid(), fData->GetPointY(index));
1203 fData->SetPointError(index, peak->CentroidErr(), fData->GetErrorY(index));
1204 fFwhm->SetPoint(index, peak->Centroid(), peak->FWHM());
1205 fFwhm->SetPointError(index, peak->CentroidErr(), peak->FWHMErr());
1207 std::cout << "Fitted peak " << low << " - " << high << " -> centroid " << peak->Centroid() << ", area " << peak->Area() << std::endl;
1208 }
1209 } else {
1211 std::cout << "We're ignoring bad fits, and this one has position " << peak->Centroid() << " +- " << peak->CentroidErr() << ", area " << peak->Area() << " +- " << peak->AreaErr() << std::endl;
1212 }
1213 }
1215 std::cout << "Ignoring peak at " << peak->Centroid() << " with negative area " << peak->Area() << std::endl;
1216 }
1217}
1218
1219void TSourceTab::RemovePoint(Int_t oldPoint)
1220{
1222 std::cout << DCYAN << __PRETTY_FUNCTION__ << ": oldPoint " << oldPoint << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
1223 }
1224 // first thing to do is check if the graph the point was removed from matches the id of this tab
1226 //auto erasedPeak = fPeaks.at(oldPoint);
1227 fPeaks.erase(fPeaks.begin() + oldPoint);
1228 auto* erasedFit = fFits.at(oldPoint);
1229 fFits.erase(fFits.begin() + oldPoint);
1230 fBadFits.push_back(erasedFit);
1231 fData->RemovePoint(oldPoint);
1232
1233 SetLineColors();
1234
1235 UpdateFits();
1236
1238}
1239
1241{
1242 /// This function sets the line colors of good and bad fits to alternating colors.
1243 /// Red and blue for the good fits, and grey and dark grey for the bad fits.
1244
1245 // set colors of good fits alternating to red or blue (peaks themselves are set to darker versions of those colors)
1246 for(size_t i = 0; i < fPeaks.size() && i < fFits.size(); ++i) {
1247 if(i % 2 == 0) {
1248 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; }
1249 fPeaks[i]->SetLineColor(kRed + 2);
1250 fFits[i]->GetFitFunction()->SetLineColor(kRed);
1251 } else {
1252 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; }
1253 fPeaks[i]->SetLineColor(kBlue + 2);
1254 fFits[i]->GetFitFunction()->SetLineColor(kBlue);
1255 }
1256 }
1257
1258 // set colors of bad fits alterating to grey and dark grey
1259 for(size_t i = 0; i < fBadFits.size(); ++i) {
1260 if(i % 2 == 0) {
1261 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << i << ": setting line colors to grey for fit function " << fBadFits[i]->GetFitFunction() << std::endl; }
1262 fBadFits[i]->GetFitFunction()->SetLineColor(kGray);
1263 } else {
1264 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << i << ": setting line colors to dark grey for fit function " << fBadFits[i]->GetFitFunction() << std::endl; }
1265 fBadFits[i]->GetFitFunction()->SetLineColor(kGray + 2);
1266 }
1267 }
1268}
1269
1270void TSourceTab::Status(const char* status, int position)
1271{
1272 fSourceStatusBar->SetText(status, position);
1273 fSourceStatusBar->MapWindow();
1274 fSourceFrame->MapSubwindows();
1275 gVirtualX->Update();
1276}
1277
1279{
1280 std::cout << "Projection " << fProjection;
1281 if(fProjection != nullptr) {
1282 std::cout << " = " << fProjection->GetName() << "/" << fProjection->GetTitle();
1283 } else {
1284 std::cout << " is null";
1285 }
1286 std::cout << std::endl;
1287
1288 std::cout << "Data " << fData;
1289 if(fData != nullptr) {
1290 std::cout << " = " << fData->GetName() << "/" << fData->GetTitle();
1291 } else {
1292 std::cout << " is null";
1293 }
1294 std::cout << std::endl;
1295
1296 std::cout << "Fwhm " << fFwhm;
1297 if(fFwhm != nullptr) {
1298 std::cout << " = " << fFwhm->GetName() << "/" << fFwhm->GetTitle();
1299 } else {
1300 std::cout << " is null";
1301 }
1302 std::cout << std::endl;
1303
1304 std::cout << fFits.size() << " good fits:";
1305 for(size_t i = 0; i < fFits.size(); ++i) {
1306 std::cout << " " << std::setw(2) << i << " - " << std::setw(8) << fFits.at(i)->GetFitFunction()->GetParameter("centroid_0");
1307 }
1308 std::cout << std::endl;
1309
1310 std::cout << fBadFits.size() << " bad fits:";
1311 for(size_t i = 0; i < fBadFits.size(); ++i) {
1312 std::cout << " " << std::setw(2) << i << " - " << std::setw(8) << fBadFits.at(i)->GetFitFunction()->GetParameter("centroid_0");
1313 }
1314 std::cout << std::endl;
1315
1316 std::cout << fPeaks.size() << " peaks:";
1317 for(size_t i = 0; i < fPeaks.size(); ++i) {
1318 std::cout << " " << std::setw(2) << i << " - " << std::setw(8) << fPeaks.at(i)->Centroid();
1319 }
1320 std::cout << std::endl;
1321}
1322
1324{
1325 std::cout << "TSourceTab " << fSourceName << " projection canvas:" << std::endl;
1326 fProjectionCanvas->GetCanvas()->Print();
1327 for(auto* obj : *(fProjectionCanvas->GetCanvas()->GetListOfPrimitives())) {
1328 obj->Print();
1329 }
1330}
1331
1333{
1334 std::cout << "TSourceTab frame:" << std::endl;
1335 fSourceFrame->Print();
1336 std::cout << "TSourceTab canvas:" << std::endl;
1337 fProjectionCanvas->Print();
1338 std::cout << "TSourceTab status bar:" << std::endl;
1339 fSourceStatusBar->Print();
1340}
1341
1342//////////////////////////////////////// TChannelTab ////////////////////////////////////////
1343TChannelTab::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)
1344 : 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))
1345{
1347 std::cout << DYELLOW << "========================================" << std::endl;
1348 std::cout << __PRETTY_FUNCTION__ << std::endl // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
1349 << " name = " << fName << ", fData = " << fData << std::endl;
1350 std::cout << "========================================" << std::endl;
1351 }
1352
1354
1356 std::cout << DYELLOW << "----------------------------------------" << std::endl;
1357 std::cout << __PRETTY_FUNCTION__ << std::endl // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
1358 << " channel " << fName << " done" << std::endl;
1359 std::cout << "----------------------------------------" << std::endl;
1360 }
1361}
1362
1364{
1365 delete fChannelFrame;
1366 delete fSourceTab;
1367 for(auto* channel : fSources) {
1368 delete channel;
1369 }
1370 delete fCanvasTab;
1371 delete fCalibrationFrame;
1372 delete fCalibrationCanvas;
1373 delete fResidualPad;
1374 delete fCalibrationPad;
1375 delete fChannelStatusBar;
1376 delete fLegend;
1377 delete fChi2Label;
1378 delete fFwhmFrame;
1379 delete fFwhmCanvas;
1380 // delete all fProjections and all fNuclei?
1381 delete fData;
1382 delete fFwhm;
1383}
1384
1386{
1387 {
1388 //std::unique_lock<std::mutex> graphicsLock(fSourceCalibration->GraphicsMutex());
1389 fChannelFrame->SetLayoutManager(new TGHorizontalLayout(fChannelFrame));
1390 }
1391
1392 fSources.resize(fNuclei.size(), nullptr);
1393 //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)
1394 for(size_t source = 0; source < fNuclei.size(); ++source) {
1395 CreateSourceTab(source);
1396 }
1397
1398 for(auto iter = fSources.begin(); iter != fSources.end(); ++iter) {
1399 if(*iter == nullptr) {
1400 fSources.erase(iter--); // erase iterator and then go to the element before this one (and then the loop goes to the next one)
1401 }
1402 }
1403
1404 {
1405 //std::unique_lock<std::mutex> graphicsLock(fSourceCalibration->GraphicsMutex());
1406 fChannelFrame->AddFrame(fSourceTab, new TGLayoutHints(kLHintsLeft | kLHintsTop | kLHintsExpandX | kLHintsExpandY, 2, 2, 2, 2));
1407
1408 fCalibrationFrame = fCanvasTab->AddTab("Calibration");
1409 fCalibrationFrame->SetLayoutManager(new TGVerticalLayout(fCalibrationFrame));
1410 fCalibrationCanvas = new TRootEmbeddedCanvas(Form("CalibrationCanvas%s", fName.c_str()), fCalibrationFrame, TSourceCalibration::PanelWidth(), TSourceCalibration::PanelHeight());
1411 fCalibrationFrame->AddFrame(fCalibrationCanvas, new TGLayoutHints(kLHintsTop | kLHintsExpandX | kLHintsExpandY, 2, 2, 2, 2));
1413 std::array<int, 3> parts = {25, 35, 40};
1414 fChannelStatusBar->SetParts(parts.data(), parts.size());
1415 fCalibrationFrame->AddFrame(fChannelStatusBar, new TGLayoutHints(kLHintsBottom | kLHintsExpandX, 2, 2, 2, 2));
1416 fFwhmFrame = fCanvasTab->AddTab("FWHM");
1417 fFwhmFrame->SetLayoutManager(new TGVerticalLayout(fFwhmFrame));
1418 fFwhmCanvas = new TRootEmbeddedCanvas(Form("FwhmCanvas%s", fName.c_str()), fFwhmFrame, TSourceCalibration::PanelWidth(), TSourceCalibration::PanelHeight());
1419 fFwhmFrame->AddFrame(fFwhmCanvas, new TGLayoutHints(kLHintsTop | kLHintsExpandX | kLHintsExpandY, 2, 2, 2, 2));
1420 //fChannelFrame->AddFrame(fCalibrationFrame, new TGLayoutHints(kLHintsRight | kLHintsBottom | kLHintsExpandX | kLHintsExpandY, 2, 2, 2, 2));
1421
1422 fInitFrame = fCanvasTab->AddTab("Initial");
1423 fInitFrame->SetLayoutManager(new TGVerticalLayout(fInitFrame));
1424 fInitCanvas = new TRootEmbeddedCanvas(Form("InitCanvas%s", fName.c_str()), fInitFrame, TSourceCalibration::PanelWidth(), TSourceCalibration::PanelHeight());
1425 fInitFrame->AddFrame(fInitCanvas, new TGLayoutHints(kLHintsTop | kLHintsExpandX | kLHintsExpandY, 2, 2, 2, 2));
1426
1427 fChannelFrame->AddFrame(fCanvasTab, new TGLayoutHints(kLHintsRight | kLHintsTop | kLHintsExpandX | kLHintsExpandY, 2, 2, 2, 2));
1428 }
1429
1430 UpdateData();
1431
1432 {
1433 //std::unique_lock<std::mutex> graphicsLock(fSourceCalibration->GraphicsMutex());
1434 fCalibrationCanvas->GetCanvas()->cd();
1435 fCalibrationPad = new TPad(Form("cal_%s", fName.c_str()), Form("calibration for %s", fName.c_str()), 0.2, 0., 1., 1.);
1436 fCalibrationPad->SetNumber(1);
1437 fCalibrationPad->Draw();
1438 fCalibrationPad->AddExec("zoom", "TChannelTab::ZoomY()");
1439
1440 fLegend = new TLegend(0.8, 0.3, 0.95, 0.3 + static_cast<double>(fNuclei.size()) * 0.05); // x1, y1, x2, y2
1441
1442 fCalibrationCanvas->GetCanvas()->cd();
1443 fResidualPad = new TPad(Form("res_%s", fName.c_str()), Form("residual for %s", fName.c_str()), 0.0, 0., 0.2, 1.);
1444 fResidualPad->SetNumber(2);
1445 fResidualPad->Draw();
1446 fResidualPad->AddExec("zoom", "TChannelTab::ZoomY()");
1447 }
1448 //Calibrate(); // also creates the residual and chi^2 label
1449
1450 // check whether we want to use this initial calibration to find more peaks
1451 //if(TSourceCalibration::UseCalibratedPeaks()) {
1452 // FindCalibratedPeaks();
1453 //}
1454
1455 // get the fwhm graphs and plot them
1456 UpdateFwhm();
1457 //{
1458 // std::unique_lock<std::mutex> graphicsLock(fSourceCalibration->GraphicsMutex());
1459 // fFwhmCanvas->GetCanvas()->cd();
1460 // fFwhm->DrawCalibration("*");
1461 // fLegend->Draw();
1462 //}
1463}
1464
1466{
1467 if(fProjections[source]->GetEntries() > 1000) {
1469 std::cout << DYELLOW << "Creating source tab " << source << ", using fSourceCalibration " << fSourceCalibration << std::endl;
1470 }
1471 TGCompositeFrame* tmpTab = nullptr;
1472 {
1473 //if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << "Creating unique lock on graphics mutex to add tab" << std::endl; }
1474 //std::unique_lock<std::mutex> graphicsLock(fSourceCalibration->GraphicsMutex());
1475 tmpTab = fSourceTab->AddTab(Form("%s_%s", fNuclei[source]->GetName(), fName.c_str()));
1476 //if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << "Releasing unique lock on graphics mutex after adding tab" << std::endl; }
1477 }
1478 fSources[source] = new TSourceTab(fSourceCalibration, this, tmpTab, fProjections[source], fNuclei[source]->GetName(), fSourceEnergies[source]);
1479 {
1480 //if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << "Creating unique lock on graphics mutex to increment status bar" << std::endl; }
1481 //std::unique_lock<std::mutex> graphicsLock(fSourceCalibration->GraphicsMutex());
1482 fProgressBar->Increment(1);
1483 //if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << "Releasing unique lock on graphics mutex after incrementing status bar" << std::endl; }
1484 }
1485 } else {
1486 fSources[source] = nullptr;
1488 std::cout << DYELLOW << "Skipping projection of source " << source << " = " << fProjections[source]->GetName() << ", only " << fProjections[source]->GetEntries() << " entries" << std::endl;
1489 }
1490 }
1492 std::cout << __PRETTY_FUNCTION__ << " source " << source << " done" << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
1493 }
1494}
1495
1497{
1498 //std::unique_lock<std::mutex> graphicsLock(fSourceCalibration->GraphicsMutex());
1499 fSourceTab->Connect("Selected(Int_t)", "TChannelTab", this, "SelectedTab(Int_t)");
1500 for(auto* source : fSources) {
1501 source->MakeConnections();
1502 }
1503 fCalibrationCanvas->GetCanvas()->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", "TChannelTab", this, "CalibrationStatus(Int_t,Int_t,Int_t,TObject*)");
1504 fCalibrationCanvas->Connect("ProcessedEvent(Event_t*)", "TChannelTab", this, "SelectCanvas(Event_t*)");
1505 if(fData != nullptr) {
1507 std::cout << DYELLOW << __PRETTY_FUNCTION__ << ": connecting fData " << fData << ":" << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
1508 fData->Print();
1509 }
1510 fData->Connect("RemovePoint(Int_t, Int_t)", "TChannelTab", this, "RemovePoint(Int_t, Int_t)");
1511 } else {
1512 std::cerr << DRED << "Failed to create connections for the data graph fData since it hasn't been created yet!" << std::endl;
1513 }
1514}
1515
1517{
1518 //std::unique_lock<std::mutex> graphicsLock(fSourceCalibration->GraphicsMutex());
1519 fSourceTab->Disconnect("Selected(Int_t)", this, "SelectedTab(Int_t)");
1520 for(auto* source : fSources) {
1521 source->Disconnect();
1522 }
1523 fCalibrationCanvas->GetCanvas()->Disconnect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", this, "CalibrationStatus(Int_t,Int_t,Int_t,TObject*)");
1524 fData->Disconnect("RemovePoint(Int_t, Int_t)", this, "RemovePoint(Int_t, Int_t)");
1525}
1526
1528{
1529 /// Simple function that enables and disables the previous and next buttons depending on which tab was selected.
1530 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)
1531 fActiveSourceTab = id;
1532 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << "active source tab " << fActiveSourceTab << RESET_COLOR << std::endl; }
1533 if(fSources[fActiveSourceTab] == nullptr) {
1534 std::cout << "Failed to get fSources[" << fActiveSourceTab << "]" << RESET_COLOR << std::endl;
1535 return;
1536 }
1537 if(fSources[fActiveSourceTab]->ProjectionCanvas() == nullptr) {
1538 std::cout << "Failed to get fSources[" << fActiveSourceTab << "]->ProjectionCanvas()" << RESET_COLOR << std::endl;
1539 return;
1540 }
1541 if(fSources[fActiveSourceTab]->ProjectionCanvas()->GetCanvas() == nullptr) {
1542 std::cout << "Failed to get fSources[" << fActiveSourceTab << "]->ProjectionCanvas()->GetCanvas()" << RESET_COLOR << std::endl;
1543 return;
1544 }
1545 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << DGREEN << "Changing gpad from \"" << gPad->GetName() << "\" to \""; }
1546 gPad = fSources[fActiveSourceTab]->ProjectionCanvas()->GetCanvas();
1547 fSources[fActiveSourceTab]->ProjectionCanvas()->GetCanvas()->cd();
1548 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << gPad->GetName() << "\"" << RESET_COLOR << std::endl; }
1549}
1550
1551void TChannelTab::SelectCanvas(Event_t* event)
1552{
1553 if(event->fType == kEnterNotify) {
1555 std::cout << "Entered channel tab => changing gPad from " << gPad->GetName();
1556 }
1557 gPad = fCalibrationCanvas->GetCanvas();
1559 std::cout << " to " << gPad->GetName() << std::endl;
1560 }
1561 }
1562}
1563
1564void TChannelTab::RemovePoint(Int_t oldGraph, Int_t oldPoint)
1565{
1567 std::cout << DCYAN << __PRETTY_FUNCTION__ << ": oldGraph " << oldGraph << ", oldPoint " << oldPoint << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
1568 }
1569 if(0 <= oldGraph && oldGraph < static_cast<int>(fSources.size())) {
1570 if(oldPoint >= 0) {
1571 fSources[oldGraph]->RemovePoint(oldPoint);
1572 return;
1573 }
1574 std::cout << "Can't remove negative point " << oldPoint << " from graph " << oldGraph << std::endl;
1575 }
1576 std::cout << "Graph the point was removed from was " << oldGraph << ", but we only have " << fSources.size() << " source tabs?" << std::endl;
1577}
1578
1579void TChannelTab::CalibrationStatus(Int_t event, Int_t px, Int_t py, TObject* selected)
1580{
1581 // 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
1582 // kButton1 = left mouse button, kButton2 = right mouse button
1583 // enum EEventType {
1584 // kNoEvent = 0,
1585 // kButton1Down = 1, kButton2Down = 2, kButton3Down = 3, kKeyDown = 4,
1586 // kWheelUp = 5, kWheelDown = 6, kButton1Shift = 7, kButton1ShiftMotion = 8,
1587 // kButton1Up = 11, kButton2Up = 12, kButton3Up = 13, kKeyUp = 14,
1588 // kButton1Motion = 21, kButton2Motion = 22, kButton3Motion = 23, kKeyPress = 24,
1589 // kArrowKeyPress = 25, kArrowKeyRelease = 26,
1590 // kButton1Locate = 41, kButton2Locate = 42, kButton3Locate = 43, kESC = 27,
1591 // kMouseMotion = 51, kMouseEnter = 52, kMouseLeave = 53,
1592 // kButton1Double = 61, kButton2Double = 62, kButton3Double = 63
1593 //};
1594 fChannelStatusBar->SetText(selected->GetName(), 0);
1595 if(event == kKeyPress) {
1596 fChannelStatusBar->SetText(Form("%c", static_cast<char>(px)), 1);
1597 } else {
1598 auto* canvas = fCalibrationCanvas->GetCanvas();
1599 if(canvas == nullptr) {
1600 fChannelStatusBar->SetText("canvas nullptr");
1601 return;
1602 }
1603 auto* pad = canvas->GetSelectedPad();
1604 if(pad == nullptr) {
1605 fChannelStatusBar->SetText("pad nullptr");
1606 return;
1607 }
1608 //if(TSourceCalibration::VerboseLevel() > EVerbosity::kLoops) {
1609 // 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;
1610 //}
1611
1612 fChannelStatusBar->SetText(Form("x = %f, y = %f", pad->AbsPixeltoX(px), pad->AbsPixeltoY(py)), 1);
1613 }
1614}
1615
1617{
1618 /// Copy data from all sources into one graph (which we use for the calib && source < fSources.size()ration).
1619 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)
1620 if(fData == nullptr) {
1621 fData = new TCalibrationGraphSet("ADC channel", "Energy [keV]");
1622 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)
1623 }
1624 fData->Clear();
1625
1626 fData->SetName(Form("data%s", fName.c_str()));
1628 std::cout << "set name of fData using " << fName.c_str() << ": " << fData->GetName() << std::endl;
1629 std::cout << "fData " << fData << ": " << (fData != nullptr ? fData->GetN() : -2) << " data points after creation" << std::endl;
1630 std::cout << "Looping over " << fNuclei.size() << " sources, fSources.size() = " << fSources.size() << std::endl;
1631 }
1632 for(size_t source = 0; source < fNuclei.size() && source < fSources.size(); ++source) {
1633 if(fSources[source]->Data() != nullptr && fSources[source]->Data()->GetN() > 0) {
1634 int index = fData->Add(fSources[source]->Data(), fNuclei[source]->GetName());
1635 if(index >= 0) {
1636 fData->SetLineColor(index, static_cast<Color_t>(source + 1)); //+1 for the color so that we start with 1 = black instead of 0 = white
1637 fData->SetMarkerColor(index, static_cast<Color_t>(source + 1));
1638 }
1639 }
1640 }
1642 std::cout << "fData " << fData << ": " << (fData != nullptr ? fData->GetN() : -1) << " data points after adding " << fNuclei.size() << " graphs" << std::endl;
1643 fData->Print();
1644 }
1645}
1646
1648{
1649 /// Copy data from all sources into one graph (which we use for the calib && source < fSources.size()ration).
1650 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)
1651 if(fFwhm == nullptr) {
1652 fFwhm = new TCalibrationGraphSet("ADC channel", "FWHM [ADC channel]");
1653 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)
1654 }
1655 fFwhm->Clear();
1656
1657 fFwhm->SetName(Form("fwhm%s", fName.c_str()));
1659 std::cout << "set name of fFwhm using " << fName.c_str() << ": " << fFwhm->GetName() << std::endl;
1660 std::cout << "fFwhm " << fFwhm << ": " << (fFwhm != nullptr ? fFwhm->GetN() : -2) << " data points after creation" << std::endl;
1661 std::cout << "Looping over " << fNuclei.size() << " sources, fSources.size() = " << fSources.size() << std::endl;
1662 }
1663 for(size_t source = 0; source < fNuclei.size() && source < fSources.size(); ++source) {
1664 if(fSources[source]->Fwhm() != nullptr && fSources[source]->Fwhm()->GetN() > 0) {
1665 int index = fFwhm->Add(fSources[source]->Fwhm(), fNuclei[source]->GetName());
1666 if(index >= 0) {
1667 fFwhm->SetLineColor(index, static_cast<Color_t>(source + 1)); //+1 for the color so that we start with 1 = black instead of 0 = white
1668 fFwhm->SetMarkerColor(index, static_cast<Color_t>(source + 1));
1669 }
1670 }
1671 }
1673 std::cout << "fFwhm " << fFwhm << ": " << (fFwhm != nullptr ? fFwhm->GetN() : -1) << " data points after adding " << fNuclei.size() << " graphs" << std::endl;
1674 fFwhm->Print();
1675 }
1676}
1677
1679{
1680 /// Uses the statis TSourceCalibration::Degree for the degree of the polynomial used to calibrate the data.
1682}
1683
1685{
1686 /// This function fit's the final data of the given channel.
1687 /// It requires all other elements to have been created already.
1688 if(TSourceCalibration::VerboseLevel() >= EVerbosity::kSubroutines) { std::cout << __PRETTY_FUNCTION__ << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
1689 TF1* calibration = new TF1("fitfunction", ::Polynomial, 0., 10000., degree + 2);
1690 calibration->FixParameter(0, degree);
1692 std::cout << DYELLOW << "fData " << fData << ": ";
1693 if(fData != nullptr) {
1694 std::cout << fData->GetN() << " data points being fit, ";
1695 }
1696 double min = 0.;
1697 double max = 0.;
1698 calibration->GetRange(min, max);
1699 std::cout << "range " << min << " - " << max << std::endl;
1700 }
1701 fData->Fit(calibration, "Q");
1702 TString text = Form("%.6f + %.6f*x", calibration->GetParameter(1), calibration->GetParameter(2));
1703 for(int i = 2; i <= degree; ++i) {
1704 text.Append(Form(" + %.6f*x^%d", calibration->GetParameter(i + 1), i));
1705 }
1706 fChannelStatusBar->SetText(text.Data(), 2);
1707 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << "re-calculating residuals and clearing fields" << std::endl; }
1708 // re-calculate the residuals
1709 fData->SetResidual(true);
1710
1711 fLegend->Clear();
1712 fCalibrationPad->cd();
1714 fLegend->Draw();
1715 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << "set chi2 label" << std::endl; }
1716 // calculate the corners of the chi^2 label from the minimum and maximum x/y-values of the graph
1717 // we position it in the top left corner about 50% of the width and 10% of the height of the graph
1718 if(fChi2Label == nullptr) {
1719 double left = fData->GetMinimumX();
1720 double right = left + (fData->GetMaximumX() - left) * 0.5;
1721 double top = fData->GetMaximumY();
1722 double bottom = top - (top - fData->GetMinimumY()) * 0.1;
1723
1724 fChi2Label = new TPaveText(left, bottom, right, top);
1726 fChi2Label->ConvertNDCtoPad();
1727 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;
1728 fData->Print("e");
1729 }
1730 } else {
1731 fChi2Label->Clear();
1732 }
1733 fChi2Label->AddText(Form("#chi^{2}/NDF = %f", calibration->GetChisquare() / calibration->GetNDF()));
1734 fChi2Label->SetFillColor(kWhite);
1735 fChi2Label->Draw();
1737 std::cout << "fChi2Label set to:" << std::endl;
1738 fChi2Label->GetListOfLines()->Print();
1739 std::cout << "Text size " << fChi2Label->GetTextSize() << std::endl;
1740 }
1741 if(fInfoLabel == nullptr) {
1742 fInfoLabel = new TPaveText(0, 0.95, 1., 1., "NDCNB");
1744 fInfoLabel->ConvertNDCtoPad();
1745 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;
1746 }
1747 } else {
1748 fInfoLabel->Clear();
1750 fInfoLabel->ConvertNDCtoPad();
1751 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;
1752 }
1753 }
1754 fInfoLabel->AddText(Form("Using %d peaks from %.0f to %.0f keV", fData->TotalGraph()->GetN(), fData->GetMinimumY(), fData->GetMaximumY()));
1755 fInfoLabel->SetFillColor(kWhite);
1756 fInfoLabel->Draw();
1757
1758 fResidualPad->cd();
1759 fData->DrawResidual("*r0");
1760
1761 fCalibrationCanvas->GetCanvas()->Modified();
1762
1763 fData->FitFunction()->SetRange(0., 10000.);
1765 double min = 0.;
1766 double max = 0.;
1767 calibration->GetRange(min, max);
1768 std::cout << DYELLOW << "calibration has range " << min << " - " << max;
1769 if(fData->FitFunction() != nullptr) {
1770 fData->FitFunction()->GetRange(min, max);
1771 std::cout << ", fit function has range " << min << " - " << max;
1772 }
1773 std::cout << std::endl;
1774 }
1775 delete calibration;
1776
1777 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)
1778}
1779
1780void TChannelTab::RecursiveRemove(const double& maxResidual)
1781{
1782 /// This function goes through the residuals of the calibration and recursively removes each point that is above the maxResidual parameter.
1783 /// The function requires an initial calibration to have residuals available.
1784 /// After each point that is removed, it re-calibrates the data, producing new residuals.
1785
1786 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)
1787
1788 if(fData == nullptr) {
1789 std::cout << "No longer have a data for channel " << Name() << "?" << std::endl;
1790 return;
1791 }
1792
1793 if(fData->TotalResidualGraph() == nullptr) {
1794 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;
1795 return;
1796 }
1797
1798 if(fData->TotalResidualGraph()->GetN() == 0) {
1799 std::cout << "Recursively removed all data points for channel " << Name() << std::endl;
1800 return;
1801 }
1802
1803 // get the maximum residual
1804 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); });
1805 auto index = std::distance(fData->TotalResidualGraph()->GetX(), maxElement);
1806
1807 if(std::abs(*maxElement) < maxResidual) {
1808 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; }
1809 return;
1810 }
1811
1812 // remove the index with the maximum residual
1813 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << DYELLOW << "std::abs(" << *maxElement << ") > " << maxResidual << ": removing index " << index << "/" << fData->TotalResidualGraph()->GetN() << std::endl; }
1814 fData->RemovePoint(index);
1815
1816 // update the data, and re-calibrate
1817 UpdateData();
1818 UpdateFwhm();
1819 Calibrate();
1820
1821 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << "RecursiveRemove(" << maxResidual << ") being called again" << RESET_COLOR << std::endl; }
1822 RecursiveRemove(maxResidual);
1823}
1824
1826{
1827 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)
1828 // write the actual parameters of the fit
1829 std::stringstream str;
1830 str << std::hex << fName;
1831 int address = 0;
1832 str >> address;
1833 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << "Got address " << hex(address, 4) << " from name " << fName << std::endl; }
1834 TF1* calibration = fData->FitFunction();
1835 if(calibration == nullptr) {
1836 std::cout << "Failed to find calibration in fData" << std::endl;
1837 fData->TotalGraph()->GetListOfFunctions()->Print();
1838 return;
1839 }
1840 std::vector<Float_t> parameters;
1841 for(int i = 0; i <= calibration->GetParameter(0); ++i) {
1842 parameters.push_back(static_cast<Float_t>(calibration->GetParameter(i + 1)));
1843 }
1844 TChannel* channel = TChannel::GetChannel(address, false);
1845 if(channel == nullptr) {
1846 std::cerr << "Failed to get channel for address " << hex(address, 4) << std::endl;
1847 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)
1848 return;
1849 }
1850 channel->SetENGCoefficients(parameters);
1851 channel->DestroyEnergyNonlinearity();
1853 double* x = fData->GetX();
1854 double* y = fData->GetY();
1855 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << "Going to add " << fData->GetN() << " points to nonlinearity graph" << std::endl; }
1856 for(int i = 0; i < fData->GetN(); ++i) {
1857 // nonlinearity expects y to be the true source energy minus the calibrated energy of the peak
1858 // the source energy is y, the peak is x
1859 channel->AddEnergyNonlinearityPoint(y[i], y[i] - calibration->Eval(x[i]));
1860 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; }
1861 }
1862 }
1863}
1864
1865void TChannelTab::Initialize(const bool& force)
1866{
1867 // loop through sources and get rough calibration for each
1868 for(auto* source : fSources) {
1870 std::cout << DYELLOW << "Rough calibration in source tab:" << std::endl;
1871 source->Print();
1872 }
1873 source->InitialCalibration(force);
1874 fProgressBar->Increment(1);
1875 }
1876 // get this rough calibration data and calibrate with linear calibration
1877 UpdateData();
1878 Calibrate(1);
1879 // copy the data to the initial data and draw it
1880 fInit = static_cast<TCalibrationGraphSet*>(fData->Clone(Form("init%s", fName.c_str())));
1882 std::cout << "Cloned data " << fData->GetName() << " to " << fInit->GetName() << ":" << std::endl;
1883 fInit->Print();
1884 }
1885 auto* oldPad = gPad;
1886 fInitCanvas->GetCanvas()->cd();
1887 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << "Switched from " << oldPad->GetName() << " to " << gPad->GetName() << std::endl; }
1889 // calculate the corners of the label from the minimum and maximum x/y-values of the graph
1890 // we position it in the top left corner about 50% of the width and 10% of the height of the graph
1891 if(fCalLabel == nullptr) {
1892 double left = fInit->GetMinimumX();
1893 double right = left + (fInit->GetMaximumX() - left) * 0.5;
1894 double top = fInit->GetMaximumY();
1895 double bottom = top - (top - fInit->GetMinimumY()) * 0.1;
1896
1897 fCalLabel = new TPaveText(left, bottom, right, top);
1899 std::cout << "fCalLabel created " << fCalLabel << " (" << left << " - " << right << ", " << bottom << " - " << top << ", from " << fInit->GetMinimumX() << "-" << fInit->GetMaximumX() << ", " << fInit->GetMinimumY() << "-" << fInit->GetMaximumY() << ") on gPad " << gPad->GetName() << std::endl;
1900 }
1901 } else {
1902 fCalLabel->Clear();
1903 }
1904 std::stringstream str;
1905 if(fInit->FitFunction()->GetNpar() == 3) {
1906 str << "Fit function: " << fInit->FitFunction()->GetParameter(1) << " + " << fInit->FitFunction()->GetParameter(2) << " x";
1907 } else {
1908 str << "Unknown fit function with " << fInit->FitFunction()->GetNpar() << " parameters";
1909 }
1911 std::cout << "Created label text \"" << str.str() << "\"" << std::endl;
1912 }
1913
1914 fCalLabel->AddText(str.str().c_str());
1915 fCalLabel->SetFillColor(kWhite);
1916 fCalLabel->Draw();
1917 oldPad->cd();
1918 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << "Switched back to " << gPad->GetName() << std::endl; }
1919}
1920
1922{
1924 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;
1925 fSourceTab->Print();
1926 for(int tab = 0; tab < fSourceTab->GetNumberOfTabs(); ++tab) {
1927 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;
1928 }
1929 }
1930 fSources[fActiveSourceTab]->FindCalibratedPeaks(fData->FitFunction());
1931 UpdateData();
1932 UpdateFwhm();
1934 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;
1935 }
1936}
1937
1939{
1941 std::cout << DYELLOW << "Finding calibrated peaks in all source tabs" << std::endl;
1942 }
1943 for(auto* source : fSources) {
1945 std::cout << DYELLOW << "Finding calibrated peaks in source tab " << source->SourceName() << std::endl;
1946 }
1947 source->FindCalibratedPeaks(fData->FitFunction());
1948 }
1949 UpdateData();
1950 UpdateFwhm();
1952 std::cout << DYELLOW << "Done finding calibrated peaks in all source tabs" << std::endl;
1953 }
1954}
1955
1957{
1958 fFwhmCanvas->GetCanvas()->cd();
1959 fFwhm->DrawCalibration("*");
1960 fLegend->Draw();
1961 for(auto* source : fSources) {
1962 source->Draw();
1963 }
1964}
1965
1966void TChannelTab::Write(TFile* output)
1967{
1968 /// Write graphs to output.
1969 output->cd();
1970 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << DYELLOW << "writing " << fName << std::endl; }
1971 fData->Write(Form("calGraph_%s", fName.c_str()), TObject::kOverwrite);
1972 fFwhm->Write(Form("fwhm_%s", fName.c_str()), TObject::kOverwrite);
1973 if(TSourceCalibration::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << DYELLOW << "wrote data " << fName << std::endl; }
1974}
1975
1977{
1978 /// 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
1979 // find the histogram in the current pad
1980 TH1* hist1 = nullptr;
1981 for(const auto&& obj : *(gPad->GetListOfPrimitives())) {
1982 if(obj->InheritsFrom(TGraph::Class())) {
1983 hist1 = static_cast<TGraph*>(obj)->GetHistogram();
1984 break;
1985 }
1986 }
1987 if(hist1 == nullptr) {
1988 std::cout << "ZoomX - Failed to find histogram for pad " << gPad->GetName() << std::endl;
1989 if(TSourceCalibration::VerboseLevel() > EVerbosity::kLoops) { gPad->GetListOfPrimitives()->Print(); }
1990 return;
1991 }
1992
1993 // extract base name and channel from pad name
1994 std::string padName = gPad->GetName();
1995 std::string baseName = padName.substr(0, 3);
1996 std::string channelName = padName.substr(4);
1997
1998 // find the corresponding partner pad to the selected one
1999 std::string newName;
2000 if(baseName == "cal") {
2001 newName = "res_" + channelName;
2002 } else if(baseName == "res") {
2003 newName = "cal_" + channelName;
2004 } else {
2005 std::cout << "Unknown combination of pad " << gPad->GetName() << " and histogram " << hist1->GetName() << std::endl;
2006 return;
2007 }
2008 auto* newObj = gPad->GetCanvas()->FindObject(newName.c_str());
2009 if(newObj == nullptr) {
2010 std::cout << "Failed to find " << newName << std::endl;
2011 return;
2012 }
2013
2014 if(newObj->IsA() != TPad::Class()) {
2015 std::cout << newObj << " = " << newObj->GetName() << ", " << newObj->GetTitle() << " is a " << newObj->ClassName() << " not a TPad" << std::endl;
2016 return;
2017 }
2018 auto* pad2 = static_cast<TPad*>(newObj);
2019 if(pad2 == nullptr) {
2020 std::cout << "Failed to find partner pad " << newName << std::endl;
2021 return;
2022 }
2023 // find the histogram in the partner pad
2024 TH1* hist2 = nullptr;
2025 for(const auto&& obj : *(pad2->GetListOfPrimitives())) {
2026 if(obj->InheritsFrom(TGraph::Class())) {
2027 hist2 = static_cast<TGraph*>(obj)->GetHistogram();
2028 break;
2029 }
2030 }
2031 if(hist2 == nullptr) {
2032 std::cout << "ZoomX - Failed to find histogram for partner pad " << newName << std::endl;
2033 if(TSourceCalibration::VerboseLevel() > EVerbosity::kLoops) { pad2->GetListOfPrimitives()->Print(); }
2034 return;
2035 }
2036
2037 TAxis* axis1 = hist1->GetXaxis();
2038 Int_t binmin = axis1->GetFirst();
2039 Int_t binmax = axis1->GetLast();
2040 Double_t xmin = axis1->GetBinLowEdge(binmin);
2041 Double_t xmax = axis1->GetBinLowEdge(binmax);
2042 TAxis* axis2 = hist2->GetXaxis();
2043 Int_t newmin = axis2->FindBin(xmin);
2044 Int_t newmax = axis2->FindBin(xmax);
2045 axis2->SetRange(newmin, newmax);
2046 pad2->Modified();
2047 pad2->Update();
2048}
2049
2051{
2052 /// 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
2053 // find the histogram in the current pad
2054 TH1* hist1 = nullptr;
2055 for(const auto&& obj : *(gPad->GetListOfPrimitives())) {
2056 if(obj->InheritsFrom(TGraph::Class())) {
2057 hist1 = static_cast<TGraph*>(obj)->GetHistogram();
2058 break;
2059 }
2060 }
2061 if(hist1 == nullptr) {
2062 std::cout << "ZoomY - Failed to find histogram for pad " << gPad->GetName() << std::endl;
2063 if(TSourceCalibration::VerboseLevel() > EVerbosity::kLoops) { gPad->GetListOfPrimitives()->Print(); }
2064 return;
2065 }
2066
2067 // extract base name and channel from pad name
2068 std::string padName = gPad->GetName();
2069 std::string baseName = padName.substr(0, 3);
2070 std::string channelName = padName.substr(4);
2071
2072 // find the corresponding partner pad to the selected one
2073 std::string newName;
2074 if(baseName == "cal") {
2075 newName = "res_" + channelName;
2076 } else if(baseName == "res") {
2077 newName = "cal_" + channelName;
2078 } else {
2079 std::cout << "Unknown combination of pad " << gPad->GetName() << " and histogram " << hist1->GetName() << std::endl;
2080 return;
2081 }
2082 auto* newObj = gPad->GetCanvas()->FindObject(newName.c_str());
2083 if(newObj == nullptr) {
2084 std::cout << "Failed to find " << newName << std::endl;
2085 return;
2086 }
2087
2088 if(newObj->IsA() != TPad::Class()) {
2089 std::cout << newObj << " = " << newObj->GetName() << ", " << newObj->GetTitle() << " is a " << newObj->ClassName() << " not a TPad" << std::endl;
2090 return;
2091 }
2092 auto* pad2 = static_cast<TPad*>(newObj);
2093 if(pad2 == nullptr) {
2094 std::cout << "Failed to find partner pad " << newName << std::endl;
2095 return;
2096 }
2097 // find the histogram in the partner pad
2098 TH1* hist2 = nullptr;
2099 for(const auto&& obj : *(pad2->GetListOfPrimitives())) {
2100 if(obj->InheritsFrom(TGraph::Class())) {
2101 hist2 = static_cast<TGraph*>(obj)->GetHistogram();
2102 if(obj->IsA() == TCalibrationGraphSet::Class()) {
2103 static_cast<TCalibrationGraphSet*>(obj)->RefreshResidualLine();
2104 }
2105 break;
2106 }
2107 }
2108 if(hist2 == nullptr) {
2109 std::cout << "ZoomY - Failed to find histogram for partner pad " << newName << std::endl;
2110 if(TSourceCalibration::VerboseLevel() > EVerbosity::kLoops) { pad2->GetListOfPrimitives()->Print(); }
2111 return;
2112 }
2113
2114 hist2->SetMinimum(hist1->GetMinimum());
2115 hist2->SetMaximum(hist1->GetMaximum());
2116
2117 pad2->Modified();
2118 pad2->Update();
2119}
2120
2122{
2123 std::cout << "TChannelTab " << Name() << " calibration canvas:" << std::endl;
2124 fCalibrationCanvas->GetCanvas()->Print();
2125 for(auto* obj : *(fCalibrationCanvas->GetCanvas()->GetListOfPrimitives())) {
2126 obj->Print();
2127 }
2128 std::cout << "TChannelTab fwhm canvas:" << std::endl;
2129 fFwhmCanvas->GetCanvas()->Print();
2130 for(auto* obj : *(fFwhmCanvas->GetCanvas()->GetListOfPrimitives())) {
2131 obj->Print();
2132 }
2133 for(auto* sourceTab : fSources) {
2134 sourceTab->PrintCanvases();
2135 }
2136}
2137
2139{
2140 std::cout << "TChannelTab frame:" << std::endl;
2141 fChannelFrame->Print();
2142 std::cout << "TChannelTab canvas:" << std::endl;
2143 fCalibrationCanvas->Print();
2144 std::cout << "TChannelTab status bar:" << std::endl;
2145 //fChannelStatusBar->Print();
2146 for(auto* sourceTab : fSources) {
2147 sourceTab->PrintLayout();
2148 }
2149}
2150
2151//////////////////////////////////////// TSourceCalibration ////////////////////////////////////////
2163bool TSourceCalibration::fFast = false;
2167std::vector<int> TSourceCalibration::fBadBins;
2168
2169TSourceCalibration::TSourceCalibration(double sigma, double threshold, int degree, int count...)
2170 : TGMainFrame(nullptr, 2 * fPanelWidth, fPanelHeight + 2 * fStatusbarHeight), fDefaultSigma(sigma), fDefaultThreshold(threshold), fDefaultDegree(degree)
2171{
2172 TH1::AddDirectory(false); // turns off warnings about multiple histograms with the same name because ROOT doesn't manage them anymore
2173
2174 va_list args;
2175 va_start(args, count); // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2176 for(int i = 0; i < count; ++i) {
2177 fMatrices.push_back(va_arg(args, TH2*)); // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2178 if(fMatrices.back() == nullptr) {
2179 std::cout << "Error, got nullptr as matrix input?" << std::endl;
2180 fMatrices.pop_back();
2181 }
2182 }
2183 va_end(args); // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2185 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)
2186 << "using " << count << "/" << fMatrices.size() << " matrices:" << std::endl;
2187 for(auto* mat : fMatrices) {
2188 std::cout << mat << std::flush << " = " << mat->GetName() << std::endl;
2189 }
2190 }
2191
2192 fOldErrorLevel = gErrorIgnoreLevel;
2193 gErrorIgnoreLevel = kError;
2194 gPrintViaErrorHandler = true; // redirects all printf's to root's normal message system
2195
2196 if(!fLogFile.empty()) {
2197 fRedirect = new TRedirect(fLogFile.c_str());
2198 std::cout << "Starting redirect to " << fLogFile << std::endl;
2199 }
2200
2201 // check matrices (# of filled bins and bin labels) and resize some vectors for later use
2202 // use the first matrix to get a reference for everything
2203 fNofBins = 0;
2204 std::map<int, int> channelToIndex; // used to be a member, but only used here
2205 for(int bin = 1; bin <= fMatrices[0]->GetNbinsX(); ++bin) {
2206 if(FilledBin(fMatrices[0], bin) && std::find(fBadBins.begin(), fBadBins.end(), bin) == fBadBins.end()) {
2207 fActualChannelId.push_back(fNofBins); // at this point fNofBins is the index at which this projection will end up
2208 fActiveBins.push_back(bin);
2209 channelToIndex[bin] = fNofBins; // this map simply stores which bin ends up at which index
2210 fChannelLabel.push_back(fMatrices[0]->GetXaxis()->GetBinLabel(bin));
2211 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << bin << ". bin: fNofBins " << fNofBins << ", channelToIndex[" << bin << "] " << channelToIndex[bin] << ", fActualChannelId.back() " << fActualChannelId.back() << std::endl; }
2212 ++fNofBins;
2213 } else {
2214 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << "skipping bin " << bin << std::endl; }
2215 }
2216 }
2217 // now loop over all other matrices and check vs. the first one
2218 for(size_t i = 1; i < fMatrices.size(); ++i) {
2219 int tmpBins = 0;
2220 for(int bin = 1; bin <= fMatrices[i]->GetNbinsX(); ++bin) {
2221 if(FilledBin(fMatrices[i], bin) && std::find(fBadBins.begin(), fBadBins.end(), bin) == fBadBins.end()) { // good bin in the current matrix
2222 // current index is tmpBins, so we check if the bin agrees with what we have
2223 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
2224 std::ostringstream error;
2225 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;
2226 throw std::invalid_argument(error.str());
2227 }
2228 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
2229 std::ostringstream error;
2230 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;
2231 throw std::invalid_argument(error.str());
2232 }
2233 ++tmpBins;
2234 } else {
2235 if(channelToIndex.find(bin) != channelToIndex.end()) { //found an empty bin that was full in the first matrix
2236 std::ostringstream error;
2237 error << "Mismatch in " << i << ". matrix (" << fMatrices[i]->GetName() << "), bin " << bin << " is empty, but should be " << channelToIndex[bin] << ". filled bin" << std::endl;
2238 throw std::invalid_argument(error.str());
2239 }
2240 }
2241 }
2242 }
2243
2244 if(fVerboseLevel > EVerbosity::kQuiet) { std::cout << fMatrices.size() << " matrices with " << fMatrices[0]->GetNbinsX() << " x-bins, fNofBins " << fNofBins << ", fActualChannelId.size() " << fActualChannelId.size() << std::endl; }
2245
2246 fOutput = new TFile("TSourceCalibration.root", "recreate");
2247 if(!fOutput->IsOpen()) {
2248 throw std::runtime_error("Unable to open output file \"TSourceCalibration.root\"!");
2249 }
2250
2251 // build the first screen
2254
2255 // Set a name to the main frame
2256 SetWindowName("SourceCalibration");
2257
2258 // Map all subwindows of main frame
2259 MapSubwindows();
2260
2261 // Initialize the layout algorithm
2262 Resize(GetDefaultSize());
2263
2264 // Map main frame
2265 MapWindow();
2266}
2267
2269{
2270 fOutput->Close();
2271 delete fOutput;
2272
2273 gErrorIgnoreLevel = fOldErrorLevel;
2274
2275 delete fRedirect;
2276}
2277
2279{
2280 /// Build initial interface with histogram <-> source assignment
2281
2282 auto* layoutManager = new TGTableLayout(this, fMatrices.size() + 1, 2, true, 5);
2283 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << "created table layout manager with 2 columns, " << fMatrices.size() + 1 << " rows" << std::endl; }
2284 SetLayoutManager(layoutManager);
2285
2286 // find the path for the .sou files
2287 std::string path = TNucleus::SourceDirectory();
2288 // The matrices and source selection boxes
2289 size_t i = 0;
2290 for(i = 0; i < fMatrices.size(); ++i) {
2291 fMatrixNames.push_back(new TGLabel(this, fMatrices[i]->GetName()));
2292 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "Text height " << fMatrixNames.back()->GetFont()->TextHeight() << std::endl; }
2293 fSource.push_back(nullptr);
2294 fSourceEnergy.emplace_back();
2295 fSourceBox.push_back(new TGComboBox(this, "Select source", kSourceBox + fSourceBox.size()));
2296
2297 int index = 0;
2298 // convert name to lower-case
2299 std::string name = fMatrices[i]->GetName();
2300 std::transform(name.begin(), name.end(), name.begin(), ::tolower);
2301 fMatrices[i]->SetName(name.c_str());
2302#if __cplusplus >= 201703L
2303 for(const auto& file : std::filesystem::directory_iterator(path)) {
2304 if(file.is_regular_file() && file.path().extension().compare(".sou") == 0) {
2305 fSourceBox.back()->AddEntry(file.path().stem().c_str(), index);
2306 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << i << "/" << index << ": Comparing matrix name " << fMatrices[i]->GetName() << " to file name " << file.path().stem().c_str() << std::endl; }
2307 if(std::strstr(fMatrices[i]->GetName(), file.path().stem().c_str()) != nullptr) {
2308 fSourceBox.back()->Select(index);
2309 SetSource(kSourceBox + fSourceBox.size() - 1, index);
2310 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << i << ": selected source " << index << std::endl; }
2312 std::cout << "matrix name " << fMatrices[i]->GetName() << " not matching " << file.path().stem().c_str() << std::endl;
2313 }
2314 ++index;
2315 }
2316 }
2317 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << i << ": created list with " << index << " sources" << std::endl; }
2318#else
2319 fSourceBox.back()->AddEntry("22Na", index);
2320 if(std::strstr(fMatrices[i]->GetName(), "22na") != nullptr) {
2321 fSourceBox.back()->Select(index);
2322 SetSource(kSourceBox + fSourceBox.size() - 1, index);
2323 }
2324 ++index;
2325 fSourceBox.back()->AddEntry("56Co", index);
2326 if(std::strstr(fMatrices[i]->GetName(), "56co") != nullptr) {
2327 fSourceBox.back()->Select(index);
2328 SetSource(kSourceBox + fSourceBox.size() - 1, index);
2329 }
2330 ++index;
2331 fSourceBox.back()->AddEntry("60Co", index);
2332 if(std::strstr(fMatrices[i]->GetName(), "60co") != nullptr) {
2333 fSourceBox.back()->Select(index);
2334 SetSource(kSourceBox + fSourceBox.size() - 1, index);
2335 }
2336 ++index;
2337 fSourceBox.back()->AddEntry("133Ba", index);
2338 if(std::strstr(fMatrices[i]->GetName(), "133ba") != nullptr) {
2339 fSourceBox.back()->Select(index);
2340 SetSource(kSourceBox + fSourceBox.size() - 1, index);
2341 }
2342 ++index;
2343 fSourceBox.back()->AddEntry("152Eu", index);
2344 if(std::strstr(fMatrices[i]->GetName(), "152eu") != nullptr) {
2345 fSourceBox.back()->Select(index);
2346 SetSource(kSourceBox + fSourceBox.size() - 1, index);
2347 }
2348 ++index;
2349 fSourceBox.back()->AddEntry("241Am", index);
2350 if(std::strstr(fMatrices[i]->GetName(), "241am") != nullptr) {
2351 fSourceBox.back()->Select(index);
2352 SetSource(kSourceBox + fSourceBox.size() - 1, index);
2353 }
2354#endif
2355
2356 fSourceBox.back()->SetMinHeight(fPanelHeight / 2);
2357
2358 fSourceBox.back()->Resize(fSourceboxWidth, fLineHeight);
2359 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; }
2360 AddFrame(fMatrixNames.back(), new TGTableLayoutHints(0, 1, i, i + 1, kLHintsRight | kLHintsCenterY));
2361 AddFrame(fSourceBox.back(), new TGTableLayoutHints(1, 2, i, i + 1, kLHintsLeft | kLHintsCenterY));
2362 }
2363
2364 // The buttons
2365 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << "Attaching start button to 0, 2, " << i << ", " << i + 1 << std::endl; }
2366 fStartButton = new TGTextButton(this, "Accept && Continue", kStartButton);
2367 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << "start button " << fStartButton << std::endl; }
2368 AddFrame(fStartButton, new TGTableLayoutHints(0, 2, i, i + 1, kLHintsCenterX | kLHintsCenterY));
2369 Layout();
2370}
2371
2373{
2374 /// Create connections for the histogram <-> source assignment interface
2375
2376 // Connect the selection of the source
2377 for(auto* box : fSourceBox) {
2378 box->Connect("Selected(Int_t, Int_t)", "TSourceCalibration", this, "SetSource(Int_t, Int_t)");
2379 }
2380
2381 //Connect the clicking of buttons
2382 fStartButton->Connect("Clicked()", "TSourceCalibration", this, "Start()");
2383}
2384
2386{
2387 /// Disconnect all signals from histogram <-> source assignment interface
2388 for(auto* box : fSourceBox) {
2389 box->Disconnect("Selected(Int_t, Int_t)", this, "SetSource(Int_t, Int_t)");
2390 }
2391
2392 fStartButton->Disconnect("Clicked()", this, "Start()");
2393}
2394
2396{
2397 HideFrame(element);
2398 RemoveFrame(element);
2399 delete element;
2400 element = nullptr;
2401}
2402
2404{
2405 for(auto* name : fMatrixNames) {
2406 DeleteElement(name);
2407 }
2408 fMatrixNames.clear();
2409 for(auto* box : fSourceBox) {
2410 DeleteElement(box);
2411 }
2412 fSourceBox.clear();
2414 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "Deleted start button " << fStartButton << std::endl; }
2415}
2416
2417void TSourceCalibration::SetSource(Int_t windowId, Int_t entryId)
2418{
2419 int index = windowId - kSourceBox;
2420 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)
2421 delete fSource[index];
2422 auto* nucleus = new TNucleus(fSourceBox[index]->GetListBox()->GetEntry(entryId)->GetTitle());
2423 TIter iter(nucleus->GetTransitionList());
2424 fSourceEnergy[index].clear();
2425 while(auto* transition = static_cast<TTransition*>(iter.Next())) {
2426 fSourceEnergy[index].emplace_back(transition->GetEnergy(), transition->GetEnergyUncertainty(), transition->GetIntensity(), transition->GetIntensityUncertainty());
2427 }
2428 fSource[index] = nucleus;
2429}
2430
2432{
2433 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)
2434 if(fEmitter == fStartButton) {
2435 SecondWindow();
2436 } else if(fEmitter == fAcceptAllButton) {
2437 for(auto& channelTab : fChannelTab) {
2438 if(channelTab == nullptr) { continue; }
2439 channelTab->UpdateChannel();
2440 channelTab->Write(fOutput);
2441 }
2443 std::cout << "Closing window" << std::endl;
2444 CloseWindow();
2445 std::cout << "all done" << std::endl;
2446 exit(0);
2447 }
2448}
2449
2451{
2452 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)
2453 if(fEmitter == nullptr) { // we only want to do this once at the beginning (after fEmitter was initialized to nullptr)
2455 TTimer::SingleShot(fWaitMs, "TSourceCalibration", this, "HandleTimer()");
2456 }
2457}
2458
2460{
2461 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << __PRETTY_FUNCTION__ << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2462 // check that all sources have been set
2463 for(size_t i = 0; i < fSource.size(); ++i) {
2464 if(fSource[i] == nullptr) {
2465 std::cerr << "Source " << i << " not set!" << std::endl;
2466 return;
2467 }
2468 }
2469 // now check that we don't have the same source twice (which wouldn't make sense)
2470 for(size_t i = 0; i < fSource.size(); ++i) {
2471 for(size_t j = i + 1; j < fSource.size(); ++j) {
2472 if(*(fSource[i]) == *(fSource[j])) {
2473 std::cerr << "Duplicate sources: " << i << " - " << fSource[i]->GetName() << " and " << j << " - " << fSource[j]->GetName() << std::endl;
2474 return;
2475 }
2476 }
2477 }
2478
2479 // disconnect signals of first screen and remove all elements
2481 RemoveAll();
2482 DeleteFirst();
2483
2484 SetLayoutManager(new TGVerticalLayout(this));
2485
2486 // create intermediate progress bar
2487 fProgressBar = new TGHProgressBar(this, TGProgressBar::kFancy, fPanelWidth);
2488 fProgressBar->SetBarColor("lightblue");
2489 fProgressBar->SetRange(0., static_cast<Float_t>(fMatrices.size() * fActiveBins.size()));
2490 fProgressBar->ShowPosition(true, false, "Creating source tabs: %.0f");
2491 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; }
2492 AddFrame(fProgressBar, new TGLayoutHints(kLHintsCenterX | kLHintsCenterY | kLHintsExpandX | kLHintsExpandY, 0, 0, 0, 0));
2493
2494 // Map all subwindows of main frame
2495 MapSubwindows();
2496
2497 // Initialize the layout algorithm
2498 Resize(TGDimension(fPanelWidth, fLineHeight));
2499
2500 // Map main frame
2501 MapWindow();
2502
2503 // create second screen and its connections
2506
2507 // remove progress bar
2509
2510 // Map all subwindows of main frame
2511 MapSubwindows();
2512
2513 // Initialize the layout algorithm
2514 Resize(TGDimension(2 * fPanelWidth, fPanelHeight + 2 * fStatusbarHeight));
2515
2516 // Map main frame
2517 MapWindow();
2518 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)
2519}
2520
2522{
2523 SetLayoutManager(new TGVerticalLayout(this));
2524 fTab = new TGTab(this, 2 * fPanelWidth, fPanelHeight + 2 * fStatusbarHeight);
2525
2526 fChannelTab.resize(fMatrices[0]->GetNbinsX());
2527 for(auto& bin : fActiveBins) {
2528 // create vector with projections of this channel for each source
2529 std::vector<GH1D*> projections(fMatrices.size());
2530 size_t index = 0;
2531 for(auto* matrix : fMatrices) {
2532 projections[index] = new GH1D(matrix->ProjectionY(Form("%s_%s", fSource[index]->GetName(), matrix->GetXaxis()->GetBinLabel(bin)), bin, bin));
2533 ++index;
2534 }
2535 auto* frame = fTab->AddTab(fMatrices[0]->GetXaxis()->GetBinLabel(bin));
2536 // create a new tab in a new thread
2537 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << DGREEN << "Creating channel for bin " << bin << std::endl; }
2538 fChannelTab[bin - 1] = new TChannelTab(this, fSource, projections, fMatrices[0]->GetXaxis()->GetBinLabel(bin), frame, fSourceEnergy, fProgressBar);
2539 //// right now the lambda uses copies of all variables, not sure if that is neccesary, or if we only need that for projections?
2540 //fFutures.emplace(std::make_pair(bin - 1, std::async(std::launch::async, [=] {
2541 // 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 << ", " << fDefaultThreshold << ", " << fDefaultDegree << ", " << fSourceEnergy.size() << " source energy vectors, " << fProgressBar << "."<< std::endl;
2542 // auto* newTab = new TChannelTab(this, fSource, projections, fMatrices[0]->GetXaxis()->GetBinLabel(bin), frame, fDefaultSigma, fDefaultThreshold, fDefaultDegree, fSourceEnergy, fProgressBar);
2543 // std::cout << "created new tab " << newTab << std::endl;
2544 // return newTab;
2545 // })));
2546 //// 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
2547 //if(fFutures.size() >= fNumberOfThreads) {
2548 // 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; }
2549 // auto tmp = fFutures.front().second.get();
2550 // fChannelTab[fFutures.front().first] = tmp;
2551 // fFutures.pop();
2552 //}
2553 }
2554 // wait for all remaining threads to finish
2555 //while(!fFutures.empty()) {
2556 // 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; }
2557 // auto tmp = fFutures.front().second.get();
2558 // fChannelTab[fFutures.front().first] = tmp;
2559 // fFutures.pop();
2560 //}
2561
2562 //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; }
2563
2564 //for(int i = 0; i < fChannelTab[0]->SourceTab()->GetNumberOfTabs(); ++i) {
2565 //if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << i << ": " << fChannelTab[0]->SourceTab()->GetTabTab(i)->GetText()->GetString() << std::endl; }
2566 //}
2567
2568 AddFrame(fTab, new TGLayoutHints(kLHintsTop | kLHintsExpandX | kLHintsExpandY, 0, 0, 0, 0));
2569
2570 // bottom frame with navigation button group, text entries, etc.
2571 fBottomFrame = new TGHorizontalFrame(this, 2 * fPanelWidth, TSourceCalibration::StatusbarHeight());
2572
2573 fLeftFrame = new TGVerticalFrame(fBottomFrame, fPanelWidth, fParameterHeight);
2574 fNavigationGroup = new TGHButtonGroup(fLeftFrame, "Navigation");
2575 fPreviousButton = new TGTextButton(fNavigationGroup, "&Previous");
2576 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << DGREEN << "prev button " << fPreviousButton << std::endl; }
2577 fPreviousButton->SetEnabled(false);
2578 fDiscardButton = new TGTextButton(fNavigationGroup, "&Discard");
2579 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "discard button " << fDiscardButton << std::endl; }
2580 fAcceptButton = new TGTextButton(fNavigationGroup, "&Accept");
2581 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "accept button " << fAcceptButton << std::endl; }
2582 fAcceptAllButton = new TGTextButton(fNavigationGroup, "Accept All && Finish");
2583 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "accept all button " << fAcceptAllButton << std::endl; }
2584 fWriteButton = new TGTextButton(fNavigationGroup, "&Write");
2585 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "write button " << fWriteButton << std::endl; }
2586 fNextButton = new TGTextButton(fNavigationGroup, "&Next");
2587 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "next button " << fNextButton << std::endl; }
2588
2589 fLeftFrame->AddFrame(fNavigationGroup, new TGLayoutHints(kLHintsTop | kLHintsExpandX, 2, 2, 2, 2));
2590
2591 fFittingGroup = new TGHButtonGroup(fLeftFrame, "Fitting/Calibrating");
2592 fFindSourcePeaksButton = new TGTextButton(fFittingGroup, "&Find Source Peaks");
2593 fFindChannelPeaksButton = new TGTextButton(fFittingGroup, "Find Channel Peaks");
2594 fFindAllPeaksButton = new TGTextButton(fFittingGroup, "Find All Peaks");
2595 fCalibrateButton = new TGTextButton(fFittingGroup, "&Calibrate");
2596 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "cal button " << fCalibrateButton << std::endl; }
2597 fRecursiveRemoveButton = new TGTextButton(fFittingGroup, "Recursive Remove");
2598 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "remove button " << fRecursiveRemoveButton << std::endl; }
2599
2600 fLeftFrame->AddFrame(fFittingGroup, new TGLayoutHints(kLHintsBottom | kLHintsExpandX, 2, 2, 2, 2));
2601
2602 fBottomFrame->AddFrame(fLeftFrame, new TGLayoutHints(kLHintsLeft, 2, 2, 2, 2));
2603
2604 fRightFrame = new TGVerticalFrame(fBottomFrame, fPanelWidth, fParameterHeight);
2605 fParameterFrame = new TGGroupFrame(fRightFrame, "Parameters", kHorizontalFrame);
2606 fParameterFrame->SetLayoutManager(new TGMatrixLayout(fParameterFrame, 0, 4, 2));
2607 fSigmaLabel = new TGLabel(fParameterFrame, "Sigma");
2608 fSigmaEntry = new TGNumberEntry(fParameterFrame, fDefaultSigma, fDigitWidth, kSigmaEntry, TGNumberFormat::EStyle::kNESRealTwo, TGNumberFormat::EAttribute::kNEAPositive);
2609 fThresholdLabel = new TGLabel(fParameterFrame, "Threshold");
2610 fThresholdEntry = new TGNumberEntry(fParameterFrame, fDefaultThreshold, fDigitWidth, kThresholdEntry, TGNumberFormat::EStyle::kNESRealThree, TGNumberFormat::EAttribute::kNEAPositive, TGNumberFormat::ELimit::kNELLimitMinMax, 0., 1.);
2611 fDegreeLabel = new TGLabel(fParameterFrame, "Degree of polynomial");
2612 fDegreeEntry = new TGNumberEntry(fParameterFrame, fDefaultDegree, 2, kDegreeEntry, TGNumberFormat::EStyle::kNESInteger, TGNumberFormat::EAttribute::kNEAPositive);
2613 fMaxResidualLabel = new TGLabel(fParameterFrame, "Max. residual (in keV)");
2614 fMaxResidualEntry = new TGNumberEntry(fParameterFrame, fDefaultMaxResidual, fDigitWidth, kMaxResidualEntry, TGNumberFormat::EStyle::kNESRealTwo, TGNumberFormat::EAttribute::kNEAPositive);
2615
2616 fParameterFrame->AddFrame(fSigmaLabel);
2617 fParameterFrame->AddFrame(fSigmaEntry);
2620 fParameterFrame->AddFrame(fDegreeLabel);
2621 fParameterFrame->AddFrame(fDegreeEntry);
2624
2625 fRightFrame->AddFrame(fParameterFrame, new TGLayoutHints(kLHintsTop | kLHintsRight | kLHintsExpandX, 2, 2, 2, 2));
2626
2627 fSettingsFrame = new TGGroupFrame(fRightFrame, "Settings", kHorizontalFrame);
2628 fSettingsFrame->SetLayoutManager(new TGMatrixLayout(fSettingsFrame, 0, 2, 2));
2629 fApplyToAll = new TGCheckButton(fSettingsFrame, "Apply to all channels", kApplyToAll);
2630 fApplyToAll->SetState(kButtonUp);
2631 fSettingsFrame->AddFrame(fApplyToAll);
2632 fWriteNonlinearities = new TGCheckButton(fSettingsFrame, "Write nonlinearities", kWriteNonlinearities);
2633 fWriteNonlinearities->SetState(kButtonUp);
2635
2636 fRightFrame->AddFrame(fSettingsFrame, new TGLayoutHints(kLHintsBottom | kLHintsRight | kLHintsExpandX, 2, 2, 2, 2));
2637
2638 fBottomFrame->AddFrame(fRightFrame, new TGLayoutHints(kLHintsBottom | kLHintsRight, 2, 2, 2, 2));
2639
2640 AddFrame(fBottomFrame, new TGLayoutHints(kLHintsBottom | kLHintsExpandX, 2, 2, 2, 2));
2641
2642 SelectedTab(0);
2643
2644 fProgressBar->Reset();
2645 fProgressBar->ShowPosition(true, false, "Finding peaks for tab #%.0f");
2646 for(auto* channelTab : fChannelTab) {
2647 if(channelTab != nullptr) {
2649 for(int i = 0; i < 20; ++i) { std::cout << std::endl; }
2650 for(int i = 0; i < 5; ++i) {
2651 std::cout << "================================================================================" << std::endl;
2652 }
2653 for(int i = 0; i < 20; ++i) { std::cout << std::endl; }
2654 }
2655 channelTab->Initialize(true);
2656 channelTab->FindAllCalibratedPeaks();
2657 channelTab->Calibrate();
2658 channelTab->Draw();
2660 channelTab->PrintCanvases();
2661 }
2662 }
2663 }
2664
2665 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << DGREEN << "Second interface done" << std::endl; }
2666}
2667
2669{
2670 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << DGREEN << __PRETTY_FUNCTION__ << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2671 fNavigationGroup->Connect("Clicked(Int_t)", "TSourceCalibration", this, "Navigate(Int_t)");
2672 fFittingGroup->Connect("Clicked(Int_t)", "TSourceCalibration", this, "Fitting(Int_t)");
2673 fTab->Connect("Selected(Int_t)", "TSourceCalibration", this, "SelectedTab(Int_t)");
2674 // we don't need to connect the sigma, threshold, and degree number entries, those are automatically read when we start the calibration
2675 for(auto* channelTab : fChannelTab) {
2676 if(channelTab != nullptr) {
2677 channelTab->MakeConnections();
2678 }
2679 }
2680}
2681
2683{
2684 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << DGREEN << __PRETTY_FUNCTION__ << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2685 fNavigationGroup->Disconnect("Clicked(Int_t)", this, "Navigate(Int_t)");
2686 fFittingGroup->Disconnect("Clicked(Int_t)", this, "Fitting(Int_t)");
2687 fTab->Disconnect("Selected(Int_t)", this, "SelectedTab(Int_t)");
2688 for(auto* channelTab : fChannelTab) {
2689 if(channelTab != nullptr) {
2690 channelTab->Disconnect();
2691 }
2692 }
2693}
2694
2696{
2697 // we get the current source tab id and use it to get the channel tab from the right source tab
2698 // 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
2699 // 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
2700 int currentChannelId = fTab->GetCurrent();
2701 int actualChannelId = fActualChannelId[currentChannelId];
2702 int nofTabs = fTab->GetNumberOfTabs();
2703 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)
2704 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; }
2705 switch(id) {
2706 case ENavigate::kPrevious: // previous
2707 fTab->SetTab(currentChannelId - 1);
2708 SelectedTab(currentChannelId - 1);
2709 break;
2710 case ENavigate::kDiscard: // discard
2711 // select the next (or if we are on the last tab, the previous) tab
2712 if(currentChannelId < nofTabs - 1) {
2713 fTab->SetTab(currentChannelId + 1);
2714 } else {
2715 fTab->SetTab(currentChannelId - 1);
2716 }
2717 // remove the original active tab
2718 fTab->RemoveTab(currentChannelId);
2719 break;
2720 case ENavigate::kAccept: // accept
2721 AcceptChannel(currentChannelId);
2722 break;
2723 case ENavigate::kAcceptAll: // accept all (no argument = -1 = all)
2724 AcceptChannel();
2725 return;
2726 case ENavigate::kWrite: // write
2728 break;
2729 case ENavigate::kNext: // next
2730 fTab->SetTab(currentChannelId + 1);
2731 SelectedTab(currentChannelId + 1);
2732 break;
2733 default:
2734 break;
2735 }
2736 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; }
2737}
2738
2740{
2741 // we get the current source tab id and use it to get the channel tab from the right source tab
2742 // 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
2743 // 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
2744 int currentChannelId = fTab->GetCurrent();
2745 int actualChannelId = fActualChannelId[currentChannelId];
2746 int nofTabs = fTab->GetNumberOfTabs();
2747 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)
2748 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; }
2749 switch(id) {
2750 case EFitting::kFindSourcePeaks: // find peaks for selected source
2752 break;
2753 case EFitting::kFindChannelPeaks: // find peaks for selected channel
2755 break;
2756 case EFitting::kFindAllPeaks: // find all peaks
2757 FindAllPeaks();
2758 break;
2759 case EFitting::kCalibrate: // calibrate
2760 Calibrate();
2761 break;
2762 case EFitting::kRecursiveRemove: // recursively remove outliers
2764 break;
2765 default:
2766 break;
2767 }
2768 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; }
2769}
2770
2772{
2773 /// Simple function that enables and disables the previous and next buttons depending on which tab was selected.
2774 /// Also sets gPad to the projection of the source selected in this tab.
2775 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)
2776 if(id < 0) { return; }
2777 if(id == 0) {
2778 fPreviousButton->SetEnabled(false);
2779 } else {
2780 fPreviousButton->SetEnabled(true);
2781 }
2782
2783 if(id >= fTab->GetNumberOfTabs() - 1) {
2784 fNextButton->SetEnabled(false);
2785 } else {
2786 fNextButton->SetEnabled(true);
2787 }
2788
2789 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; }
2790
2791 if(fChannelTab[fActiveBins[id] - 1] == nullptr) {
2792 std::cout << "Failed to get fChannelTab[" << fActiveBins[id] - 1 << " = fActiveBins[" << id << "]-1]" << RESET_COLOR << std::endl;
2793 return;
2794 }
2795 if(fChannelTab[fActiveBins[id] - 1]->SelectedSourceTab() == nullptr) {
2796 std::cout << "Failed to get fChannelTab[" << fActiveBins[id] - 1 << " = fActiveBins[" << id << "]-1]->SelectedSourceTab()" << RESET_COLOR << std::endl;
2797 return;
2798 }
2799 if(fChannelTab[fActiveBins[id] - 1]->SelectedSourceTab()->ProjectionCanvas() == nullptr) {
2800 std::cout << "Failed to get fChannelTab[" << fActiveBins[id] - 1 << " = fActiveBins[" << id << "]-1]->SelectedSourceTab()->ProjectionCanvas()" << RESET_COLOR << std::endl;
2801 return;
2802 }
2803 if(fChannelTab[fActiveBins[id] - 1]->SelectedSourceTab()->ProjectionCanvas()->GetCanvas() == nullptr) {
2804 std::cout << "Failed to get fChannelTab[" << fActiveBins[id] - 1 << " = fActiveBins[" << id << "]-1]->SelectedSourceTab()->ProjectionCanvas()->GetCanvas()" << RESET_COLOR << std::endl;
2805 return;
2806 }
2807 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << DGREEN << "Changing gpad from \"" << gPad->GetName() << "\" to \""; }
2808 gPad = fChannelTab[fActiveBins[id] - 1]->SelectedSourceTab()->ProjectionCanvas()->GetCanvas();
2809 fChannelTab[fActiveBins[id] - 1]->SelectedSourceTab()->ProjectionCanvas()->GetCanvas()->cd();
2810 if(TSourceCalibration::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << gPad->GetName() << "\"" << RESET_COLOR << std::endl; }
2811}
2812
2813void TSourceCalibration::AcceptChannel(const int& channelId)
2814{
2815 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)
2816
2817 // select the next (or if we are on the last tab, the previous) tab
2818 int nofTabs = fTab->GetNumberOfTabs();
2819 int minChannel = 0;
2820 int maxChannel = nofTabs - 1;
2821 if(channelId >= 0) {
2822 minChannel = channelId;
2823 maxChannel = channelId;
2824 }
2825
2826 // we need to loop backward, because removing the first channel would make the second one the new first and so on
2827 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)
2828 for(int currentChannelId = maxChannel; currentChannelId >= minChannel; --currentChannelId) {
2829 int actualChannelId = fActualChannelId[currentChannelId];
2830 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)
2831
2832 if(minChannel == maxChannel) { // we don't need to select the tab if we close all
2833 if(currentChannelId < nofTabs) {
2834 fTab->SetTab(currentChannelId + 1);
2835 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)
2836 } else {
2837 fTab->SetTab(currentChannelId - 1);
2838 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)
2839 }
2840 }
2841
2842 // remove the original active tab
2843 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)
2844 fActualChannelId.erase(fActualChannelId.begin() + currentChannelId);
2845 fActiveBins.erase(fActiveBins.begin() + currentChannelId);
2846 }
2847 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)
2848 // 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!)
2849 SelectedTab(fTab->GetCurrent());
2850 // if this was also the last source vector we initiate the last screen
2851 if(fActualChannelId.empty() || fActiveBins.empty()) {
2852 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "last channel tab done - writing files now" << std::endl; }
2854 TTimer::SingleShot(fWaitMs, "TSourceCalibration", this, "HandleTimer()");
2855 }
2856 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << RESET_COLOR << std::flush; }
2857}
2858
2860{
2861 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << DGREEN << __PRETTY_FUNCTION__ << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2862 if(fApplyToAll->IsDown()) {
2863 auto* window = CreateProgressbar("Finding current source peaks in channel #%.0f");
2864
2865 for(auto& bin : fActiveBins) {
2866 if(fChannelTab[bin - 1] != nullptr) {
2867 fChannelTab[bin - 1]->FindCalibratedPeaks();
2868 fChannelTab[bin - 1]->Calibrate();
2869 }
2870 fProgressBar->Increment(1);
2871 }
2872 window->CloseWindow();
2873 delete fProgressBar;
2874 } else if(fChannelTab[fActiveBins[fTab->GetCurrent()] - 1] != nullptr) {
2875 fChannelTab[fActiveBins[fTab->GetCurrent()] - 1]->FindCalibratedPeaks();
2876 fChannelTab[fActiveBins[fTab->GetCurrent()] - 1]->Calibrate();
2877 } else {
2878 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)
2879 }
2880 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << RESET_COLOR << std::flush; }
2881}
2882
2884{
2885 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << DGREEN << __PRETTY_FUNCTION__ << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2886 if(fChannelTab[fActiveBins[fTab->GetCurrent()] - 1] != nullptr) {
2887 fChannelTab[fActiveBins[fTab->GetCurrent()] - 1]->FindAllCalibratedPeaks();
2888 fChannelTab[fActiveBins[fTab->GetCurrent()] - 1]->Calibrate();
2889 } else {
2890 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)
2891 }
2892 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << RESET_COLOR << std::flush; }
2893}
2894
2896{
2897 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << DGREEN << __PRETTY_FUNCTION__ << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2898 auto* window = CreateProgressbar("Finding all source peaks for channel #%.0f");
2899
2900 for(auto& bin : fActiveBins) {
2901 if(fChannelTab[bin - 1] != nullptr) {
2902 fChannelTab[bin - 1]->FindAllCalibratedPeaks();
2903 fChannelTab[bin - 1]->Calibrate();
2904 } else {
2905 std::cout << __PRETTY_FUNCTION__ << ": fChannelTab[" << bin << "] is a nullptr!" << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2906 }
2907 fProgressBar->Increment(1);
2908 }
2909 window->CloseWindow();
2910 delete fProgressBar;
2911 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << RESET_COLOR << std::flush; }
2912}
2913
2915{
2916 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << DGREEN << __PRETTY_FUNCTION__ << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2917 if(fApplyToAll->IsDown()) {
2918 auto* window = CreateProgressbar("Calibrating channel #%.0f");
2919
2920 for(auto& bin : fActiveBins) {
2921 if(fChannelTab[bin - 1] != nullptr) {
2922 fChannelTab[bin - 1]->Calibrate();
2923 }
2924 fProgressBar->Increment(1);
2925 }
2926 window->CloseWindow();
2927 delete fProgressBar;
2928 } else if(fChannelTab[fActiveBins[fTab->GetCurrent()] - 1] != nullptr) {
2929 fChannelTab[fActiveBins[fTab->GetCurrent()] - 1]->Calibrate();
2930 } else {
2931 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)
2932 }
2933 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << RESET_COLOR << std::flush; }
2934}
2935
2937{
2938 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << DGREEN << __PRETTY_FUNCTION__ << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2939 if(fApplyToAll->IsDown()) {
2940 auto* window = CreateProgressbar("Recursively removing peaks for channel #%.0f");
2941
2942 for(auto& bin : fActiveBins) {
2943 if(fChannelTab[bin - 1] != nullptr) {
2944 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << "Starting recursive remove for bin " << bin << ", fChannelTab[bin - 1] " << fChannelTab[bin - 1] << std::endl; }
2945 fChannelTab[bin - 1]->RecursiveRemove(MaxResidual());
2946 }
2947 fProgressBar->Increment(1);
2948 }
2949 window->CloseWindow();
2950 delete fProgressBar;
2951 } else if(fChannelTab[fActiveBins[fTab->GetCurrent()] - 1] != nullptr) {
2952 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; }
2953 fChannelTab[fActiveBins[fTab->GetCurrent()] - 1]->RecursiveRemove(MaxResidual());
2954 } else {
2955 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)
2956 }
2957 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << RESET_COLOR << std::flush; }
2958}
2959
2961{
2962 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << DGREEN << __PRETTY_FUNCTION__ << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2963 std::ostringstream fileName;
2964 for(auto* source : fSource) {
2965 fileName << source->GetName();
2966 }
2967 fileName << ".cal";
2968 TChannel::WriteCalFile(fileName.str());
2969 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)
2971}
2972
2973TGTransientFrame* TSourceCalibration::CreateProgressbar(const char* format)
2974{
2975 auto* window = new TGTransientFrame(gClient->GetRoot(), this, fPanelWidth, fStatusbarHeight);
2976 fProgressBar = new TGHProgressBar(window, TGProgressBar::kFancy, fPanelWidth);
2977 fProgressBar->SetBarColor("lightblue");
2978 fProgressBar->SetRange(0., static_cast<Float_t>(fActiveBins.size()));
2979 fProgressBar->ShowPosition(true, false, format);
2980 window->AddFrame(fProgressBar, new TGLayoutHints(kLHintsCenterX | kLHintsCenterY | kLHintsExpandX | kLHintsExpandY, 0, 0, 0, 0));
2981 window->CenterOnParent();
2982 window->MapSubwindows();
2983 window->MapWindow();
2984 window->Resize(window->GetDefaultSize());
2985
2986 return window;
2987}
2988
2990{
2991 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << DGREEN << __PRETTY_FUNCTION__ << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
2992
2993 Print();
2994 if(fBottomFrame != nullptr) { fBottomFrame->Print(); }
2995 if(fLeftFrame != nullptr) { fLeftFrame->Print(); }
2996 if(fRightFrame != nullptr) { fRightFrame->Print(); }
2997 if(fTab != nullptr) { fTab->Print(); }
2998 std::cout << "TSourceCalibration all " << fChannelTab.size() << " channel tabs:" << std::endl;
2999 for(auto* channelTab : fChannelTab) {
3000 if(channelTab != nullptr) { channelTab->PrintLayout(); }
3001 }
3002}
3003
3005{
3006 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << DGREEN << __PRETTY_FUNCTION__ << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
3007
3008 for(auto* channelTab : fChannelTab) {
3009 if(channelTab != nullptr) { channelTab->PrintCanvases(); }
3010 }
3011}
#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
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
TGNumberEntry * fThresholdEntry
void SetSource(Int_t windowId, Int_t entryId)
std::vector< int > fActiveBins
TGTextButton * fDiscardButton
std::vector< TNucleus * > fSource
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.
static double FitRange()
TGTextButton * fStartButton
double fDefaultThreshold
The default threshold used for the peak finding algorithm, can be changed later. Co-56 source needs a...
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)