GamR  0.0.0
GammaROOT
Loading...
Searching...
No Matches
TimeWalk.cc
Go to the documentation of this file.
1#include <iomanip>
2#include <iostream>
3#include <sstream>
4
5#include <TMinuit.h>
6#include <TPaveText.h>
7#include <TRandom3.h>
8
9#include <toolkit/MonteCarlo.hh>
10
11#include "TimeWalk.hh"
12
13// ignoring npar, gin, iflag in GlobalChiSquare::doit
14// TMinuit expects these extra variables that we don't use
15#pragma GCC diagnostic ignored "-Wunused-parameter"
16
17namespace GamR {
18namespace Processor {
19
26TimeWalk::TimeWalk(std::vector<Int_t> ids) : fDetIds(ids)
27{
28 // Fit Function String
29 // x = E_j
30 // [0] = a_j, [1] = b_j, [2] = f_j
31 // [3] = a_i, [4] = b_i, [5] = f_i
32 // [6] = E_i,ref
33
34 std::string funcstr("([0]/sqrt(x+[1])+[2]) - ([3]/sqrt([6]+[4])+[5])");
35 for (const auto &i : fDetIds) {
36 for (const auto &j : fDetIds) {
37 if (i == j) {
38 continue;
39 }
40 std::string pairing = std::to_string(i) + "_" + std::to_string(j);
41 fFits[{i, j}] = std::make_unique<TF1>(("TCF_" + pairing).c_str(), funcstr.c_str(), 0, 4096);
42 }
43 }
44
45} // TimeWalk::TimeWalk
46
47TimeWalk::~TimeWalk() {}
48
57Int_t TimeWalk::AddDataSet(std::string fn)
58{
59 TFile f(fn.c_str());
60 if (!f.IsOpen() || f.IsZombie()) {
61 return 1;
62 }
63
64 GamR::Nucleus::Photon *ref;
65 f.GetObject("TC_REF", ref);
66 f.Close();
67
68 if (ref) {
69 fReference[fn] = std::unique_ptr<GamR::Nucleus::Photon>(ref);
70 return 0;
71 } else {
72 return 2;
73 }
74}
75
83void TimeWalk::Draw(TCanvas *c, Bool_t fits, std::string saveas)
84{
85 c->Clear();
86 c->Divide(fDetIds.size(), fDetIds.size());
87 c->SetTitle("Time Walk Calibration Graphs");
88 int i = 0;
89 int j = 0;
90 for (auto &idi : fDetIds) {
91 for (auto &idj : fDetIds) {
92 if (idi == idj) {
93 ++j;
94 continue;
95 }
96 c->cd(i * fDetIds.size() + j + 1);
97 // auto mg = std::make_shared<TMultiGraph>("mg", "");
98 Bool_t first = kTRUE;
99 for (auto &dataset : fReference) {
100 auto datafn = dataset.first;
101 if (first) {
102 fGraphs[datafn][{idi, idj}]->Draw("AP");
103 first = kFALSE;
104 } else {
105 fGraphs[datafn][{idi, idj}]->Draw("P SAME");
106 }
107 }
108 // c -> BuildLegend();
109 c->Update();
110 ++j;
111 }
112 ++i;
113 j = 0;
114 }
115 std::vector<Bool_t> equation(fDetIds.size());
116 if (fits) {
117 i = 0;
118 j = 0;
119 for (auto &idi : fDetIds) {
120 for (auto &idj : fDetIds) {
121 if (idi == idj) {
122 ++j;
123 continue;
124 }
125 if (!equation[i]) {
126 c->cd(i * fDetIds.size() + i + 1);
127 TPaveText *pavetext = new TPaveText(0.1, 0.1, 0.9, 0.9);
128 pavetext->SetTextFont(53);
129 pavetext->SetTextSize(12);
130 pavetext->SetTextAlign(11);
131 std::stringstream ss;
132 ss.precision(10);
133 ss << "Detector " << idi;
134 pavetext->AddText(ss.str().c_str());
135 ss.str("");
136 ss << "a_{" << idi << "} = " << fFitResults[{idi, 0}];
137 pavetext->AddText(ss.str().c_str());
138 ss.str("");
139 ss << "b_{" << idi << "} = " << fFitResults[{idi, 1}];
140 pavetext->AddText(ss.str().c_str());
141 ss.str("");
142 ss << "f_{" << idi << "} = " << fFitResults[{idi, 2}];
143 pavetext->AddText(ss.str().c_str());
144 pavetext->Draw();
145 equation[i] = kTRUE;
146 }
147 c->cd(i * fDetIds.size() + j + 1);
148 // auto mg = std::make_shared<TMultiGraph>("mg", "");
149 for (auto &dataset : fReference) {
150 double xmin;
151 double xmax;
152 fFits[{idi, idj}]->GetRange(xmin, xmax);
153 xmin = fFitResults[{idj, 1}];
154 xmin = xmin < 0 ? -1 * xmin : 0;
155 fFits[{idi, idj}]->SetRange(xmin + 0.1, xmax);
156 fFits[{idi, idj}]->SetLineColorAlpha(1, 0.8);
157 fFits[{idi, idj}]->SetLineWidth(1);
158 fFits[{idi, idj}]->SetParameter(6, dataset.second->GetEnergy());
159 fFits[{idi, idj}]->DrawCopy("SAME");
160 }
161 // c -> BuildLegend();
162 c->Update();
163 ++j;
164 }
165 ++i;
166 j = 0;
167 }
168 }
169 if (!saveas.empty()) {
170 c->SaveAs(saveas.c_str());
171 }
172}
173
188Int_t TimeWalk::AddPeak(std::string dataset, GamR::Nucleus::Photon photon, Double_t time)
189{
190 if (fReference.count(dataset) == 0) {
191 int error = AddDataSet(dataset);
192 if (error)
193 return error;
194 } else {
195 fRelative[dataset].push_back({photon, time}); // peak, bg, time correction
196 }
197 return 0;
198}
199
211void TimeWalk::MakeGraphs(Int_t ntoys, Bool_t showmc)
212{
213 // Setting up
214 TH1::SetDefaultSumw2(kTRUE);
215 TCanvas *c1 = new TCanvas("canvas", "First Order Moment", 1280, 720);
216 TH2D *pkhist;
217 TH2D *bghist;
218 std::vector<Double_t> energy, tdiff, energyerr, tdifferr;
219 int colour = 0;
220
221 auto mcmean = [](TH1 *h, TRandom3 &rnd) {
222 int low = h->GetXaxis()->FindBin(h->GetMean() - 5 * h->GetStdDev());
223 int high = h->GetXaxis()->FindBin(h->GetMean() + 5 * h->GetStdDev());
224 h->GetXaxis()->SetRange(low, high);
225 // higher order iterations
226 for (int order = 2; order <= 3; ++order) {
227 low = h->GetXaxis()->FindBin(h->GetMean() - order * 5 * h->GetStdDev());
228 high = h->GetXaxis()->FindBin(h->GetMean() + order * 5 * h->GetStdDev());
229 h->GetXaxis()->SetRange(low, high);
230 }
231 for (int x = low; x <= high; ++x) {
232 Double_t error = h->GetBinError(x);
233 Double_t count = h->GetBinContent(x);
234 h->SetBinContent(x, rnd.Gaus(count, error));
235 h->SetBinError(x, 0);
236 }
237 double mean = h->GetMean();
238 double median;
239 double q = 0.5;
240 h->GetQuantiles(1, &median, &q);
241 return std::pair<double, double>(median, mean);
242 };
243
244 // Graph for each data set
245 for (const auto &dataset : fReference) {
246 auto &datafn = dataset.first;
247 auto &refpeak = dataset.second;
248 TFile datafile(datafn.c_str());
249
250 Double_t width_refpk = refpeak->GetGate()->GetWidth();
251 Double_t width_refbg = refpeak->GetGateBG()->GetWidth();
252 std::stringstream ss;
253
254 // Each pairing
255 for (auto &i : fDetIds) {
256 for (auto &j : fDetIds) {
257 if (i == j) {
258 continue;
259 }
260 energy.clear();
261 energyerr.clear();
262 tdiff.clear();
263 tdifferr.clear();
264 std::string pairing(std::to_string(i) + "_" + std::to_string(j));
265
266 datafile.GetObject(("P" + pairing).c_str(), pkhist);
267 datafile.GetObject(("B" + pairing).c_str(), bghist);
268
269 std::cout << "Generating data points for pairing " << i << " " << j << std::endl;
270 for (auto &relative : fRelative[dataset.first]) {
271 auto &photon = relative.first;
272 auto &correc = relative.second;
273 ss.str("");
274 ss << pairing << "_" << static_cast<int>(photon.GetEnergy());
275
276 TH1D *pp = photon.GetGate()->ApplyX(pkhist, ss.str().c_str());
277 TH1D *bp = photon.GetGate()->ApplyX(bghist, "bg_pk");
278 TH1D *pb = photon.GetGateBG()->ApplyX(pkhist, "pk_bg");
279 TH1D *bb = photon.GetGateBG()->ApplyX(bghist, "bg_bg");
280
281 Double_t width_pk = photon.GetGate()->GetWidth();
282 Double_t width_bg = photon.GetGateBG()->GetWidth();
283
284 // Full subtraction
285 pp->Add(pb, -1 * width_refpk / width_bg);
286 pp->Add(bp, -1 * width_refbg / width_pk);
287 pp->Add(bb, +2 * width_refbg / width_bg);
288
289 // If low stats, ignore it
290 if (pp->Integral() < 15) {
291 continue;
292 }
293
294 // Fill the data
295 // Energy
296 energy.push_back(photon.GetEnergy());
297 energyerr.push_back(0);
298
299 // HistMean(pp, ntoys)
300 c1->Clear();
301 auto mcresult = GamR::TK::MonteCarlo(mcmean, pp, ntoys, showmc, kFALSE, c1);
302
303 // Time
304 tdiff.push_back(mcresult->GetMean(1) - correc);
305 tdifferr.push_back(mcresult->GetStdDev(1));
306
307 // Print back
308 std::cout << "Energy:\t" << energy.back();
309 std::cout << "\tTime:\t" << tdiff.back() << " +/- " << tdifferr.back();
310 std::cout << std::endl;
311
312 } // got data
313
314 // Make the graph
315 fGraphs[datafn][{i, j}] = std::make_unique<TGraphErrors>(energy.size(), energy.data(), tdiff.data(),
316 energyerr.data(), tdifferr.data());
317 // Title
318 ss.str("");
319 ss << "TCG_" << pairing;
320 fGraphs[datafn][{i, j}]->SetName((ss.str() + "_" + datafn).c_str());
321
322 // Title and xis labels
323 ss << " Time Walk Calibration E_{" << i << "}=" << refpeak->GetEnergy() << " [keV]";
324 ss << ";Gamma Energy, E_{" << j << "} [keV]";
325 ss << ";$\\Delta t = t_{" << i << "} - t_{" << j << "}$ [ps]";
326 fGraphs[datafn][{i, j}]->SetTitle(ss.str().c_str());
327
328 // Better Markers
329 colour += 1;
330 fGraphs[datafn][{i, j}]->SetLineColorAlpha((colour % 9) + 40, 0.9);
331 fGraphs[datafn][{i, j}]->SetLineWidth(1);
332 fGraphs[datafn][{i, j}]->SetMarkerColorAlpha((colour % 9) + 40, 0.9);
333 fGraphs[datafn][{i, j}]->SetMarkerColor(kOpenSquare);
334 fGraphs[datafn][{i, j}]->SetMarkerSize(2);
335
336 } // det j
337 } // det i
338 datafile.Close();
339 } // data set
340
341 this->Draw(c1, kFALSE);
342
343 return;
344} // TimeWalk::MakeGraphs
345
346class TimeWalk::GlobalChiSquare {
347private:
348 static TimeWalk *p;
349
350public:
351 GlobalChiSquare(TimeWalk *parent) { p = parent; }
352
353 static void minuit(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
354 {
355 f = 0;
356 for (const auto &dataset : p->fReference) {
357 std::string data = dataset.first;
358
359 for (size_t i = 0; i < p->fDetIds.size(); ++i) {
360 for (size_t j = 0; j < p->fDetIds.size(); ++j) {
361 if (i == j) {
362 continue;
363 }
364 auto &idi = p->fDetIds[i];
365 auto &idj = p->fDetIds[j];
366 p->fFits[{idi, idj}]->SetParameter(0, par[j * 3 + 0]); // A_j
367 p->fFits[{idi, idj}]->SetParameter(1, par[j * 3 + 1]); // B_j
368 p->fFits[{idi, idj}]->SetParameter(2, par[j * 3 + 2]); // f_j
369
370 p->fFits[{idi, idj}]->SetParameter(3, par[i * 3 + 0]); // A_i
371 p->fFits[{idi, idj}]->SetParameter(4, par[i * 3 + 1]); // B_i
372 p->fFits[{idi, idj}]->SetParameter(5, par[i * 3 + 2]); // f_i
373
374 p->fFits[{idi, idj}]->SetParameter(6, dataset.second->GetEnergy());
375
376 f += p->fGraphs[data][{idi, idj}]->Chisquare(p->fFits[{idi, idj}].get());
377 }
378 }
379 }
380 } // end operator()
381};
382TimeWalk *TimeWalk::GlobalChiSquare::p = nullptr;
383
389void TimeWalk::Fit()
390{
391 Int_t npars = 3 * fDetIds.size();
392 auto gMinuit = new TMinuit(npars);
393
394 GlobalChiSquare func(this);
395 gMinuit->SetFCN(func.minuit);
396
397 Double_t arglist[10];
398 arglist[0] = 1; // ???
399 Int_t ierflg = 0; // some error flag ???
400
401 gMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
402 arglist[0] = 2;
403 gMinuit->mnexcm("SET STR", arglist, 1, ierflg);
404
405 /* Ye olde stringstream */
406 std::stringstream ss;
407
408 TRandom3 rnd;
409
410 /* SET INITIAL VALUES */
411 Int_t parnum = 0;
412 ss.str("");
413 for (auto &detid : fDetIds) {
414 /* Curvature: A */
415 ss << "a" << detid;
416 gMinuit->mnparm(parnum, ss.str(), 2500, 1, -1e4, 1e4, ierflg);
417 ss.str("");
418 ++parnum;
419
420 /* Bend point: B */
421 ss << "b" << detid;
422 gMinuit->mnparm(parnum, ss.str(), -10, 0.01, -1e4, 1e4, ierflg);
423 ss.str("");
424 ++parnum;
425
426 /* Constant flight time: F */
427 ss << "f" << detid;
428 gMinuit->mnparm(parnum, ss.str(), 0, 0.1, -1e3, 1e3, ierflg);
429 ss.str("");
430 ++parnum;
431
432 /* Linear correction: L */
433 // ss << "l" << detid;
434 // gMinuit -> mnparm(parnum, ss.str(), 0, 0.0001, -0.01, 0.01, ierflg);
435 // ss.str(""); ++parnum;
436 }
437 gMinuit->FixParameter(2); // Fix Flight of Detector 1
438
439 /* MINIMIZE */
440 arglist[0] = 100000; // Max Iterations
441 arglist[1] = 1; // Number of sigma max step
442 gMinuit->mnexcm("MINIGRAD", arglist, 2, ierflg);
443 gMinuit->mnimpr(); // Improve fit by searching for other minima
444 gMinuit->mnexcm("HESSE", arglist, 2, ierflg);
445 gMinuit->mnexcm("MINOS", arglist, 2, ierflg);
446
447 TCanvas *c1 = new TCanvas("cTCAL", "Time Walk Calibration Fit Results", 1280, 720);
448 this->Draw(c1, kTRUE);
449
450 // Print results
451 std::cout << "\nFit Results\n";
452 double param, err;
453 for (size_t i = 0; i < fDetIds.size(); ++i) {
454 std::cout << std::endl << "Detector " << fDetIds[i] << std::endl;
455 // Curvature
456 gMinuit->GetParameter(i * 3 + 0, param, err);
457 fFitResults[{fDetIds[i], 0}] = param;
458 std::cout << "(([1]+x)-((" << param << "+(y*" << err << "))";
459
460 // Bend point
461 gMinuit->GetParameter(i * 3 + 1, param, err);
462 fFitResults[{fDetIds[i], 1}] = param;
463 std::cout << "/sqrt([0]+(" << param << "+(y*" << err << ")))";
464
465 // Linear correction
466 // gMinuit -> GetParameter(i*4+3, param, err);
467 // std::cout << "+[0]*(" << param << "+(y*" << err << "))";
468
469 // Flight
470 gMinuit->GetParameter(i * 3 + 2, param, err);
471 fFitResults[{fDetIds[i], 2}] = param;
472 std::cout << "+(" << param << "+(y*" << err << "))))";
473 std::cout << std::endl;
474 }
475
476 std::cout << "All Done" << std::endl;
477 return;
478} // TimeWalk::Fit
479
480} // namespace Processor
481} // namespace GamR
bp
GamR::TK::BPeak * bp(TCanvas *canvas, Option_t *foption, Option_t *option)
Definition Abbreviations.cc:176
ss
void ss(std::vector< int > indexes, TCanvas *canvas, Option_t *option)
Definition Abbreviations.cc:205
GamR::Processor::TimeWalk::GlobalChiSquare::minuit
static void minuit(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
Definition TimeWalk.cc:353
GamR::Processor::TimeWalk::Draw
void Draw(TCanvas *c, Bool_t fits=kFALSE, std::string saveas="")
Draws Time Walk graphs.
Definition TimeWalk.cc:83
GamR::Processor::TimeWalk::MakeGraphs
void MakeGraphs(Int_t ntoys, Bool_t showmc=kTRUE)
Iterating through calibration peaks for each detector pair: calculate the time difference spectrum (w...
Definition TimeWalk.cc:211
GamR::Processor::TimeWalk::TimeWalk
TimeWalk(std::vector< Int_t > ids)
TimeWalk Calibration Processor.
Definition TimeWalk.cc:26
GamR::Processor::TimeWalk::AddPeak
Int_t AddPeak(std::string dataset, GamR::Nucleus::Photon photon, Double_t time)
Add a calibration peak to be used to create a data point in the time walk curve.
Definition TimeWalk.cc:188
GamR::Processor::TimeWalk::Fit
void Fit()
Conducts a global fit to solve for an absolute time calibration for each detector.
Definition TimeWalk.cc:389
GamR::TK::MonteCarlo
std::shared_ptr< TH2D > MonteCarlo(std::function< std::pair< double, double >(TH1 *, TRandom3 &)> func, TH1 *system, size_t ntoys, Bool_t draw, const Bool_t weight, TCanvas *c1)
Definition MonteCarlo.cc:9
GamR
Definition Gain.cc:19