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