GamR  0.0.0
GammaROOT
Loading...
Searching...
No Matches
BackSub.cc
Go to the documentation of this file.
1#include <TBox.h>
2#include <TMarker.h>
3
4#include <utils/Utilities.hh>
5
6#include <spect/BackSub.hh>
7
8namespace GamR {
9 namespace Spect {
10
11 int BackSub2D::SetPeak(TCanvas *canvas) {
12 canvas->SetCrosshair();
13
14 canvas->Update();
15
16 GamR::Utils::Clicker click;
17 std::vector<std::string> messages = {"Click for corner 1, press any key to exit...", "Click for corner 2, press any key to exit..."};
18
19 int retval = click.GetClicks(canvas, 2, messages);
20
21 if (retval > 0) { // quit prematurely
22 std::cout << "gating cancelled" << std::endl;
23 canvas->SetCrosshair(0);
24 return retval;
25 }
26
27 double lowX = click.xs[0];
28 double lowY = click.ys[0];
29 double highX = click.xs[1];
30 double highY = click.ys[1];
31
32 if (highX < lowX) {
33 double temp = highX;
34 highX = lowX;
35 lowX = temp;
36 }
37
38 if (highY < lowY) {
39 double temp = highY;
40 highY = lowY;
41 lowY = temp;
42 }
43
44 TBox *box = new TBox(lowX, lowY, highX, highY);
45 box->SetLineColor(kRed);
46 box->SetFillStyle(0);
47 box->SetLineStyle(1);
48 box->SetLineWidth(2);
49 box->Draw("same");
50
51 //make gates
52 GamR::TK::Gate x(lowX, highX);
53 GamR::TK::Gate y(lowY, highY);
54 fPeakX = x;
55 fPeakY = y;
56
57 canvas->SetCrosshair(0);
58 return retval;
59
60 }
61
62 int BackSub2D::AddBackX(TCanvas *canvas) {
63 canvas->SetCrosshair();
64
65 canvas->Update();
66
67 GamR::Utils::Clicker click;
68 std::vector<std::string> messages = {"Click for corner 1, press any key to exit...", "Click for corner 2, press any key to exit..."};
69
70 int retval = click.GetClicks(canvas, 2, messages);
71
72 if (retval > 0) { // quit prematurely
73 std::cout << "gating cancelled" << std::endl;
74 canvas->SetCrosshair(0);
75 return retval;
76 }
77
78 double lowX = click.xs[0];
79 double lowY = click.ys[0];
80 double highX = click.xs[1];
81 double highY = click.ys[1];
82
83 if (highX < lowX) {
84 double temp = highX;
85 highX = lowX;
86 lowX = temp;
87 }
88
89 if (highY < lowY) {
90 double temp = highY;
91 highY = lowY;
92 lowY = temp;
93 }
94
95 canvas->SetCrosshair(0);
96 //std::cout << *this << std::endl;
97 //make gates
98 GamR::TK::Gate x(lowX, highX);
99 //GamR::TK::Gate y(lowY, highY);
100 fBackX.push_back(x);
101 //fPeakY = y;
102
103 TBox *box = new TBox(lowX, fPeakY.GetLow(), highX, fPeakY.GetHigh());
104 box->SetLineColor(kBlue);
105 box->SetFillStyle(0);
106 box->SetLineStyle(1);
107 box->SetLineWidth(2);
108 box->Draw("same");
109 return retval;
110
111 }
112
113 int BackSub2D::AddBackY(TCanvas *canvas) {
114 canvas->SetCrosshair();
115
116 canvas->Update();
117
118 GamR::Utils::Clicker click;
119 std::vector<std::string> messages = {"Click for corner 1, press any key to exit...", "Click for corner 2, press any key to exit..."};
120
121 int retval = click.GetClicks(canvas, 2, messages);
122 if (retval > 0) { // quit prematurely
123 std::cout << "gating cancelled" << std::endl;
124 canvas->SetCrosshair(0);
125 return retval;
126 }
127
128 double lowX = click.xs[0];
129 double lowY = click.ys[0];
130 double highX = click.xs[1];
131 double highY = click.ys[1];
132
133 if (highX < lowX) {
134 double temp = highX;
135 highX = lowX;
136 lowX = temp;
137 }
138
139 if (highY < lowY) {
140 double temp = highY;
141 highY = lowY;
142 lowY = temp;
143 }
144
145 canvas->SetCrosshair(0);
146 //std::cout << *this << std::endl;
147 //make gates
148 //GamR::TK::Gate x(lowX, highX);
149 GamR::TK::Gate y(lowY, highY);
150 fBackY.push_back(y);
151 //fBackY.push_back(y);
152
153 TBox *box = new TBox(fPeakX.GetLow(), lowY, fPeakX.GetHigh(), highY);
154 box->SetLineColor(kGreen);
155 box->SetFillStyle(0);
156 box->SetLineStyle(1);
157 box->SetLineWidth(2);
158 box->Draw("same");
159
160 return retval;
161
162 }
163
164 int BackSub2D::AddBackDiag(TCanvas *canvas) {
165 canvas->SetCrosshair();
166
167 canvas->Update();
168
169 GamR::Utils::Clicker click;
170 std::vector<std::string> messages = {"Click for corner 1, press any key to exit...", "Click for corner 2, press any key to exit..."};
171
172 int retval = click.GetClicks(canvas, 2, messages);
173 if (retval > 0) { // quit prematurely
174 std::cout << "gating cancelled" << std::endl;
175 canvas->SetCrosshair(0);
176 return retval;
177 }
178
179 double lowX = click.xs[0];
180 double lowY = click.ys[0];
181 double highX = click.xs[1];
182 double highY = click.ys[1];
183
184 if (highX < lowX) {
185 double temp = highX;
186 highX = lowX;
187 lowX = temp;
188 }
189
190 if (highY < lowY) {
191 double temp = highY;
192 highY = lowY;
193 lowY = temp;
194 }
195
196 canvas->SetCrosshair(0);
197 //std::cout << *this << std::endl;
198 //make gates
199 GamR::TK::Gate x(lowX, highX);
200 GamR::TK::Gate y(lowY, highY);
201 fBackDiagX.push_back(x);
202 fBackDiagY.push_back(y);
203
204 TBox *box = new TBox(lowX, lowY, highX, highY);
205 box->SetLineColor(kMagenta);
206 box->SetFillStyle(0);
207 box->SetLineStyle(1);
208 box->SetLineWidth(2);
209 box->Draw("same");
210
211 return retval;
212 }
213
214 int BackSub2D::AddBackBack(TCanvas *canvas) {
215 canvas->SetCrosshair();
216
217 canvas->Update();
218
219 GamR::Utils::Clicker click;
220 std::vector<std::string> messages = {"Click for corner 1, press any key to exit...", "Click for corner 2, press any key to exit..."};
221
222 int retval = click.GetClicks(canvas, 2, messages);
223 if (retval > 0) { // quit prematurely
224 std::cout << "gating cancelled" << std::endl;
225 canvas->SetCrosshair(0);
226 return retval;
227 }
228
229 double lowX = click.xs[0];
230 double lowY = click.ys[0];
231 double highX = click.xs[1];
232 double highY = click.ys[1];
233
234 if (highX < lowX) {
235 double temp = highX;
236 highX = lowX;
237 lowX = temp;
238 }
239
240 if (highY < lowY) {
241 double temp = highY;
242 highY = lowY;
243 lowY = temp;
244 }
245
246 canvas->SetCrosshair(0);
247 //std::cout << *this << std::endl;
248 //make gates
249 GamR::TK::Gate x(lowX, highX);
250 GamR::TK::Gate y(lowY, highY);
251 fBackBackX.push_back(x);
252 fBackBackY.push_back(y);
253
254 TBox *box = new TBox(lowX, lowY, highX, highY);
255 box->SetLineColor(kOrange);
256 box->SetFillStyle(0);
257 box->SetLineStyle(1);
258 box->SetLineWidth(2);
259 box->Draw("same");
260 return retval;
261 }
262
263 TH2D * BackSub2D::Subtract(TH2D *hist) {
264 int startBinX = hist->GetXaxis()->FindBin(fPeakX.GetLow());
265 int startBinY = hist->GetYaxis()->FindBin(fPeakY.GetLow());
266
267 int stopBinX = hist->GetXaxis()->FindBin(fPeakX.GetHigh());
268 int stopBinY = hist->GetYaxis()->FindBin(fPeakY.GetHigh());
269
270 subtracted = new TH2D("sub_hist", "subtracted histogram",
271 stopBinX - startBinX + 1,
272 hist->GetXaxis()->GetBinLowEdge(startBinX),
273 hist->GetXaxis()->GetBinLowEdge(stopBinX+1),
274 stopBinY - startBinY + 1,
275 hist->GetYaxis()->GetBinLowEdge(startBinY),
276 hist->GetYaxis()->GetBinLowEdge(stopBinY+1));
277
278 unsubtracted = new TH2D("unsub_hist", "unsubtracted histogram",
279 stopBinX - startBinX + 1,
280 hist->GetXaxis()->GetBinLowEdge(startBinX),
281 hist->GetXaxis()->GetBinLowEdge(stopBinX+1),
282 stopBinY - startBinY + 1,
283 hist->GetYaxis()->GetBinLowEdge(startBinY),
284 hist->GetYaxis()->GetBinLowEdge(stopBinY+1));
285
286 for (int i=startBinX; i<=stopBinX; ++i) {
287 for (int j=startBinY; j<=stopBinY; ++j) {
288 subtracted->SetBinContent(i-startBinX+1, j-startBinY+1, hist->GetBinContent(i, j));
289 subtracted->SetBinError(i-startBinX+1, j-startBinY+1, hist->GetBinError(i, j));
290
291 unsubtracted->SetBinContent(i-startBinX+1, j-startBinY+1, hist->GetBinContent(i, j));
292 unsubtracted->SetBinError(i-startBinX+1, j-startBinY+1, hist->GetBinError(i, j));
293 }
294 }
295
296 int nBinsY = 0;
297 y_sub = new TH1D("y_sub", "Y subtraction",
298 stopBinX - startBinX + 1,
299 hist->GetXaxis()->GetBinLowEdge(startBinX),
300 hist->GetXaxis()->GetBinLowEdge(stopBinX+1));
301
302 for (auto &gate : fBackY) {
303 int start = hist->GetYaxis()->FindBin(gate.GetLow());
304 int stop = hist->GetYaxis()->FindBin(gate.GetHigh());
305
306 for (int j=start; j<=stop; ++j) {
307 nBinsY += 1;
308 for (int i=startBinX; i<=stopBinX; ++i) {
309 y_sub->SetBinContent(i-startBinX+1,
310 y_sub->GetBinContent(i-startBinX+1) + hist->GetBinContent(i, j));
311 y_sub->SetBinError(i-startBinX+1,
312 std::sqrt(std::pow(y_sub->GetBinError(i-startBinX+1), 2) +
313 std::pow(hist->GetBinError(i,j), 2)));
314 }
315 }
316 }
317
318 int nBinsX = 0;
319 x_sub = new TH1D("x_sub", "X subtraction",
320 stopBinY - startBinY + 1,
321 hist->GetYaxis()->GetBinLowEdge(startBinY),
322 hist->GetYaxis()->GetBinLowEdge(stopBinY+1));
323
324 for (auto &gate : fBackX) {
325 int start = hist->GetXaxis()->FindBin(gate.GetLow());
326 int stop = hist->GetXaxis()->FindBin(gate.GetHigh());
327
328 for (int i=start; i<=stop; ++i) {
329 nBinsX += 1;
330 for (int j=startBinY; j<=stopBinY; ++j) {
331 x_sub->SetBinContent(j-startBinY+1,
332 x_sub->GetBinContent(j-startBinY+1) + hist->GetBinContent(i, j));
333 x_sub->SetBinError(j-startBinY+1,
334 std::sqrt(std::pow(x_sub->GetBinError(j-startBinY+1), 2) +
335 std::pow(hist->GetBinError(i,j), 2)));
336 }
337 }
338 }
339
340 //background-background
341 int nBinsBackground = 0;
342
343 double backback = 0;
344 double backback_err = 0;
345
346 for (int k=0; k<fBackBackX.size(); ++k) {
347 auto gateX = fBackBackX[k];
348 auto gateY = fBackBackY[k];
349
350 int startX = hist->GetXaxis()->FindBin(gateX.GetLow());
351 int stopX = hist->GetXaxis()->FindBin(gateX.GetHigh());
352 int startY = hist->GetYaxis()->FindBin(gateY.GetLow());
353 int stopY = hist->GetYaxis()->FindBin(gateY.GetHigh());
354
355 for (int i=startX; i<=stopX; ++i) {
356 for (int j=startY; j<=stopY; ++j) {
357 nBinsBackground += 1;
358 backback += hist->GetBinContent(i, j);
359 backback_err = std::sqrt(std::pow(backback_err, 2) +
360 std::pow(hist->GetBinError(i, j), 2));
361 }
362 }
363 }
364
365 if (nBinsBackground > 0) {
366 backback = backback/(double)nBinsBackground;
367 backback_err = backback_err/(double)nBinsBackground;
368 }
369
370 std::vector<int> nBinsSlice(stopBinX + stopBinY - startBinX - startBinY + 1);
371 std::vector<double> sumSlice(stopBinX + stopBinY - startBinX - startBinY + 1);
372 std::vector<double> sumErrSlice(stopBinX + stopBinY - startBinX - startBinY + 1);
373
374 int startSum = startBinX + startBinY;
375 int stopSum = stopBinX + stopBinY;
376
377 //diagonal bits
378 for (int k=0; k<fBackDiagX.size(); ++k) {
379 auto gateX = fBackDiagX[k];
380 auto gateY = fBackDiagY[k];
381
382 int startX = hist->GetXaxis()->FindBin(gateX.GetLow());
383 int stopX = hist->GetXaxis()->FindBin(gateX.GetHigh());
384 int startY = hist->GetYaxis()->FindBin(gateY.GetLow());
385 int stopY = hist->GetYaxis()->FindBin(gateY.GetHigh());
386
387 //go through each diagonal (constant sum) slice
388 for (int s=startSum; s<=stopSum; ++s) {
389 for (int i=startX; i<=stopX; ++i) {
390 int j = s - i; //so that i+j = sum
391 if (j >= startY && j <= stopY) {
392 nBinsSlice[s-startSum] += 1;
393 sumSlice[s-startSum] += hist->GetBinContent(i, j);
394 sumErrSlice[s-startSum] = std::sqrt(std::pow(sumErrSlice[s-startSum], 2) +
395 std::pow(hist->GetBinError(i, j), 2));
396 }
397 }
398 }
399 }
400
401 for (int k = 0; k<nBinsSlice.size(); ++k) {
402 if (nBinsSlice[k] > 0) {
403 sumSlice[k] = sumSlice[k]/(double)nBinsSlice[k];
404 sumErrSlice[k] = sumErrSlice[k]/(double)nBinsSlice[k];
405 }
406 }
407
408 int peakWidthY = stopBinY - startBinY + 1;
409 int peakWidthX = stopBinX - startBinX + 1;
410
411 if (nBinsX > 0) {
412 x_sub->Scale((double)1.0/(double)nBinsX);
413 }
414 if (nBinsY > 0) {
415 y_sub->Scale((double)1.0/(double)nBinsY);
416 }
417
418 for (int i=1; i<=peakWidthX; ++i) {
419 for (int j=1; j<=peakWidthY; ++j) {
420 double counts = subtracted->GetBinContent(i, j);
421 double error = subtracted->GetBinError(i, j);
422
423 if (nBinsSlice[i+j-2]>0) {
424 subtracted->SetBinContent(i, j, counts - x_sub->GetBinContent(j) - y_sub->GetBinContent(i) - sumSlice[i+j-2] + 2.0*backback);
425 subtracted->SetBinError(i, j, std::sqrt(std::pow(error, 2) +
426 std::pow(x_sub->GetBinError(j), 2) +
427 std::pow(y_sub->GetBinError(i), 2) +
428 std::pow(sumErrSlice[i+j-2], 2) +
429 std::pow(2.0*backback_err, 2)));
430 }
431 else {
432 subtracted->SetBinContent(i, j, counts - x_sub->GetBinContent(j) - y_sub->GetBinContent(i) + backback);
433 subtracted->SetBinError(i, j, std::sqrt(std::pow(error, 2) +
434 std::pow(x_sub->GetBinError(j), 2) +
435 std::pow(y_sub->GetBinError(i), 2) +
436 std::pow(backback_err, 2)));
437 }
438 }
439 }
440 return subtracted;
441 }
442
443 double BackSub2D::GetPeakCounts() {
444 double sum = 0;
445 for (int i=1; i<=subtracted->GetNbinsX(); ++i) {
446 for (int j=1; j<=subtracted->GetNbinsY(); ++j) {
447 sum += subtracted->GetBinContent(i,j);
448 }
449 }
450 return sum;
451 }
452
453 double BackSub2D::GetPeakCountsError() {
454 double err = 0;
455 for (int i=1; i<=subtracted->GetNbinsX(); ++i) {
456 for (int j=1; j<=subtracted->GetNbinsY(); ++j) {
457 err += std::pow(subtracted->GetBinError(i,j), 2);
458 }
459 }
460 return std::sqrt(err);
461 }
462
463 BackSub2D::BackSub2D(TH2D *hist) {
464 TCanvas *c1 = new TCanvas();
465
466 hist->Draw("colz2");
467
468 SetPeak(c1);
469
470 while (true) {
471 std::cout << "Select area(s) for X subtraction" << std::endl;
472 int retval = AddBackX(c1);
473 if (retval < 0) { break; }
474 }
475 while (true) {
476 std::cout << "Select area(s) for Y subtraction" << std::endl;
477 int retval = AddBackY(c1);
478 if (retval < 0) { break; }
479 }
480 while (true) {
481 std::cout << "Select area(s) for background-background region" << std::endl;
482 int retval = AddBackBack(c1);
483 if (retval < 0) { break; }
484 }
485 while (true) {
486 std::cout << "Select area(s) for diagonal subtraction" << std::endl;
487 int retval = AddBackDiag(c1);
488 if (retval < 0) { break; }
489 }
490 Subtract(hist);
491
492 std::cout << std::endl;
493 std::cout << "Peak Area: " << GetPeakCounts() << " +/- " << GetPeakCountsError() << std::endl;
494 return;
495 }
496
497
498 }
499}
GamR::Spect::BackSub2D::SetPeak
int SetPeak(GamR::TK::Gate x, GamR::TK::Gate y)
Definition BackSub.hh:35
GamR::Spect::BackSub2D::AddBackBack
int AddBackBack(GamR::TK::Gate x, GamR::TK::Gate y)
Definition BackSub.hh:79
GamR::Spect::BackSub2D::fBackDiagY
std::vector< GamR::TK::Gate > fBackDiagY
Definition BackSub.hh:23
GamR::Spect::BackSub2D::AddBackX
int AddBackX(GamR::TK::Gate x)
Definition BackSub.hh:47
GamR::Spect::BackSub2D::fPeakX
GamR::TK::Gate fPeakX
Definition BackSub.hh:15
GamR::Spect::BackSub2D::BackSub2D
BackSub2D(TH2D *hist)
Definition BackSub.cc:463
GamR::Spect::BackSub2D::fBackY
std::vector< GamR::TK::Gate > fBackY
Definition BackSub.hh:18
GamR::Spect::BackSub2D::AddBackY
int AddBackY(GamR::TK::Gate y)
Definition BackSub.hh:52
GamR::Spect::BackSub2D::fBackBackY
std::vector< GamR::TK::Gate > fBackBackY
Definition BackSub.hh:27
GamR::Spect::BackSub2D::fBackBackX
std::vector< GamR::TK::Gate > fBackBackX
Definition BackSub.hh:26
GamR::Spect::BackSub2D::GetPeakCountsError
double GetPeakCountsError()
Definition BackSub.cc:453
GamR::Spect::BackSub2D::Subtract
TH2D * Subtract(TH2D *hist)
Definition BackSub.cc:263
GamR::Spect::BackSub2D::fBackDiagX
std::vector< GamR::TK::Gate > fBackDiagX
Definition BackSub.hh:22
GamR::Spect::BackSub2D::fPeakY
GamR::TK::Gate fPeakY
Definition BackSub.hh:16
GamR::Spect::BackSub2D::fBackX
std::vector< GamR::TK::Gate > fBackX
Definition BackSub.hh:19
GamR::Spect::BackSub2D::AddBackDiag
int AddBackDiag(GamR::TK::Gate x, GamR::TK::Gate y)
Definition BackSub.hh:67
GamR::Utils::Clicker::xs
std::vector< double > xs
Definition Utilities.hh:38
GamR::Utils::Clicker::ys
std::vector< double > ys
Definition Utilities.hh:39
GamR::Utils::Clicker::GetClicks
int GetClicks(TVirtualPad *canvas, int n, std::vector< std::string > &messages, int draw=0, int print=0)
Definition Utilities.cc:127
GamR
Definition Gain.cc:19