GamR  0.0.0
GammaROOT
Loading...
Searching...
No Matches
Peak.cc
Go to the documentation of this file.
1#include <TMatrixD.h>
2#include <TGraphErrors.h>
3#include <TMarker.h>
4#include <HFitInterface.h>
5#include <Math/WrappedMultiTF1.h>
6#include <Fit/BinData.h>
7#include <Fit/Fitter.h>
8
9#include <spect/Cut.hh>
10#include "Peak.hh"
11
12namespace GamR {
13 namespace TK {
14 BPeak::BPeak() { linear = new TF1(); backlow = new TF1(); backhigh = new TF1(); };
15 BPeak::~BPeak() { delete linear; delete backlow; delete backhigh; };
16
17 void BPeak::Set(TH1 *hist, GamR::TK::Gate peak, TF1 *background) {
18 int iStart = peak.GetBinLow(hist);
19 int iStop = peak.GetBinHigh(hist);
20
21 double sumX = 0;
22 double sumXY = 0;
23 double sumY2E2 = 0;
24 double sumY = 0;
25 double sumY2 = 0;
26 double sumYerr = 0;
27 int nData = 0;
28 double max = 0;
29 double sumBinErr = 0;
30 double areaErr = 0;
31 TH1D *h = new TH1D("h", "h", iStop-iStart+1, hist->GetBinLowEdge(iStart), hist->GetBinLowEdge(iStop)+hist->GetBinWidth(iStop));
32 for (int i=iStart; i<=iStop; ++i) {
33 double binWidth = hist->GetBinWidth(i);
34 double binCent = hist->GetBinCenter(i);
35 double binErr = hist->GetBinError(i);
36 double backErr = 0;
37 if (background->GetNpar() < 4) { //crude to suppress error message for now
38 backErr = background->IntegralError(hist->GetBinLowEdge(i), hist->GetBinLowEdge(i)+binWidth)/binWidth;
39 }
40 binErr = std::sqrt(binErr*binErr + backErr*backErr);
41 sumBinErr = std::sqrt(sumBinErr*sumBinErr + binErr*binErr);
42 areaErr = std::sqrt(areaErr*areaErr + binErr*binErr*binWidth*binWidth);
43 /*
44 ++nData;
45 double binCont = (hist->GetBinContent(i) - background->Eval(binCent))*binWidth;
46*binWidth;
47 double sigc = background->GetParError(0);
48 double sigm = background->GetParError(1);
49
50 std::cout << binErr << " " << backErr << std::endl;
51 binErr = std::sqrt(binErr*binErr + backErr*backErr);
52 sumYerr = std::sqrt(sumYerr*sumYerr + binErr*binErr);
53 sumX += binCent;
54 sumY += binCont;
55 sumY2 += binCont*binCont;
56 sumXY += binCont*binCent;
57 sumY2E2 += binCont*binCont*binErr*binErr;
58 if (binCont/binWidth > max) { max = binCont/binWidth; }
59 */
60 h->SetBinContent(i-iStart+1, hist->GetBinContent(i)-background->Eval(binCent));
61 h->SetBinError(i-iStart+1, binErr);
62 }
63
64 //double mean = sumXY/sumY;
65 Cent = h->GetMean();
66 //stddev/neff
67 //double neff = sumY*sumY/sumY2;
68
69 //CentError = std::sqrt(sumY2E2/(sumY*sumY));
70 /*
71 double sumYdiff2 = 0;
72 double sumdiff2 = 0;
73 for (int i=iStart; i<=iStop; ++i) {
74 double binWidth = hist->GetBinWidth(i);
75 double binCent = hist->GetBinCenter(i);
76 double binCont = (hist->GetBinContent(i) - background->Eval(binCent))*binWidth;
77 sumYdiff2 += binCont*std::pow(binCent - mean, 2);
78 sumdiff2 += std::pow(binCent - mean, 2);
79 }
80 */
81
82 //double stddev = std::sqrt(1.0/(iStop-iStart+1)*sumdiff2);
83 CentError = h->GetMeanError();//stddev/neff;
84 Width = 2.3548*h->GetStdDev();//2.3548*std::sqrt(sumY/(sumY*sumY - sumY2)*sumYdiff2); //assuming something vaguely gaussian to get a FWHM here
85 WidthError = 2.3548*h->GetStdDevError();
86 Amp = h->GetMaximum();
87 AmpError = h->GetBinError(h->GetMaximumBin());
88 Area = h->Integral("width");
89 AreaError = areaErr;
90 Counts = h->Integral();
91 CountsError = sumBinErr;
92
93 h->Delete();
94 }
95
96 void BPeak::Set(TH1 *hist, double x1, double y1, double x2, double y2) {
97 TGraph *graph = new TGraph();
98 graph->SetPoint(0, x1, y1);
99 graph->SetPoint(1, x2, y2);
100
101 delete linear;
102 linear = new TF1("linear", "[0]+[1]*x", x1, x2);
103 double m = (y2-y1)/(x2-x1);
104 double c = y2 - m*x2;
105
106 linear->SetParameter(0, c);
107 linear->SetParameter(1, m);
108 graph->Fit(linear, "RQN");
109
110 GamR::TK::Gate peak(x1, x2);
111 mPeak = peak;
112
113 Set(hist, peak, linear);
114
115 linear->Draw("same");
116
117 this->Draw("same");
118
119 graph->Delete();
120 }
121
122 void BPeak::Set(TH1 *hist, GamR::TK::Gate peak, std::vector<GamR::TK::Gate > background, Option_t *foption, Option_t *option) {
123 //make TGraphErrors;
124 if (background.size()==0) { std::cout << "Must have at least one background region" << std::endl; return; }
125
126 TString opts(option);
127 opts.ToLower();
128
129 if (opts.Contains("s")) {
130 SetCont(hist, peak, background, foption, option);
131 return;
132 }
133 if (opts.Contains("p")) {
134 SetPoint(hist, peak, background, foption, option);
135 return;
136 }
137 if (opts.Contains("x")) {
138 SetStep(hist, peak, background, foption, option);
139 return;
140 }
141
142 SetCont(hist, peak, background, foption, option);
143 return;
144 }
145
146 void BPeak::SetCont(TH1 *hist, GamR::TK::Gate peak, std::vector<GamR::TK::Gate > background, Option_t *foption, Option_t *option) { //where a single continuous background should be fitted
147
148 TString opts(option);
149 opts.ToLower();
150 if (opts.Contains("v")) {
151 std::cout << "Setting peak with regions:" << std::endl;
152 std::cout << peak.GetLow() << " " << peak.GetHigh() << " ";
153 for (auto &b : background) {
154 std::cout << b.GetLow() << " " << b.GetHigh() << " ";
155 }
156 std::cout << std::endl;
157 }
158
159 mPeak = peak;
160 mBackgrounds = background;
161
162 TGraphErrors *back = new TGraphErrors();
163 double lowest=1e9;
164 double highest=-1e9;
165 for (int b=0; b<background.size(); ++b) {
166 int iStart = background[b].GetBinLow(hist);
167 int iStop = background[b].GetBinHigh(hist);
168 double last_err = -1.0;
169 for (int i=iStart; i<=iStop; ++i) {
170 int nPoint = back->GetN();
171 double x = hist->GetBinCenter(i);
172 double y = hist->GetBinContent(i);
173 if (hist->GetBinError(i) > 0.0) {
174 back->SetPoint(nPoint, x, y);
175 back->SetPointError(nPoint, 0, hist->GetBinError(i));
176 last_err = hist->GetBinError(i);
177 }
178 else {
179 //hacky and possibly statistically dubious?
180 //necessary for cases with low counts where bins with no counts should be
181 //taken into account
182 if (last_err > 0.0) {
183 back->SetPoint(nPoint, x, y);
184 back->SetPointError(nPoint, 0, last_err);
185 }
186 }
187 if (x < lowest) { lowest = x; }
188 if (x > highest) {highest = x; }
189 }
190 }
191
192 delete linear;
193
194 if (opts.Contains("c")) {
195 linear = new TF1("linear", "[0]", lowest, highest);
196 double y1 = hist->GetBinContent(hist->FindBin(lowest));
197 double y2 = hist->GetBinContent(hist->FindBin(highest));
198 linear->SetParameter(0, (y1+y2)/2.0);
199 }
200 else { //default behaviour
201 linear = new TF1("linear", "[0]+[1]*x", lowest, highest);
202 double x1 = hist->GetBinCenter(hist->FindBin(lowest));
203 double x2 = hist->GetBinCenter(hist->FindBin(highest));
204 double y1 = hist->GetBinContent(hist->FindBin(lowest));
205 double y2 = hist->GetBinContent(hist->FindBin(highest));
206
207 double w = x2 - x1;
208
209 double m = (y2-y1)/(x2-x1);
210 double c = y2 - m*x2;
211
212 linear->SetParameter(0, c);
213 linear->SetParameter(1, m);
214 }
215
216 back->Fit(linear, ("RQN"+(std::string)foption).c_str());
217
218 Set(hist, peak, linear);
219
220 linear->Draw("same");
221
222 this->Draw("same");
223
224 back->Delete();
225 }
226
227 void BPeak::SetPoint(TH1 *hist, GamR::TK::Gate peak, std::vector<GamR::TK::Gate > background, Option_t *foption, Option_t *option) { //where upper and lower background regions are fitted and extrpolated to the boundaries of the peak area
228
229 TString opts(option);
230 opts.ToLower();
231 if (opts.Contains("v")) {
232 std::cout << "Setting peak with regions:" << std::endl;
233 std::cout << peak.GetLow() << " " << peak.GetHigh() << " ";
234 for (auto &b : background) {
235 std::cout << b.GetLow() << " " << b.GetHigh() << " ";
236 }
237 std::cout << std::endl;
238 }
239
240 mPeak = peak;
241 mBackgrounds = background;
242
243 TGraphErrors *back_low = new TGraphErrors();
244 TGraphErrors *back_high = new TGraphErrors();
245 TGraphErrors *backs[2] = {back_low, back_high};
246 int back_ind = -1;
247 double lowest[2]={1e9, 1e9};
248 double highest[2]={-1e9, -1e9};
249 for (int b=0; b<background.size(); ++b) {
250 float av_pos = (background[b].GetLow() + background[b].GetHigh())/2.0;
251 if (av_pos < peak.GetLow()) { back_ind = 0;}
252 else if (av_pos > peak.GetHigh()) { back_ind = 1; }
253 else { std::cout << "Invalid background region, ignoring" << std::endl; continue; }
254
255 int iStart = background[b].GetBinLow(hist);
256 int iStop = background[b].GetBinHigh(hist);
257 double last_err = -1.0;
258 for (int i=iStart; i<=iStop; ++i) {
259 int nPoint = backs[back_ind]->GetN();
260 double x = hist->GetBinCenter(i);
261 double y = hist->GetBinContent(i);
262 if (hist->GetBinError(i) > 0.0) {
263 backs[back_ind]->SetPoint(nPoint, x, y);
264 backs[back_ind]->SetPointError(nPoint, 0, hist->GetBinError(i));
265 last_err = hist->GetBinError(i);
266 }
267 else {
268 //hacky and possibly statistically dubious?
269 //necessary for cases with low counts where bins with no counts should be
270 //taken into account
271 if (last_err > 0.0) {
272 backs[back_ind]->SetPoint(nPoint, x, y);
273 backs[back_ind]->SetPointError(nPoint, 0, last_err);
274 }
275 }
276 if (x < lowest[back_ind]) { lowest[back_ind] = x; }
277 if (x > highest[back_ind]) {highest[back_ind] = x; }
278 }
279 }
280
281 Lowest = lowest[0];
282 Highest = highest[1];
283
284 delete linear;
285 delete backlow;
286 delete backhigh;
287
288 if (opts.Contains("c")) {
289 backlow = new TF1("backlow", "[0]", lowest[0], highest[0]);
290 double y1 = hist->GetBinContent(hist->FindBin(lowest[0]));
291 double y2 = hist->GetBinContent(hist->FindBin(highest[0]));
292 backlow->SetParameter(0, (y1+y2)/2.0);
293 backs[0]->Fit(backlow, ("RQN"+(std::string)foption).c_str());
294
295 backhigh = new TF1("backlow", "[0]", lowest[1], highest[1]);
296 y1 = hist->GetBinContent(hist->FindBin(lowest[1]));
297 y2 = hist->GetBinContent(hist->FindBin(highest[1]));
298 backhigh->SetParameter(0, (y1+y2)/2.0);
299 backs[1]->Fit(backhigh, ("RQN"+(std::string)foption).c_str());
300
301 }
302 else { //default behaviour
303 backlow = new TF1("backlow", "[0]+[1]*x", lowest[0], highest[0]);
304 double x1 = hist->GetBinCenter(hist->FindBin(lowest[0]));
305 double x2 = hist->GetBinCenter(hist->FindBin(highest[0]));
306 double y1 = hist->GetBinContent(hist->FindBin(lowest[0]));
307 double y2 = hist->GetBinContent(hist->FindBin(highest[0]));
308
309 double w = x2 - x1;
310
311 double m = (y2-y1)/(x2-x1);
312 double c = y2 - m*x2;
313
314 backlow->SetParameter(0, c);
315 backlow->SetParameter(1, m);
316 backs[0]->Fit(backlow, ("RQN"+(std::string)foption).c_str());
317
318 backhigh = new TF1("backhigh", "[0]+[1]*x", lowest[1], highest[1]);
319 x1 = hist->GetBinCenter(hist->FindBin(lowest[1]));
320 x2 = hist->GetBinCenter(hist->FindBin(highest[1]));
321 y1 = hist->GetBinContent(hist->FindBin(lowest[1]));
322 y2 = hist->GetBinContent(hist->FindBin(highest[1]));
323
324 w = x2 - x1;
325
326 m = (y2-y1)/(x2-x1);
327 c = y2 - m*x2;
328
329 backhigh->SetParameter(0, c);
330 backhigh->SetParameter(1, m);
331 backs[1]->Fit(backhigh, ("RQN"+(std::string)foption).c_str());
332
333 }
334
335 linear = new TF1("linear", "[0] + [1]*x", peak.GetLow(), peak.GetHigh());
336 double y1 = backlow->Eval(hist->GetBinCenter(hist->FindBin(peak.GetLow())));
337 double y2 = backhigh->Eval(hist->GetBinCenter(hist->FindBin(peak.GetHigh())));
338 double x1 = peak.GetLow();
339 double x2 = peak.GetHigh();
340
341 double m = (y2-y1)/(x2-x1);
342 double c = y2 - m*x2;
343
344 linear->SetParameter(0, c);
345 linear->SetParameter(1, m);
346
347 TGraphErrors *back_graph = new TGraphErrors();
348 int i = hist->FindBin(x1);
349 double binWidth = hist->GetBinWidth(i);
350 double backErr = backlow->IntegralError(hist->GetBinLowEdge(i), hist->GetBinLowEdge(i)+binWidth)/binWidth;
351 back_graph->SetPoint(0, x1, y1);
352 back_graph->SetPointError(0, 0, backErr);
353
354 i = hist->FindBin(x2);
355 binWidth = hist->GetBinWidth(i);
356 backErr = backhigh->IntegralError(hist->GetBinLowEdge(i), hist->GetBinLowEdge(i)+binWidth)/binWidth;
357 back_graph->SetPoint(0, x2, y2);
358 back_graph->SetPointError(0, 0, backErr);
359
360 back_graph->Fit(linear, ("RQN"+(std::string)foption).c_str());
361
362 Set(hist, peak, linear);
363
364 linear->Draw("same");
365
366 this->Draw("same");
367
368 backs[0]->Delete();
369 backs[1]->Delete();
370 }
371
372 void BPeak::SetStep(TH1 *hist, GamR::TK::Gate peak, std::vector<GamR::TK::Gate > background, Option_t *foption, Option_t *option) { //where upper and lower background regions are fitted and extrpolated to the boundaries of the peak area, a step is created between them with width and centroid of the peak
373 SetPoint(hist, peak, background, foption, option);
374
375 TString opts(option);
376 opts.ToLower();
377
378 //now get centroid and width
379 double cent = GetCent();
380 double cent_err = GetCentError();
381 double sigma = GetWidth()/2.3548;
382 double sigma_err = GetWidthError()/2.3548;
383
384 delete linear;
385
386 if (opts.Contains("c")) {
387 linear = new TF1("linear", "0.5 * ROOT::Math::erfc((x-[0])/(sqrt(2)*[1])) * ([2]) + (1.0 - 0.5 * ROOT::Math::erfc((x-[0])/(sqrt(2)*[1]))) * ([3])", Lowest, Highest);
388
389 linear->SetParameter(0, cent);
390 linear->SetParameter(1, sigma);
391 linear->SetParameter(2, backlow->GetParameter(0));
392 linear->SetParameter(3, backhigh->GetParameter(0));
393
394 linear->SetParError(0, cent_err);
395 linear->SetParError(1, sigma_err);
396 linear->SetParError(2, backlow->GetParError(0));
397 linear->SetParError(3, backhigh->GetParError(0));
398 }
399 else {
400 linear = new TF1("linear", "0.5 * ROOT::Math::erfc((x-[0])/(sqrt(2)*[1])) * ([2] + [3]*x) + (1.0 - 0.5 * ROOT::Math::erfc((x-[0])/(sqrt(2)*[1]))) * ([4] + [5]*x)", Lowest, Highest);
401
402 linear->SetParameter(0, cent);
403 linear->SetParameter(1, sigma);
404 linear->SetParameter(2, backlow->GetParameter(0));
405 linear->SetParameter(3, backlow->GetParameter(1));
406 linear->SetParameter(4, backhigh->GetParameter(0));
407 linear->SetParameter(5, backhigh->GetParameter(1));
408
409 linear->SetParError(0, cent_err);
410 linear->SetParError(1, sigma_err);
411 linear->SetParError(2, backlow->GetParError(0));
412 linear->SetParError(3, backlow->GetParError(1));
413 linear->SetParError(4, backhigh->GetParError(0));
414 linear->SetParError(5, backhigh->GetParError(1));
415 }
416 Set(hist, peak, linear);
417
418 //do this again
419 for (int i=0; i<3; ++i) {
420 cent = GetCent();
421 cent_err = GetCentError();
422 sigma = GetWidth()/2.3548;
423 sigma_err = GetWidthError()/2.3548;
424
425 if (opts.Contains("c")) {
426 linear->SetParameter(0, cent);
427 linear->SetParameter(1, sigma);
428 linear->SetParameter(2, backlow->GetParameter(0));
429 linear->SetParameter(3, backhigh->GetParameter(0));
430
431 linear->SetParError(0, cent_err);
432 linear->SetParError(1, sigma_err);
433 linear->SetParError(2, backlow->GetParError(0));
434 linear->SetParError(3, backhigh->GetParError(0));
435 }
436 else {
437 linear->SetParameter(0, cent);
438 linear->SetParameter(1, sigma);
439 linear->SetParameter(2, backlow->GetParameter(0));
440 linear->SetParameter(3, backlow->GetParameter(1));
441 linear->SetParameter(4, backhigh->GetParameter(0));
442 linear->SetParameter(5, backhigh->GetParameter(1));
443
444 linear->SetParError(0, cent_err);
445 linear->SetParError(1, sigma_err);
446 linear->SetParError(2, backlow->GetParError(0));
447 linear->SetParError(3, backlow->GetParError(1));
448 linear->SetParError(4, backhigh->GetParError(0));
449 linear->SetParError(5, backhigh->GetParError(1));
450 }
451
452 Set(hist, peak, linear);
453 }
454
455 linear->Draw("same");
456
457 this->Draw("same");
458 }
459
460 void BPeak::Paint(Option_t *option /*="same"*/) {
461 mPeak.SetLineColor(kBlue);
462 mPeak.SetFillStyle(3003);
463 mPeak.SetFillColor(kBlue);
464 mPeak.Paint("");
465
466 for (int i=0; i<mBackgrounds.size(); ++i) {
467 mBackgrounds[i].SetLineColor(kRed);
468 mBackgrounds[i].SetFillStyle(3003);
469 mBackgrounds[i].SetFillColor(kRed);
470 mBackgrounds[i].Paint("");
471 }
472 // TAttLine::SetLineColor(kRed);
473 // TAttLine::Modify();
474
475 // double x1[2];
476 // double y1[2];
477 // double x2[2];
478 // double y2[2];
479
480 // // std::vector<double[2]> bx1;
481 // // std::vector<double[2]> bx2;
482
483 // x1[0] = mPeak.GetLow();
484 // x1[1] = mPeak.GetLow();
485
486 // y1[0] = gPad->GetUymin();
487 // y1[1] = gPad->GetUymax();
488
489 // x2[0] = mPeak.GetHigh();
490 // x2[1] = mPeak.GetHigh();
491
492 // y2[0] = gPad->GetUymin();
493 // y2[1] = gPad->GetUymax();
494
495 // gPad->PaintPolyLine(2, &x1[0], &y1[0], option);
496 // gPad->PaintPolyLine(2, &x2[0], &y2[0], option);
497
498 // for (int i=0; i<mBackgrounds.size(); ++i) {
499 // TAttLine::SetLineColor(kBlue);
500 // TAttLine::Modify();
501 // x1[0] = mBackgrounds[i].GetLow();
502 // x1[1] = mBackgrounds[i].GetLow();
503
504 // x2[0] = mBackgrounds[i].GetHigh();
505 // x2[1] = mBackgrounds[i].GetHigh();
506
507 // gPad->PaintPolyLine(2, &x1[0], &y1[0], option);
508 // gPad->PaintPolyLine(2, &x2[0], &y2[0], option);
509 // }
510 }
511
512 Int_t BPeak::DistancetoPrimitive(Int_t px, Int_t py) {
513 (void)py;
514 double lowDist = std::abs(gPad->XtoPixel(mPeak.GetLow()) - px);
515 double highDist = std::abs(gPad->XtoPixel(mPeak.GetHigh()) - px);
516
517 for (int i=0; i<mBackgrounds.size(); ++i) {
518 double low = mBackgrounds[i].GetLow();
519 double high = mBackgrounds[i].GetHigh();
520
521 low = std::abs(gPad->XtoPixel(low) - px);
522 high = std::abs(gPad->XtoPixel(high) - px);
523
524 if (low < lowDist) { lowDist = low; }
525 if (high < highDist) { highDist = high; }
526 }
527
528 return fmin(lowDist, highDist);
529 }
530
531 void BPeak::ExecuteEvent(Int_t event, Int_t px, Int_t py) {
532 (void)py;
533 if (!gPad) { return; }
534 int kMaxDiff = 20;
535
536 if (!gPad->IsEditable()) { return; }
537
538 switch(event) {
539 case kMouseMotion:
540 if (DistancetoPrimitive(px,py) < kMaxDiff) {
541 gPad->SetCursor(kPointer);
542 return;
543 }
544 break;
545 }
546 }
547
548 void BPeak::Set(TH1 *hist, GamR::TK::Gate peak, GamR::TK::Gate background, Option_t *foption, Option_t *option) {
549 std::vector<GamR::TK::Gate > bg;
550 bg.push_back(background);
551 Set(hist, peak, bg, foption, option);
552 }
553
554 void BPeak::Set(TH1 *hist, Option_t *foption, Option_t *option) {
555 GamR::TK::Gate peak;
556 std::cout << "Select peak region: " << std::endl;
557 peak.SetGate(gPad->GetCanvas(), "x");
558 std::cout << "Select background regions:" << std::endl;
559 std::vector<GamR::TK::Gate > background;
560 while(true) {
561 GamR::TK::Gate bg;
562 int retval = bg.SetGate(gPad->GetCanvas(), "x");
563 if (retval>0) { break; }
564 background.push_back(bg);
565 }
566 Set(hist, peak, background, foption, option);
567 }
568
569 void BPeak::Set(TH2 *hist, Option_t *foption, Option_t *option) {
570 GamR::TK::Gate prompt;
571 std::cout << "Select prompt region: " << std::endl;
572 TH1D *projx = hist->ProjectionX();
573 projx->Draw();
574 prompt.SetGate(gPad->GetCanvas(), "x");
575 std::cout << "Select non-prompt regions: " << std::endl;
576 std::vector<GamR::TK::Gate > nonprompt;
577 while(true) {
578 GamR::TK::Gate np;
579 int retval = np.SetGate(gPad->GetCanvas(), "x");
580 if (retval>0) { break; }
581 nonprompt.push_back(np);
582 }
583 if (nonprompt.size() == 0) { std::cout << "Must have at least one non-prompt region!" << std::endl; return; }
584 //make prompt and non-prompt histograms
585 TH1D *hPrompt = prompt.ApplyX(hist, "hPrompt");
586 double p_width = prompt.GetBinWidth(projx);
587 TH1D *hNonPrompt = nonprompt[0].ApplyX(hist, "hNonPrompt0");
588 double np_width = nonprompt[0].GetBinWidth(projx);
589 for (int i=1; i<nonprompt.size(); ++i) {
590 hNonPrompt->Add(nonprompt[i].ApplyX(hist, ("hNonPrompt"+std::to_string(i)).c_str()));
591 np_width += nonprompt[i].GetBinWidth(projx);
592 }
593 hNonPrompt->Scale(p_width/np_width);
594
595 //now plot them both (?)
596 hPrompt->Draw("hist");
597 hNonPrompt->SetLineColor(kRed);
598 hNonPrompt->Draw("hist same");
599 //make LFit
600 BPeak::LFit lfit(hNonPrompt);
601
602 //select peak and background regions
603 GamR::TK::Gate peak;
604 std::cout << "Select peak region: " << std::endl;
605 peak.SetGate(gPad->GetCanvas(), "x");
606 std::cout << "Select background regions:" << std::endl;
607 std::vector<GamR::TK::Gate > background;
608 double lowest = 1e9;
609 double highest = -1e9;
610 while(true) {
611 GamR::TK::Gate bg;
612 int retval = bg.SetGate(gPad->GetCanvas(), "x");
613 if (retval>0) { break; }
614 background.push_back(bg);
615 if (bg.GetLow() < lowest) {
616 lowest = bg.GetLow();
617 }
618 if (bg.GetHigh() > highest) {
619 highest = bg.GetHigh();
620 }
621 }
622 if (background.size() == 0) { std::cout << "Must have at least one background region!" << std::endl; return; }
623
624 ROOT::Fit::DataRange range(background[0].GetLow(), background[0].GetHigh());
625 for (int i=1; i<background.size(); ++i) {
626 range.AddRange(0, background[i].GetLow(), background[i].GetHigh());
627 }
628
629 delete linear;
630 linear = new TF1("linear", lfit, lowest, highest, 2);
631 //set initial guesses in linear here
632
633 double x1 = hPrompt->GetBinCenter(hPrompt->FindBin(lowest));
634 double x2 = hPrompt->GetBinCenter(hPrompt->FindBin(highest));
635 double y1 = hPrompt->GetBinContent(hPrompt->FindBin(lowest)) - hNonPrompt->GetBinContent(hNonPrompt->FindBin(lowest));
636 double y2 = hPrompt->GetBinContent(hPrompt->FindBin(highest)) - hNonPrompt->GetBinContent(hNonPrompt->FindBin(highest));
637
638 double w = x2 - x1;
639
640 double m = (y2-y1)/(x2-x1);
641 double c = y2 - m*x2;
642
643 linear->SetParameter(0, c);
644 linear->SetParameter(1, m);
645
646 ROOT::Fit::DataOptions opt;
647 opt.fUseEmpty = true;
648 ROOT::Fit::BinData data(opt, range);
649 ROOT::Fit::FillData(data, hPrompt);
650 ROOT::Math::WrappedMultiTF1 fitFunc(*linear, linear->GetNdim());
651 ROOT::Fit::Fitter fitter;
652 fitter.SetFunction(fitFunc, false);
653 fitter.LikelihoodFit(data);
654
655 fitter.Result().Print(std::cout);
656 //linear->Draw("same");
657
658 Set(hPrompt, peak, linear);
659
660 TH2D* hBackSub = (TH2D*)hPrompt->Clone("hBackSub");
661 hBackSub->Add(hNonPrompt, -1.0);
662 hBackSub->Draw("hist");
663
664 TF1 *lindraw = new TF1("lindraw", "[0] + [1]*x", lowest, highest);
665 lindraw->SetParameter(0, linear->GetParameter(0));
666 lindraw->SetParameter(1, linear->GetParameter(1));
667 lindraw->Draw("same");
668
669
670 }
671
672 void BPeak::Print() {
673 TString sPrintString;
674 sPrintString.Form("Centroid Height FWHM Area Counts ");
675 std::printf("%s\n", sPrintString.Data());
676
677 char line[120];
678 *line = '\0';
679 GamR::Utils::wrresult(line, GetCent(),
680 GetCentError(), 15);
681 GamR::Utils::wrresult(line, GetAmp(),
682 GetAmpError(), 30);
683 GamR::Utils::wrresult(line, GetWidth(),
684 GetWidthError(), 45);
685 GamR::Utils::wrresult(line, GetArea(),
686 GetAreaError(), 60);
687 GamR::Utils::wrresult(line, GetCounts(),
688 GetCountsError(), 75);
689 std::printf("%s\n", line);
690 }
691
692 GaussianPeak::GaussianPeak(double Low, double High) : FPeak(Low, High, Gaussian) {
693 fPeakFunc = new TF1("GaussianPeak", "[0]*exp(-pow(x-[1], 2)/(2*pow([2],2)))",
694 Low, High);
695 fPeakFunc->SetParName(0, "Amp");
696 fPeakFunc->SetParName(1, "Cent");
697 fPeakFunc->SetParName(2, "Sigma");
698 }
699
700 StepGaussianPeak::StepGaussianPeak(double Low, double High) : FPeak(Low, High, StepGaussian) {
701 fPeakFunc = new TF1("StepGaussianPeak", "[0]*exp(-pow(x-[1], 2)/(2*pow([2],2))) + [0]*([3]/100)*std::erfc((x-[1])/(sqrt(2)*[2]))",
702 Low, High);
703 fPeakFunc->SetParName(0, "Amp");
704 fPeakFunc->SetParName(1, "Cent");
705 fPeakFunc->SetParName(2, "Sigma");
706 fPeakFunc->SetParName(3, "H");
707
708 }
709
710 OneTailGaussianPeak::OneTailGaussianPeak(double Low, double High) : FPeak(Low, High, OneTailGaussian) {
711 fPeakFunc = new TF1("OneTailGaussianPeak", "[0]*(1 - [3]/100)*exp(-pow(x-[1],2)/(2*pow([2],2))) + "
712 "[0]*([3]/100)*exp((x-[1])/[4])*std::erfc((x-[1])/"
713 "(sqrt(2)*[2]) + [2]/(sqrt(2)*[4]))",
714 Low, High);
715 fPeakFunc->SetParName(0, "Amp");
716 fPeakFunc->SetParName(1, "Cent");
717 fPeakFunc->SetParName(2, "Sigma");
718 fPeakFunc->SetParName(3, "R");
719 fPeakFunc->SetParName(4, "Beta");
720
721 }
722
723 double OneTailGaussianPeak::GetAreaError() {
724 TMatrixD covMatrix(5,5);
725 for (int i=0; i<5; i++) {
726 for (int j=0; j<5; j++) {
727 covMatrix[i][j] = 0;
728 }
729 }
730
731 for(int i=0; i<5; i++) {
732 covMatrix[i][i] = fPeakFunc->GetParError(i) * fPeakFunc->GetParError(i);
733 }
734
735 double areaError = fPeakFunc->IntegralError(fLow, fHigh, fPeakFunc->GetParameters(), covMatrix.GetMatrixArray());
736 return areaError;
737 }
738
739 OneTailStepGaussianPeak::OneTailStepGaussianPeak(double Low, double High) : FPeak(Low, High, OneTailStepGaussian) {
740 fPeakFunc = new TF1("OneTailStepGaussianPeak", "[0]*(1 - [3]/100)*exp(-pow(x-[1],2)/(2*pow([2],2))) + "
741 "[0]*([3]/100)*exp((x-[1])/[4])*std::erfc((x-[1])/"
742 "(sqrt(2)*[2]) + [2]/(sqrt(2)*[4])) + "
743 "[0]*([5]/100)*std::erfc((x-[1])/(sqrt(2)*[2]))",
744 Low, High);
745 fPeakFunc->SetParName(0, "Amp");
746 fPeakFunc->SetParName(1, "Cent");
747 fPeakFunc->SetParName(2, "Sigma");
748 fPeakFunc->SetParName(3, "R");
749 fPeakFunc->SetParName(4, "Beta");
750 fPeakFunc->SetParName(5, "H");
751 }
752
753 double OneTailStepGaussianPeak::GetArea() {
754 //construct OneTailGaussianPeak and use that to get area
755 GamR::TK::OneTailGaussianPeak purePeak(fLow, fHigh);
756
757 for (int i=0; i<5; ++i) {
758 purePeak.GetFunc()->SetParameter(i, GetFunc()->GetParameter(i));
759 purePeak.GetFunc()->SetParError(i, GetFunc()->GetParError(i));
760 }
761
762 double area = purePeak.GetArea();
763 return area;
764 }
765
766 double OneTailStepGaussianPeak::GetAreaError() {
767 //construct OneTailGaussianPeak and use that to get area
768 GamR::TK::OneTailGaussianPeak purePeak(fLow, fHigh);
769
770 for (int i=0; i<5; ++i) {
771 purePeak.GetFunc()->SetParameter(i, GetFunc()->GetParameter(i));
772 purePeak.GetFunc()->SetParError(i, GetFunc()->GetParError(i));
773 }
774
775 double area = purePeak.GetAreaError();
776 return area;
777 }
778
779 TwoTailGaussianPeak::TwoTailGaussianPeak(double Low, double High) : FPeak(Low, High, TwoTailGaussian) {
780 fPeakFunc = new TF1("TwoTailGaussianPeak",
781 "[0]*(1 -([3]/100 + "
782 "[5]/100))*exp(-pow(x-[1],2)/(2*pow([2],2))) + "
783 "[0]*([3]/100)*exp((x-[1])/[4])*std::erfc((x-[1])/"
784 "(sqrt(2)*[2]) + [2]/(sqrt(2)*[4])) + "
785 "[0]*([5]/100)*exp((x-[1])/[6])*std::erfc((x-[1])/"
786 "(sqrt(2)*[2]) + [2]/(sqrt(2)*[6]))",
787 Low, High);
788 fPeakFunc->SetParName(0, "Amp");
789 fPeakFunc->SetParName(1, "Cent");
790 fPeakFunc->SetParName(2, "Sigma");
791 fPeakFunc->SetParName(3, "R1");
792 fPeakFunc->SetParName(4, "Beta1");
793 fPeakFunc->SetParName(5, "R2");
794 fPeakFunc->SetParName(6, "Beta2");
795 }
796
797 double TwoTailGaussianPeak::GetAreaError() {
798 TMatrixD covMatrix(7,7);
799 for (int i=0; i<7; i++) {
800 for (int j=0; j<7; j++) {
801 covMatrix[i][j] = 0;
802 }
803 }
804
805 for(int i=0; i<7; i++) {
806 covMatrix[i][i] = fPeakFunc->GetParError(i) * fPeakFunc->GetParError(i);
807 }
808
809 double areaError = fPeakFunc->IntegralError(fLow, fHigh, fPeakFunc->GetParameters(), covMatrix.GetMatrixArray());
810 return areaError;
811 }
812
813 TwoTailStepGaussianPeak::TwoTailStepGaussianPeak(double Low, double High) : FPeak(Low, High, TwoTailStepGaussian) {
814 fPeakFunc = new TF1("TwoTailStepGaussianPeak",
815 "[0]*(1 -([3]/100 + "
816 "[5]/100))*exp(-pow(x-[1],2)/(2*pow([2],2))) + "
817 "[0]*([3]/100)*exp((x-[1])/[4])*std::erfc((x-[1])/"
818 "(sqrt(2)*[2]) + [2]/(sqrt(2)*[4])) + "
819 "[0]*([5]/100)*exp((x-[1])/[6])*std::erfc((x-[1])/"
820 "(sqrt(2)*[2]) + [2]/(sqrt(2)*[6])) + "
821 "[0]*([7]/100)*std::erfc((x-[1])/(sqrt(2)*[2]))",
822 Low, High);
823 fPeakFunc->SetParName(0, "Amp");
824 fPeakFunc->SetParName(1, "Cent");
825 fPeakFunc->SetParName(2, "Sigma");
826 fPeakFunc->SetParName(3, "R1");
827 fPeakFunc->SetParName(4, "Beta1");
828 fPeakFunc->SetParName(5, "R2");
829 fPeakFunc->SetParName(6, "Beta2");
830 fPeakFunc->SetParName(7, "H");
831 }
832
833 double TwoTailStepGaussianPeak::GetArea() {
834 //construct TwoTailGaussianPeak and use that to get area
835 GamR::TK::TwoTailGaussianPeak purePeak(fLow, fHigh);
836
837 for (int i=0; i<7; ++i) {
838 purePeak.GetFunc()->SetParameter(i, GetFunc()->GetParameter(i));
839 purePeak.GetFunc()->SetParError(i, GetFunc()->GetParError(i));
840 }
841
842 double area = purePeak.GetArea();
843 return area;
844 }
845
846 double TwoTailStepGaussianPeak::GetAreaError() {
847 //construct TwoTailGaussianPeak and use that to get area
848 GamR::TK::TwoTailGaussianPeak purePeak(fLow, fHigh);
849
850 for (int i=0; i<7; ++i) {
851 purePeak.GetFunc()->SetParameter(i, GetFunc()->GetParameter(i));
852 purePeak.GetFunc()->SetParError(i, GetFunc()->GetParError(i));
853 }
854
855 double area = purePeak.GetAreaError();
856 return area;
857 }
858
859 BPeak *FitBPeak(TCanvas *canvas, Option_t *foption, Option_t *option) {
860 TH1D *hist = GamR::Utils::GetHist1D(canvas);
861 if (!hist) { return NULL; }
862
863 BPeak *peak = new BPeak();
864 peak->Set(hist, foption, option);
865
866 peak->Print();
867 return peak;
868 }
869
870 BPeak *ClickBPeak(TCanvas *canvas) {
871 TH1D *hist = GamR::Utils::GetHist1D(canvas);
872 if (!hist) { return NULL; }
873
874 BPeak *peak = new BPeak();
875 //get two (x,y) coordinates
876
877 canvas->SetCrosshair();
878
879 canvas->Update();
880
881 GamR::Utils::Clicker click;
882 std::vector<std::string> messages = {"Click for lower bound of region, press any key to exit...", "Click for upper bound of region, press any key to exit..."};
883
884 int retval = click.GetClicks(canvas, 2, messages);
885
886 if (retval > 0) { // quit prematurely
887 std::cout << "gating cancelled" << std::endl;
888 canvas->SetCrosshair(0);
889 return peak;
890 }
891
892 double lowX = click.xs[0];
893 double lowY = click.ys[0];
894 double highX = click.xs[1];
895 double highY = click.ys[1];
896
897 canvas->SetCrosshair(0);
898
899 peak->Set(hist, lowX, lowY, highX, highY);
900
901 peak->Print();
902 return peak;
903 }
904
905
906 OneTailGaussianPeak::OneTailGaussianPeak() {}
907 OneTailGaussianPeak::~OneTailGaussianPeak() {}
908
909 }
910}
py
void py(TVirtualPad *canvas)
Definition Abbreviations.cc:195
px
void px(TVirtualPad *canvas)
Definition Abbreviations.cc:194
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::TK::OneTailStepGaussianPeak::GetAreaError
double GetAreaError()
Definition Peak.cc:766
GamR::TK::BPeak::GetAmp
double GetAmp()
Definition Peak.hh:100
GamR::TK::Gate::Paint
void Paint(Option_t *option="same")
Definition Gate.cc:721
GamR::TK::Gate::GetLow
double GetLow() const
Definition Gate.hh:70
GamR::TK::Gate::GetBinWidth
int GetBinWidth(TH1 *hist) const
Definition Gate.hh:75
GamR::TK::FPeak::fPeakFunc
TF1 * fPeakFunc
Definition Peak.hh:131
GamR::TK::BPeak::GetCent
double GetCent()
Definition Peak.hh:102
GamR::TK::BPeak::GetCentError
double GetCentError()
Definition Peak.hh:103
GamR::TK::OneTailStepGaussianPeak::OneTailStepGaussianPeak
OneTailStepGaussianPeak()
Definition Peak.hh:269
GamR::TK::FPeak::fHigh
double fHigh
Definition Peak.hh:133
GamR::TK::BPeak::GetAreaError
double GetAreaError()
Definition Peak.hh:97
GamR::TK::Gate::GetBinLow
int GetBinLow(TH1 *hist) const
Definition Gate.hh:73
GamR::TK::OneTailGaussianPeak::GetArea
double GetArea()
Definition Peak.hh:239
GamR::TK::BPeak::SetCont
void SetCont(TH1 *hist, GamR::TK::Gate peak, std::vector< GamR::TK::Gate > background, Option_t *foption="", Option_t *option="")
Definition Peak.cc:146
GamR::TK::OneTailGaussianPeak::OneTailGaussianPeak
OneTailGaussianPeak()
Definition Peak.cc:906
GamR::TK::TwoTailGaussianPeak::GetArea
double GetArea()
Definition Peak.hh:311
GamR::TK::TwoTailStepGaussianPeak::GetAreaError
double GetAreaError()
Definition Peak.cc:846
GamR::TK::Gate::ApplyX
TH1D * ApplyX(TH2 *hist) const
Definition Gate.cc:218
GamR::TK::BPeak::SetStep
void SetStep(TH1 *hist, GamR::TK::Gate peak, std::vector< GamR::TK::Gate > background, Option_t *foption="", Option_t *option="")
Definition Peak.cc:372
GamR::TK::OneTailStepGaussianPeak::GetArea
double GetArea()
Definition Peak.cc:753
GamR::TK::StepGaussianPeak::StepGaussianPeak
StepGaussianPeak()
Definition Peak.hh:193
GamR::TK::BPeak::GetWidthError
double GetWidthError()
Definition Peak.hh:105
GamR::TK::TwoTailStepGaussianPeak::TwoTailStepGaussianPeak
TwoTailStepGaussianPeak()
Definition Peak.hh:342
GamR::TK::Gate::SetGate
int SetGate()
Definition Gate.cc:76
GamR::TK::BPeak::GetCounts
double GetCounts()
Definition Peak.hh:98
GamR::TK::FPeak::FPeak
FPeak()
Definition Peak.hh:139
GamR::TK::TwoTailStepGaussianPeak::GetArea
double GetArea()
Definition Peak.cc:833
GamR::TK::OneTailGaussianPeak::GetAreaError
double GetAreaError()
Definition Peak.cc:723
GamR::TK::BPeak::GetAmpError
double GetAmpError()
Definition Peak.hh:101
GamR::TK::BPeak::GetWidth
double GetWidth()
Definition Peak.hh:104
GamR::TK::BPeak::Print
void Print()
Definition Peak.cc:672
GamR::TK::BPeak::BPeak
BPeak()
Definition Peak.cc:14
GamR::TK::OneTailGaussianPeak::~OneTailGaussianPeak
~OneTailGaussianPeak()
Definition Peak.cc:907
GamR::TK::Gate::GetBinHigh
int GetBinHigh(TH1 *hist) const
Definition Gate.hh:74
GamR::TK::ClickBPeak
BPeak * ClickBPeak(TCanvas *canvas)
Definition Peak.cc:870
GamR::TK::BPeak::GetCountsError
double GetCountsError()
Definition Peak.hh:99
GamR::TK::GaussianPeak::GaussianPeak
GaussianPeak()
Definition Peak.hh:158
GamR::TK::Gate::GetHigh
double GetHigh() const
Definition Gate.hh:71
GamR::TK::BPeak::SetPoint
void SetPoint(TH1 *hist, GamR::TK::Gate peak, std::vector< GamR::TK::Gate > background, Option_t *foption="", Option_t *option="")
Definition Peak.cc:227
GamR::TK::BPeak::GetArea
double GetArea()
Definition Peak.hh:96
GamR::TK::TwoTailGaussianPeak::TwoTailGaussianPeak
TwoTailGaussianPeak()
Definition Peak.hh:308
GamR::TK::FPeak::GetFunc
TF1 * GetFunc()
Definition Peak.hh:143
GamR::TK::BPeak::~BPeak
~BPeak()
Definition Peak.cc:15
GamR::TK::TwoTailGaussianPeak::GetAreaError
double GetAreaError()
Definition Peak.cc:797
GamR::TK::FPeak::fLow
double fLow
Definition Peak.hh:132
GamR::TK::FitBPeak
BPeak * FitBPeak(TCanvas *canvas, Option_t *foption, Option_t *option)
Definition Peak.cc:859
GamR::TK::OneTailStepGaussian
@ OneTailStepGaussian
Definition Peak.hh:33
GamR::TK::StepGaussian
@ StepGaussian
Definition Peak.hh:33
GamR::TK::Gaussian
@ Gaussian
Definition Peak.hh:33
GamR::TK::TwoTailGaussian
@ TwoTailGaussian
Definition Peak.hh:33
GamR::TK::OneTailGaussian
@ OneTailGaussian
Definition Peak.hh:33
GamR::TK::TwoTailStepGaussian
@ TwoTailStepGaussian
Definition Peak.hh:33
GamR::Utils::GetHist1D
TH1D * GetHist1D(TVirtualPad *canvas)
Definition Utilities.cc:182
GamR::Utils::wrresult
int wrresult(char *out, float value, float err, int minlen)
Definition Utilities.cc:241
GamR
Definition Gain.cc:19