GRSISort "v4.0.0.5"
An extension of the ROOT analysis Framework
Loading...
Searching...
No Matches
TCalibrateDescant.cxx
Go to the documentation of this file.
1#include "TCalibrateDescant.h"
2
3#include <iomanip>
4
5double MaximumRecoilEnergy(double gammaEnergy)
6{
7 return gammaEnergy - gammaEnergy / (1. + 2. * gammaEnergy / 510.998928);
8}
9
11{
12 switch(source) {
14 return MaximumRecoilEnergy(510.99895000);
16 return MaximumRecoilEnergy(1274.537);
18 return MaximumRecoilEnergy(1368.626);
20 return MaximumRecoilEnergy(2754.007);
22 return (MaximumRecoilEnergy(1173.228) + MaximumRecoilEnergy(1332.492)) / 2.;
24 return MaximumRecoilEnergy(661.657);
26 return MaximumRecoilEnergy(121.7817);
28 return MaximumRecoilEnergy(344.2785);
30 return MaximumRecoilEnergy(1408.006);
32 return MaximumRecoilEnergy(59.5409);
33 case TCalibrateDescant::ESourceType::k241Am: // photopeak, not Compton edge!
34 return 59.5409;
35 }
36 return 0.;
37}
38
40{
41 switch(source) {
43 return (MaximumRecoilEnergy(510.99895000) - MaximumRecoilEnergy(510.99895000)) / 2.;
45 return (MaximumRecoilEnergy(1274.537 + 0.007) - MaximumRecoilEnergy(1274.537 - 0.007)) / 2.;
47 return (MaximumRecoilEnergy(1368.626 + 0.005) - MaximumRecoilEnergy(1368.626 - 0.005)) / 2.;
49 return (MaximumRecoilEnergy(2754.007 + 0.011) - MaximumRecoilEnergy(2754.007 - 0.011)) / 2.;
50 case TCalibrateDescant::ESourceType::k60Co: // only taking one peaks uncertainty into account ...
51 return (MaximumRecoilEnergy(1332.492 + 0.004) + MaximumRecoilEnergy(1332.492 - 0.004)) / 2.;
53 return (MaximumRecoilEnergy(661.657 + 0.003) - MaximumRecoilEnergy(661.657 - 0.003)) / 2.;
55 return (MaximumRecoilEnergy(121.7817 + 0.0003) - MaximumRecoilEnergy(121.7817 - 0.0003)) / 2.;
57 return (MaximumRecoilEnergy(344.2785 + 0.0012) - MaximumRecoilEnergy(344.2785 - 0.0012)) / 2.;
59 return (MaximumRecoilEnergy(1408.006 + 0.003) - MaximumRecoilEnergy(1408.006 - 0.003)) / 2.;
61 return (MaximumRecoilEnergy(59.5409 + 0.0001) - MaximumRecoilEnergy(59.5409 - 0.0001)) / 2.;
62 case TCalibrateDescant::ESourceType::k241Am: // photopeak, not Compton edge!
63 return 0.0001;
64 }
65 return 0.;
66}
67
68double FullEdge(double* x, double* par) // NOLINT(performance-no-int-to-ptr, readability-non-const-parameter)
69{
70 // 0 - amplitude, 1 - position, 2 - sigma of the upper part (low x), 3 - dSigma of the lower part (high x)
71 // 4 - amplitude of gaussian peak, 5 - difference of peak position from edge position (par[2]), 6 - sigma of gaussian
72 // 7 - amplitude of noise gaussian, 8 - position of noise gaussian, 9 - sigma of noise gaussian, 10 - threshold
73 // 11 - threshold sigma, 12 - bg constant, 13 - bg amplitude, 14 - bg decay constant, 15 - cut off
74
75 // all parameters should be positive - except maybe the relative peak position???
76 for(int i = 0; i < 16; ++i) {
77 par[i] = TMath::Abs(par[i]);
78 }
79
80 if(x[0] < par[15]) { return 0.; }
81
82 double thresholdFactor = (1. + TMath::Erf((x[0] - par[10]) / par[11])) / 2.;
83
84 // create a smooth transition from sigma to sigma + dSigma
85 double sigma = par[2] + par[3] / 2. * (1. + TMath::Erf((x[0] - par[1]) / (par[2] + par[3])));
86
87 double edge = par[0] / 2. * TMath::Erfc((x[0] - par[1]) / sigma);
88 double peak = par[4] * TMath::Exp(-0.5 * TMath::Power((x[0] - (par[1] - par[5])) / par[6], 2));
89 double noise = par[7] * TMath::Exp(-0.5 * TMath::Power((x[0] - par[8]) / par[9], 2));
90
91 // Erfc goes from +2 for -inf to 0 for +inf, divide by 2 to make it go from 1 to 0
92 return thresholdFactor * (edge + peak + noise + par[12] + par[13] * TMath::Exp(-par[14] * x[0]));
93}
94
95TGHorizontalFrame* TParameterInput::Build(const std::string& name, const Int_t& baseId, const Double_t& xmin, const Double_t& xmax)
96{
97 fBaseId = baseId;
98
99 // create entries and sliders
100 // TGNumberEntry(parent window, value, digit width, id, style, attribute, limits, min, max)
101 // style - Integer, Hex, Real, Degree, etc.
102 // attribute - AnyNumber, NonNegative, or Positive
103 // limits - NoLimits, LimitMin, LimitMax, LimitMinMax
104 fLabel = new TGLabel(this, name.c_str());
105 fSlider = new TGTripleHSlider(this, 200, kDoubleScaleBoth, fBaseId, kHorizontalFrame);
106 fEntryLow = new TGNumberEntry(this, 0, 6, fBaseId, TGNumberFormat::kNESRealOne, TGNumberFormat::kNEANonNegative, TGNumberFormat::kNELLimitMinMax, xmin, xmax);
107 fEntry = new TGNumberEntry(this, 0, 6, fBaseId + 2, TGNumberFormat::kNESRealOne, TGNumberFormat::kNEANonNegative, TGNumberFormat::kNELLimitMinMax, xmin, xmax);
108 fEntryHigh = new TGNumberEntry(this, 0, 6, fBaseId + 1, TGNumberFormat::kNESRealOne, TGNumberFormat::kNEANonNegative, TGNumberFormat::kNELLimitMinMax, xmin, xmax);
109
110 // hints: hints, pad left, pad right, pad top, pad bottom
111 AddFrame(fLabel, new TGLayoutHints(kLHintsTop | kLHintsLeft, 1, 1, 1, 1));
112 AddFrame(fSlider, new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX, 1, 1, 1, 1));
113 // these go in reverse order because they are all right aligned!
114 AddFrame(fEntryHigh, new TGLayoutHints(kLHintsTop | kLHintsRight, 1, 1, 1, 1));
115 AddFrame(fEntry, new TGLayoutHints(kLHintsTop | kLHintsRight, 1, 1, 1, 1));
116 AddFrame(fEntryLow, new TGLayoutHints(kLHintsTop | kLHintsRight, 1, 1, 1, 1));
117
118 return this;
119}
120
121void TParameterInput::Set(double val)
122{
123 std::cout << __PRETTY_FUNCTION__ << ": " << val << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
124 PrintStatus("Set single initial");
125 fEntry->SetNumber(val);
126 UpdateSlider();
127 PrintStatus("Set single final");
128}
129
130void TParameterInput::Set(double val, double low, double high)
131{
132 std::cout << __PRETTY_FUNCTION__ << ": " << val << ", " << low << ", " << high << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
133 PrintStatus("Set initial");
134 fEntry->SetNumber(val);
135 fEntry->SetLimits(TGNumberFormat::kNELLimitMinMax, low, high);
136 fEntryLow->SetNumber(low);
137 fEntryLow->SetLimits(TGNumberFormat::kNELLimitMinMax, 0.1 * low, high);
138 fEntryHigh->SetNumber(high);
139 fEntryHigh->SetLimits(TGNumberFormat::kNELLimitMinMax, low, 10. * high);
140 UpdateSlider();
141 PrintStatus("Set final");
142}
143
145{
146 PrintStatus("UpdateSlider initial");
147 fSlider->SetPointerPosition(Value());
148 PrintStatus("UpdateSlider intermittent");
149 fSlider->SetPosition(LowLimit(), HighLimit());
150 fSlider->SetRange(0.1 * LowLimit(), 2. * HighLimit());
151 fSlider->PositionChanged();
152 PrintStatus("UpdateSlider final");
153}
154
156{
157 PrintStatus("UpdateEntries initial");
158 fEntry->SetNumber(fSlider->GetPointerPosition());
159 fEntryLow->SetNumber(fSlider->GetMinPosition());
160 fEntryHigh->SetNumber(fSlider->GetMaxPosition());
161 PrintStatus("UpdateEntries final");
162}
163
165{
166 fSlider->Connect("PointerPositionChanged()", "TParameterInput", this, "UpdateEntries()");
167 // fSlider->Connect("PointerPositionChanged()", "TCalibrateDescant", parent, "UpdateInitialFunction()");//creates loop on initialization? Or just set's one parameter at a time to zero and then prints all settings
168 fSlider->Connect("PositionChanged()", "TParameterInput", this, "UpdateEntries()");
169
170 fEntry->Connect("ValueSet(Long_t)", "TParameterInput", this, "UpdateSlider()");
171 fEntry->Connect("ValueChanged(Long_t)", "TParameterInput", this, "UpdateSlider()");
172 fEntry->Connect("ValueSet(Long_t)", "TCalibrateDescant", parent, "UpdateInitialFunction()");
173 fEntry->Connect("ValueChanged(Long_t)", "TCalibrateDescant", parent, "UpdateInitialFunction()");
174 fEntryLow->Connect("ValueSet(Long_t)", "TParameterInput", this, "UpdateSlider()");
175 fEntryLow->Connect("ValueChanged(Long_t)", "TParameterInput", this, "UpdateSlider()");
176 fEntryHigh->Connect("ValueSet(Long_t)", "TParameterInput", this, "UpdateSlider()");
177 fEntryHigh->Connect("ValueChanged(Long_t)", "TParameterInput", this, "UpdateSlider()");
178}
179
180Bool_t TParameterInput::ProcessMessage(Long_t msg, Long_t parameter1, Long_t parameter2)
181{
182 /// This functions deals with changes in the text fields of the TGNumberEntry as those don't seem to emit signals?
183 std::cout << __PRETTY_FUNCTION__ << ": msg " << msg << ", parameter 1 " << parameter1 << ", parameter 2 " << parameter2 << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
184 switch(GET_MSG(msg)) {
185 case kC_TEXTENTRY:
186 switch(GET_SUBMSG(msg)) {
187 case kTE_TEXTCHANGED:
188 case kTE_ENTER:
189 case kTE_TAB:
190 UpdateSlider();
191 Emit("ValueChanged()");
192 break;
193 default:
194 break;
195 }
196 break;
197 default:
198 break;
199 }
200
201 return true;
202}
203
204void TParameterInput::PrintStatus(const char* function)
205{
206 std::cout << Name() << " - " << std::setw(40) << function << ": entries - " << Value() << ", " << LowLimit() << ", " << HighLimit() << ", sliders - " << fSlider->GetPointerPosition() << ", " << fSlider->GetMinPosition() << ", " << fSlider->GetMaxPosition() << std::endl;
207}
208
210 : TGMainFrame(nullptr, 100, 100, kMainFrame | kHorizontalFrame), fMatrix(hist), fSource(source)
211{
216}
217
219{
220 // create the left frame and right frame
221 fLeftFrame = new TGVerticalFrame(this, 400, 400);
222 fRightFrame = new TGVerticalFrame(this, 400, 400);
223
224 // create canvas widgets
225 fFitCanvas = new TRootEmbeddedCanvas("FitCanvas", fLeftFrame, 200, 200);
226
227 fCalibrationCanvas = new TRootEmbeddedCanvas("CalibrationCanvas", fLeftFrame, 200, 200);
228
229 // maybe add splitter between canvases?
230
231 // create status bar
232 fStatusBar = new TGStatusBar(fLeftFrame, 400, 10, kHorizontalFrame);
233 std::array<Int_t, 3> parts = {25, 25, 50};
234 fStatusBar->SetParts(parts.data(), parts.size());
235
236 // build parameter entries
238 fAmplitude->Build("Amplitude: ", kAmplitude, 0., 1.);
240 fPosition->Build("Position: ", kPosition, 0., 1.);
242 fSigma->Build("Sigma: ", kSigma, 0., 1.);
244 fDSigma->Build("#DeltaSigma: ", kDSigma, 0., 1.);
246 fPeakAmp->Build("Peak Amplitude: ", kPeakAmp, 0., 1.);
248 fPeakPos->Build("Peak Position: ", kPeakPos, 0., 1.);
250 fPeakSigma->Build("Peak Sigma: ", kPeakSigma, 0., 1.);
252 fNoiseAmp->Build("Noise Amplitude: ", kNoiseAmp, 0., 1.);
254 fNoisePos->Build("Noise Position: ", kNoisePos, 0., 1.);
256 fNoiseSigma->Build("Noise Sigma: ", kNoiseSigma, 0., 1.);
258 fThreshold->Build("Threshold: ", kThreshold, 0., 1.);
260 fThresholdSigma->Build("Threshold Sigma: ", kThresholdSigma, 0., 1.);
262 fBgConst->Build("Background Constant: ", kBgConst, 0., 1.);
264 fBgAmp->Build("Background Amplitude: ", kBgAmp, 0., 1.);
266 fBgDecayConst->Build("Background Decay Constant: ", kBgDecayConst, 0., 1.);
268 fCutoff->Build("Cutoff: ", kCutoff, 0., 1.);
269
270 // build button frame
271 fTopButtonFrame = new TGHorizontalFrame(fRightFrame, 400, 400);
272 fBottomButtonFrame = new TGHorizontalFrame(fRightFrame, 400, 400);
273
274 // buttons
275 fPreviousButton = new TGTextButton(fTopButtonFrame, "&Previous");
276 fFitButton = new TGTextButton(fTopButtonFrame, "&Fit");
277 fNextButton = new TGTextButton(fTopButtonFrame, "&Next");
278 fUpdateInitialButton = new TGTextButton(fBottomButtonFrame, "&Update Initial Parameters");
279 fResetFitButton = new TGTextButton(fBottomButtonFrame, "&Reset Fit");
280 fSaveButton = new TGTextButton(fBottomButtonFrame, "&Save");
281
282 fTopButtonFrame->AddFrame(fPreviousButton, new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX, 1, 1, 1, 1));
283 fTopButtonFrame->AddFrame(fFitButton, new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX, 1, 1, 1, 1));
284 fTopButtonFrame->AddFrame(fNextButton, new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX, 1, 1, 1, 1));
285 fBottomButtonFrame->AddFrame(fUpdateInitialButton, new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX, 1, 1, 1, 1));
286 fBottomButtonFrame->AddFrame(fResetFitButton, new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX, 1, 1, 1, 1));
287 fBottomButtonFrame->AddFrame(fSaveButton, new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX, 1, 1, 1, 1));
288
289 // build left and right frame
290 fLeftFrame->AddFrame(fFitCanvas, new TGLayoutHints(kLHintsTop | kLHintsCenterX | kLHintsExpandX | kLHintsExpandY, 1, 1, 1, 2));
291 fLeftFrame->AddFrame(fStatusBar, new TGLayoutHints(kLHintsTop | kLHintsCenterX | kLHintsExpandX, 1, 1, 1, 2));
292 fLeftFrame->AddFrame(fCalibrationCanvas, new TGLayoutHints(kLHintsBottom | kLHintsCenterX | kLHintsExpandX | kLHintsExpandY, 1, 1, 2, 1));
293
294 fRightFrame->AddFrame(fAmplitude, new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX, 1, 1, 1, 1));
295 fRightFrame->AddFrame(fPosition, new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX, 1, 1, 1, 1));
296 fRightFrame->AddFrame(fSigma, new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX, 1, 1, 1, 1));
297 fRightFrame->AddFrame(fDSigma, new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX, 1, 1, 1, 1));
298 fRightFrame->AddFrame(fPeakAmp, new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX, 1, 1, 1, 1));
299 fRightFrame->AddFrame(fPeakPos, new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX, 1, 1, 1, 1));
300 fRightFrame->AddFrame(fPeakSigma, new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX, 1, 1, 1, 1));
301 fRightFrame->AddFrame(fNoiseAmp, new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX, 1, 1, 1, 1));
302 fRightFrame->AddFrame(fNoisePos, new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX, 1, 1, 1, 1));
303 fRightFrame->AddFrame(fNoiseSigma, new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX, 1, 1, 1, 1));
304 fRightFrame->AddFrame(fThreshold, new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX, 1, 1, 1, 1));
305 fRightFrame->AddFrame(fThresholdSigma, new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX, 1, 1, 1, 1));
306 fRightFrame->AddFrame(fBgConst, new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX, 1, 1, 1, 1));
307 fRightFrame->AddFrame(fBgAmp, new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX, 1, 1, 1, 1));
308 fRightFrame->AddFrame(fBgDecayConst, new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX, 1, 1, 1, 1));
309 fRightFrame->AddFrame(fCutoff, new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX, 1, 1, 1, 1));
310
311 // reverse order since we start these from the bottom
312 fRightFrame->AddFrame(fBottomButtonFrame, new TGLayoutHints(kLHintsBottom | kLHintsLeft | kLHintsExpandX | kLHintsExpandY, 1, 1, 1, 1));
313 fRightFrame->AddFrame(fTopButtonFrame, new TGLayoutHints(kLHintsBottom | kLHintsLeft | kLHintsExpandX | kLHintsExpandY, 1, 1, 1, 1));
314
315 AddFrame(fLeftFrame, new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX | kLHintsExpandY, 1, 1, 1, 1));
316 AddFrame(fRightFrame, new TGLayoutHints(kLHintsTop | kLHintsRight | kLHintsExpandX | kLHintsExpandY, 1, 1, 1, 1));
317
318 SetWindowName("Descant Calibration");
319
320 MapSubwindows();
321 Resize(GetDefaultSize());
322 MapWindow();
323}
324
326{
327 // connect canvases zoom function
328 fFitCanvas->GetCanvas()->Connect("RangeChanged()", "TCalibrateDescant", this, "FitCanvasZoomed()");
329 fCalibrationCanvas->GetCanvas()->Connect("RangeChanged()", "TCalibrateDescant", this, "CalibrationCanvasZoomed()");
330
331 // connect status bar info to canvases
332 fFitCanvas->GetCanvas()->Connect("ProcessedEvent(Int_t, Int_t, Int_t, TObject*)", "TCalibrateDescant", this, "Status(Int_t, Int_t, Int_t, TObject*)");
333 fCalibrationCanvas->GetCanvas()->Connect("ProcessedEvent(Int_t, Int_t, Int_t, TObject*)", "TCalibrateDescant", this, "Status(Int_t, Int_t, Int_t, TObject*)");
334
335 // connect buttons
336 fPreviousButton->Connect("Clicked()", "TCalibrateDescant", this, "Previous()");
337 fFitButton->Connect("Clicked()", "TCalibrateDescant", this, "Fit()");
338 fNextButton->Connect("Clicked()", "TCalibrateDescant", this, "Next()");
339 fUpdateInitialButton->Connect("Clicked()", "TCalibrateDescant", this, "UpdateInitialParameters()");
340 fResetFitButton->Connect("Clicked()", "TCalibrateDescant", this, "ResetFit()");
341 fSaveButton->Connect("Clicked()", "TCalibrateDescant", this, "Save()");
342
343 // connect input parameters
344 fAmplitude->Connect(this);
345 fPosition->Connect(this);
346 fSigma->Connect(this);
347 fDSigma->Connect(this);
348 fPeakAmp->Connect(this);
349 fPeakPos->Connect(this);
350 fPeakSigma->Connect(this);
351 fNoiseAmp->Connect(this);
352 fNoisePos->Connect(this);
353 fNoiseSigma->Connect(this);
354 fThreshold->Connect(this);
356 fBgConst->Connect(this);
357 fBgAmp->Connect(this);
358 fBgDecayConst->Connect(this);
359 fCutoff->Connect(this);
360}
361
362Bool_t TCalibrateDescant::ProcessMessage(Long_t msg, Long_t parameter1, Long_t parameter2)
363{
364 /// This functions deals with changes in the text fields of the TGNumberEntry as those don't seem to emit signals?
365 std::cout << __PRETTY_FUNCTION__ << ": msg " << msg << ", parameter 1 " << parameter1 << ", parameter 2 " << parameter2 << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
366 switch(GET_MSG(msg)) {
367 case kC_TEXTENTRY:
368 switch(GET_SUBMSG(msg)) {
369 case kTE_TEXTCHANGED:
370 case kTE_ENTER:
371 case kTE_TAB:
372 // UpdateSlider();
374 break;
375 default:
376 break;
377 }
378 break;
379 default:
380 break;
381 }
382
383 return true;
384}
385
387{
388 fProjections.resize(fMatrix->GetXaxis()->GetNbins());
389 fCalibrations.resize(fMatrix->GetXaxis()->GetNbins());
390 for(int i = 0; i < static_cast<int>(fProjections.size()); ++i) {
391 fProjections[i] = fMatrix->ProjectionY(Form("%s_py%d", fMatrix->GetName(), i + 1), i + 1, i + 1);
392 fProjections[i]->SetStats(false);
393 fCalibrations[i] = new TGraphErrors;
394 }
395 double xmin = fMatrix->GetYaxis()->GetBinLowEdge(1);
396 double xmax = fMatrix->GetYaxis()->GetBinLowEdge(fMatrix->GetYaxis()->GetNbins() + 1);
397 fInitial = new TF1("initial", FullEdge, xmin, xmax, 16);
398 fInitial->SetNpx(10000);
399 fInitial->SetLineColor(1);
400 fInitial->SetLineStyle(2);
401 fFit = new TF1("fit", FullEdge, xmin, xmax, 16);
402 fFit->SetNpx(10000);
403 fFit->SetLineColor(2);
404 fFit->SetLineStyle(1);
405 fEdge = new TF1("edge", FullEdge, xmin, xmax, 16);
406 fEdge->SetNpx(10000);
407 fEdge->SetLineColor(3);
408 fPeak = new TF1("peak", FullEdge, xmin, xmax, 16);
409 fPeak->SetNpx(10000);
410 fPeak->SetLineColor(4);
411 fNoise = new TF1("noise", FullEdge, xmin, xmax, 16);
412 fNoise->SetNpx(10000);
413 fNoise->SetLineColor(6);
414 fBg = new TF1("bg", FullEdge, xmin, xmax, 16);
415 fBg->SetNpx(10000);
416 fBg->SetLineColor(7);
417}
418
420{
421 std::cout << "Updating interface for current projection " << fCurrentProjection << std::endl;
424 // fProjections[fCurrentProjection]->GetListOfFunctions()->Clear();
425 // fProjections[fCurrentProjection]->GetListOfFunctions()->Add(fInitial);
426 // draw projection and calibration graph
427 fFitCanvas->GetCanvas()->cd();
429 fCalibrationCanvas->GetCanvas()->cd();
431 // force redrawing?
432 fFitCanvas->GetCanvas()->Modified();
433 fFitCanvas->GetCanvas()->Update();
434 fCalibrationCanvas->GetCanvas()->Modified();
435 fCalibrationCanvas->GetCanvas()->Update();
436 MapWindow();
437}
438
440{
441 std::cout << __PRETTY_FUNCTION__ << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
442 fInitial->FixParameter(0, fAmplitude->Value());
443 fInitial->FixParameter(1, fPosition->Value());
444 fInitial->FixParameter(2, fSigma->Value());
445 fInitial->FixParameter(3, fDSigma->Value());
446 fInitial->FixParameter(4, fPeakAmp->Value());
447 fInitial->FixParameter(5, fPeakPos->Value());
448 fInitial->FixParameter(6, fPeakSigma->Value());
449 fInitial->FixParameter(7, fNoiseAmp->Value());
450 fInitial->FixParameter(8, fNoisePos->Value());
451 fInitial->FixParameter(9, fNoiseSigma->Value());
452 fInitial->FixParameter(10, fThreshold->Value());
453 fInitial->FixParameter(11, fThresholdSigma->Value());
454 fInitial->FixParameter(12, fBgConst->Value());
455 fInitial->FixParameter(13, fBgAmp->Value());
456 fInitial->FixParameter(14, fBgDecayConst->Value());
457 fInitial->FixParameter(15, fCutoff->Value());
458 fInitial->Copy(*fFit);
459 fFit->SetLineColor(2);
460 fFit->SetLineStyle(1);
461 for(int i = 0; i < 15; ++i) {
462 std::cout << fInitial->GetParName(i) << ": " << fInitial->GetParameter(i) << std::endl;
463 }
464 if(fProjections[fCurrentProjection]->GetListOfFunctions()->GetSize() == 0) {
465 fProjections[fCurrentProjection]->GetListOfFunctions()->Add(fInitial);
466 } else {
467 std::cout << "Already " << fProjections[fCurrentProjection]->GetListOfFunctions()->GetSize() << " functions added to histogram:" << std::endl;
468 fProjections[fCurrentProjection]->GetListOfFunctions()->Print();
469 }
470 fFitCanvas->GetCanvas()->cd();
472 fFitCanvas->GetCanvas()->Modified();
473 fFitCanvas->GetCanvas()->Update();
474}
475
477{
478 /// initialize parameters base on our current histogram
479 // only use zoomed in region
480 int firstBin = fProjections[fCurrentProjection]->GetXaxis()->GetFirst();
481 int nBins = fProjections[fCurrentProjection]->GetXaxis()->GetLast();
482
483 double xmin = fProjections[fCurrentProjection]->GetXaxis()->GetBinLowEdge(firstBin);
484 double xmax = fProjections[fCurrentProjection]->GetXaxis()->GetBinLowEdge(nBins + 1);
485
486 // first find the threshold (aka first bin with more than 2 counts)
487 int threshold = 0;
488 for(threshold = firstBin; threshold < nBins; ++threshold) {
489 if(fProjections[fCurrentProjection]->GetBinContent(threshold) > 2) {
490 break;
491 }
492 }
493 // find the rough position of the edge (first bin from the top above 10 times the average)
494 // if we have a maximum bin that is not in the top or bottom 10 % of the range
495 // we use the area around it to estimate an average, otherwise we take the average of the whole range
496 double average = fProjections[fCurrentProjection]->Integral() / nBins;
497 int maxBin = fProjections[fCurrentProjection]->GetMaximumBin();
498 if(xmin + (xmax - xmin) * 0.1 < fProjections[fCurrentProjection]->GetXaxis()->GetBinCenter(maxBin) &&
499 fProjections[fCurrentProjection]->GetXaxis()->GetBinCenter(maxBin) < xmax - (xmax - xmin) * 0.1) {
500 average = fProjections[fCurrentProjection]->Integral(static_cast<Int_t>(maxBin - (xmax - xmin) * 0.1), static_cast<Int_t>(maxBin + (xmax - xmin) * 0.1)) / ((xmax - xmin) * 0.2);
501 }
502 int roughBin = 0;
503 for(roughBin = nBins; roughBin >= firstBin; --roughBin) {
504 if(fProjections[fCurrentProjection]->GetBinContent(roughBin) > 10. * average) {
505 break;
506 }
507 }
508 double roughPos = fProjections[fCurrentProjection]->GetBinCenter(roughBin);
509
510 std::cout << "Initializing parameters for current projection " << fCurrentProjection << " based on threshold bin " << threshold << ", average " << average << ", rough bin " << roughBin << ", and rough position " << roughPos << std::endl;
511
512 fInitial->FixParameter(0, 0.9 * fProjections[fCurrentProjection]->GetBinContent((roughBin - threshold) / 2)); // amplitude
513 fInitial->FixParameter(1, roughPos); // position
514 fInitial->FixParameter(2, 0.1 * roughPos); // sigma
515 fInitial->FixParameter(3, 0.2 * roughPos); // d sigma
516 fInitial->FixParameter(4, 0.3 * fProjections[fCurrentProjection]->GetBinContent((roughBin - threshold) / 2)); // peak amp
517 fInitial->FixParameter(5, 0.3 * roughPos); // peak pos
518 fInitial->FixParameter(6, 0.2 * roughPos); // peak sigma
519 fInitial->FixParameter(7, 0.1 * fProjections[fCurrentProjection]->Integral(threshold, threshold + 10) / TMath::Exp(-0.001 * threshold)); // noise amp
520 fInitial->FixParameter(8, 1.); // noise pos
521 fInitial->FixParameter(9, 1000.); // noise sigma
522 fInitial->FixParameter(10, fProjections[fCurrentProjection]->GetBinCenter(threshold)); // threshold
523 fInitial->FixParameter(11, 1.); // thres sigma
524 fInitial->FixParameter(12, 10.); // bg const
525 fInitial->FixParameter(13, 1.); // bg amp
526 fInitial->FixParameter(14, 0.0001); // bg decay const
527 fInitial->FixParameter(15, fProjections[fCurrentProjection]->GetBinLowEdge(threshold)); // cutoff
528
529 fInitial->SetParName(0, fAmplitude->Name());
530 fInitial->SetParName(1, fPosition->Name());
531 fInitial->SetParName(2, fSigma->Name());
532 fInitial->SetParName(3, fDSigma->Name());
533 fInitial->SetParName(4, fPeakAmp->Name());
534 fInitial->SetParName(5, fPeakPos->Name());
535 fInitial->SetParName(6, fPeakSigma->Name());
536 fInitial->SetParName(7, fNoiseAmp->Name());
537 fInitial->SetParName(8, fNoisePos->Name());
538 fInitial->SetParName(9, fNoiseSigma->Name());
539 fInitial->SetParName(10, fThreshold->Name());
540 fInitial->SetParName(11, fThresholdSigma->Name());
541 fInitial->SetParName(12, fBgConst->Name());
542 fInitial->SetParName(13, fBgAmp->Name());
543 fInitial->SetParName(14, fBgDecayConst->Name());
544 fInitial->SetParName(15, fCutoff->Name());
545
546 fAmplitude->Set(fInitial->GetParameter(0), xmin, xmax);
547 fPosition->Set(fInitial->GetParameter(1), xmin, xmax);
548 fSigma->Set(fInitial->GetParameter(2), xmin, xmax);
549 fDSigma->Set(fInitial->GetParameter(3), xmin, xmax);
550 fPeakAmp->Set(fInitial->GetParameter(4), xmin, xmax);
551 fPeakPos->Set(fInitial->GetParameter(5), xmin, xmax);
552 fPeakSigma->Set(fInitial->GetParameter(6), xmin, xmax);
553 fNoiseAmp->Set(fInitial->GetParameter(7), xmin, xmax);
554 fNoisePos->Set(fInitial->GetParameter(8), xmin, xmax);
555 fNoiseSigma->Set(fInitial->GetParameter(9), xmin, xmax);
556 fThreshold->Set(fInitial->GetParameter(10), xmin, xmax);
557 fThresholdSigma->Set(fInitial->GetParameter(11), xmin, xmax);
558 fBgConst->Set(fInitial->GetParameter(12), xmin, xmax);
559 fBgAmp->Set(fInitial->GetParameter(13), xmin, xmax);
560 fBgDecayConst->Set(fInitial->GetParameter(14), xmin, xmax);
561 fCutoff->Set(fInitial->GetParameter(15), xmin, xmax);
562
563 fInitial->Copy(*fFit);
564 fFit->SetLineColor(2);
565 fFit->SetLineStyle(1);
566}
567
569{
570 std::cout << __PRETTY_FUNCTION__ << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
573}
574
576{
577 std::cout << __PRETTY_FUNCTION__ << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
578 if(fCurrentProjection + 1 < static_cast<int>(fProjections.size())) { ++fCurrentProjection; }
580}
581
583{
584 fFitCanvas->GetCanvas()->cd();
585
586 // set parameters from TParameterInput
587 fFit->SetParameter(0, fAmplitude->Value());
588 fFit->SetParLimits(0, fAmplitude->LowLimit(), fAmplitude->HighLimit());
589 fFit->SetParameter(1, fPosition->Value());
590 fFit->SetParLimits(1, fPosition->LowLimit(), fPosition->HighLimit());
591 fFit->SetParameter(2, fSigma->Value());
592 fFit->SetParLimits(2, fSigma->LowLimit(), fSigma->HighLimit());
593 fFit->SetParameter(3, fDSigma->Value());
594 fFit->SetParLimits(3, fDSigma->LowLimit(), fDSigma->HighLimit());
595 fFit->SetParameter(4, fPeakAmp->Value());
596 fFit->SetParLimits(4, fPeakAmp->LowLimit(), fPeakAmp->HighLimit());
597 fFit->SetParameter(5, fPeakPos->Value());
598 fFit->SetParLimits(5, fPeakPos->LowLimit(), fPeakPos->HighLimit());
599 fFit->SetParameter(6, fPeakSigma->Value());
600 fFit->SetParLimits(6, fPeakSigma->LowLimit(), fPeakSigma->HighLimit());
601 fFit->SetParameter(7, fNoiseAmp->Value());
602 fFit->SetParLimits(7, fNoiseAmp->LowLimit(), fNoiseAmp->HighLimit());
603 fFit->SetParameter(8, fNoisePos->Value());
604 fFit->SetParLimits(8, fNoisePos->LowLimit(), fNoisePos->HighLimit());
605 fFit->SetParameter(9, fNoiseSigma->Value());
606 fFit->SetParLimits(9, fNoiseSigma->LowLimit(), fNoiseSigma->HighLimit());
607 fFit->SetParameter(10, fThreshold->Value());
608 fFit->SetParLimits(10, fThreshold->LowLimit(), fThreshold->HighLimit());
609 fFit->SetParameter(11, fThresholdSigma->Value());
610 fFit->SetParLimits(11, fThresholdSigma->LowLimit(), fThresholdSigma->HighLimit());
611 fFit->SetParameter(12, fBgConst->Value());
612 fFit->SetParLimits(12, fBgConst->LowLimit(), fBgConst->HighLimit());
613 fFit->SetParameter(13, fBgAmp->Value());
614 fFit->SetParLimits(13, fBgAmp->LowLimit(), fBgAmp->HighLimit());
615 fFit->SetParameter(14, fBgDecayConst->Value());
616 fFit->SetParLimits(14, fBgDecayConst->LowLimit(), fBgDecayConst->HighLimit());
617 fFit->SetParameter(15, fCutoff->Value());
618 fFit->SetParLimits(15, fCutoff->LowLimit(), fCutoff->HighLimit());
619
620 std::cout << "done setting parameters" << std::endl;
621
622 // perform the fit
623 fProjections[fCurrentProjection]->Fit(fFit, "RQN"); // R - use function range, Q - quiet, N - do not store/draw
624
625 std::cout << "done fitting" << std::endl;
626
627 // get parameters and parameter errors
628
629 // update the component functions
630 for(int i = 0; i < 15; ++i) {
631 std::cout << "parameter " << i << " - " << fFit->GetParName(i) << ": " << fFit->GetParameter(i) << " +- " << fFit->GetParError(i) << std::endl;
632 double par = fFit->GetParameter(i);
633 std::cout << fEdge << ", " << fPeak << ", " << fNoise << ", " << fBg << std::endl;
634 fEdge->Print();
635 fPeak->Print();
636 fNoise->Print();
637 fBg->Print();
638 switch(i) {
639 case 0:
640 fEdge->SetParameter(i, par);
641 fPeak->SetParameter(i, 0.);
642 fNoise->SetParameter(i, 0.);
643 fBg->SetParameter(i, 0.);
644 break;
645 case 4:
646 fEdge->SetParameter(i, 0.);
647 fPeak->SetParameter(i, par);
648 fNoise->SetParameter(i, 0.);
649 fBg->SetParameter(i, 0.);
650 break;
651 case 7:
652 fEdge->SetParameter(i, 0.);
653 fPeak->SetParameter(i, 0.);
654 fNoise->SetParameter(i, par);
655 fBg->SetParameter(i, 0.);
656 break;
657 case 12:
658 fEdge->SetParameter(i, 0.);
659 fPeak->SetParameter(i, 0.);
660 fNoise->SetParameter(i, 0.);
661 fBg->SetParameter(i, par);
662 break;
663 case 13:
664 fEdge->SetParameter(i, 0.);
665 fPeak->SetParameter(i, 0.);
666 fNoise->SetParameter(i, 0.);
667 fBg->SetParameter(i, par);
668 break;
669 default:
670 fEdge->SetParameter(i, par);
671 fPeak->SetParameter(i, par);
672 fNoise->SetParameter(i, par);
673 fBg->SetParameter(i, par);
674 break;
675 };
676 }
677
678 std::cout << "done updating components" << std::endl;
679
680 // add all functions to histogram (after clearing all functions from it)
681 if(fProjections[fCurrentProjection]->GetListOfFunctions()->GetSize() < 5) {
682 fProjections[fCurrentProjection]->GetListOfFunctions()->Add(fBg);
683 fProjections[fCurrentProjection]->GetListOfFunctions()->Add(fNoise);
684 fProjections[fCurrentProjection]->GetListOfFunctions()->Add(fPeak);
685 fProjections[fCurrentProjection]->GetListOfFunctions()->Add(fEdge);
686 fProjections[fCurrentProjection]->GetListOfFunctions()->Add(fFit);
687 } else {
688 std::cout << "Already " << fProjections[fCurrentProjection]->GetListOfFunctions()->GetSize() << " functions added to histogram:" << std::endl;
689 fProjections[fCurrentProjection]->GetListOfFunctions()->Print();
690 }
691
692 std::cout << "done adding functions" << std::endl;
693
694 // add new point to calibration graph
695 AddCalibrationPoint(fFit->GetParameter(1), fFit->GetParError(1));
696
697 fFitCanvas->GetCanvas()->Modified();
698 fFitCanvas->GetCanvas()->Update();
699 fCalibrationCanvas->GetCanvas()->Modified();
700 fCalibrationCanvas->GetCanvas()->Update();
701
703
704 std::cout << "done" << std::endl;
705}
706
708{
709 std::cout << __PRETTY_FUNCTION__ << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
710 fAmplitude->Set(fFit->GetParameter(0));
711 fPosition->Set(fFit->GetParameter(1));
712 fSigma->Set(fFit->GetParameter(2));
713 fDSigma->Set(fFit->GetParameter(3));
714 fPeakAmp->Set(fFit->GetParameter(4));
715 fPeakPos->Set(fFit->GetParameter(5));
716 fPeakSigma->Set(fFit->GetParameter(6));
717 fNoiseAmp->Set(fFit->GetParameter(7));
718 fNoisePos->Set(fFit->GetParameter(8));
719 fNoiseSigma->Set(fFit->GetParameter(9));
720 fThreshold->Set(fFit->GetParameter(10));
721 fThresholdSigma->Set(fFit->GetParameter(11));
722 fBgConst->Set(fFit->GetParameter(12));
723 fBgAmp->Set(fFit->GetParameter(13));
724 fBgDecayConst->Set(fFit->GetParameter(14));
725 fCutoff->Set(fFit->GetParameter(15));
726}
727
729{
730 std::cout << __PRETTY_FUNCTION__ << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
731 fInitial->Copy(*fFit);
732 fFit->SetLineColor(2);
733 fFit->SetLineStyle(1);
734 fFitCanvas->GetCanvas()->Modified();
735 fFitCanvas->GetCanvas()->Update();
736}
737
739{
740 std::cout << __PRETTY_FUNCTION__ << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
741}
742
743void TCalibrateDescant::AddCalibrationPoint(double value, double uncertainty)
744{
745 // check if point already exists for this energy and update it if so
746 int i = 0;
747 Double_t* x = fCalibrations[fCurrentProjection]->GetX();
748 for(i = 0; i < fCalibrations[fCurrentProjection]->GetN(); ++i) {
749 if(x[i] == SourceEnergy(fSource)) {
750 fCalibrations[fCurrentProjection]->SetPoint(i, x[i], value);
751 fCalibrations[fCurrentProjection]->SetPointError(i, fCalibrations[fCurrentProjection]->GetErrorX(i), uncertainty);
752 return;
753 }
754 }
755 // create new point
757 fCalibrations[fCurrentProjection]->SetPointError(i, SourceEnergyUncertainty(fSource), uncertainty);
758}
759
761{
762 std::cout << __PRETTY_FUNCTION__ << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
763 // update range of fit functions
764 Double_t xmin = 0.;
765 Double_t ymin = 0.;
766 Double_t xmax = 0.;
767 Double_t ymax = 0.;
768 fFitCanvas->GetCanvas()->GetRange(xmin, ymin, xmax, ymax);
769
770 std::cout << "updating ranges to " << xmin << " - " << xmax << std::endl;
771 fInitial->SetRange(xmin, xmax);
772 fFit->SetRange(xmin, xmax);
773 fEdge->SetRange(xmin, xmax);
774 fPeak->SetRange(xmin, xmax);
775 fNoise->SetRange(xmin, xmax);
776 fBg->SetRange(xmin, xmax);
777
778 // InitializeParameters();// or UpdateInterface()?
779 // fProjections[fCurrentProjection]->GetListOfFunctions()->Clear();
780 // fProjections[fCurrentProjection]->GetListOfFunctions()->Add(fInitial);
781}
782
784{
785 std::cout << __PRETTY_FUNCTION__ << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
786 // nothing to do for this one?
787}
788
789void TCalibrateDescant::Status(Int_t px, Int_t py, Int_t, TObject* selected)
790{
791 // std::cout<<__PRETTY_FUNCTION__<<std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
792 fStatusBar->SetText(selected->GetName(), 0);
793 fStatusBar->SetText(selected->GetTitle(), 1);
794 fStatusBar->SetText(selected->GetObjectInfo(px, py), 2);
795}
TH1D * hist
Definition UserFillObj.h:3
TParameterInput * fPeakPos
TParameterInput * fNoiseAmp
TParameterInput * fPosition
TCalibrateDescant(TH2 *hist, const ESourceType &source=ESourceType::k137Cs)
TParameterInput * fPeakSigma
TGStatusBar * fStatusBar
TGTextButton * fFitButton
TGVerticalFrame * fLeftFrame
void AddCalibrationPoint(double value, double uncertainty)
TParameterInput * fSigma
TGTextButton * fPreviousButton
TGHorizontalFrame * fTopButtonFrame
TParameterInput * fBgDecayConst
TGHorizontalFrame * fBottomButtonFrame
TParameterInput * fBgAmp
TParameterInput * fThresholdSigma
std::vector< TH1D * > fProjections
TGTextButton * fSaveButton
TGTextButton * fUpdateInitialButton
TParameterInput * fBgConst
TParameterInput * fAmplitude
TRootEmbeddedCanvas * fFitCanvas
TGTextButton * fResetFitButton
TParameterInput * fDSigma
std::vector< TGraphErrors * > fCalibrations
TParameterInput * fPeakAmp
Bool_t ProcessMessage(Long_t msg, Long_t parameter1, Long_t parameter2) override
TGVerticalFrame * fRightFrame
TRootEmbeddedCanvas * fCalibrationCanvas
TGTextButton * fNextButton
TParameterInput * fNoisePos
TParameterInput * fCutoff
void Status(Int_t px, Int_t py, Int_t pz, TObject *selected)
TParameterInput * fThreshold
TParameterInput * fNoiseSigma
TGNumberEntry * fEntry
TGNumberEntry * fEntryLow
TGNumberEntry * fEntryHigh
TGTripleHSlider * fSlider
double HighLimit() const
const char * Name() const
void Connect(TCalibrateDescant *parent)
TGHorizontalFrame * Build(const std::string &name, const Int_t &baseId, const Double_t &xmin, const Double_t &xmax)
Bool_t ProcessMessage(Long_t msg, Long_t parameter1, Long_t parameter2) override
double Value() const
double LowLimit() const
void PrintStatus(const char *)
void Set(double val)
double SourceEnergyUncertainty(const TCalibrateDescant::ESourceType &source)
double FullEdge(double *x, double *par)
double SourceEnergy(const TCalibrateDescant::ESourceType &source)
double MaximumRecoilEnergy(double gammaEnergy)