GRSISort "v4.1.1.0"
An extension of the ROOT analysis Framework
Loading...
Searching...
No Matches
TCalibrationGraph.cxx
Go to the documentation of this file.
1#include "TCalibrationGraph.h"
2
3#include <iostream>
4#include <algorithm>
5
6#include "TMath.h"
7#include "TPad.h"
8#include "TF1.h"
9#include "TBuffer.h"
10#include "TH1F.h"
11
13{
14 if(TCalibrationGraphSet::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << __PRETTY_FUNCTION__ << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
15 Int_t px = gPad->GetEventX();
16 Int_t py = gPad->GetEventY();
17 if(TCalibrationGraphSet::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << "px, py " << px << ", " << py << " on gPad " << gPad->GetName() << std::endl; }
18 if(fIsResidual) { return fParent->RemoveResidualPoint(px, py); }
19 return fParent->RemovePoint(px, py);
20}
21
22#if ROOT_VERSION_CODE < ROOT_VERSION(6, 26, 0)
23void TCalibrationGraph::Scale(const double& scale)
24{
25 Double_t* y = GetY();
26 Double_t* ey = GetEY();
27
28 for(int i = 0; i < GetN(); ++i) {
29 y[i] *= scale;
30 ey[i] *= scale;
31 }
32}
33#endif
34
36
37TCalibrationGraphSet::TCalibrationGraphSet(TGraphErrors* graph, const std::string& label)
38 : fTotalGraph(new TGraphErrors), fTotalResidualGraph(new TGraphErrors)
39{
40 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << __PRETTY_FUNCTION__ << " fTotalGraph " << fTotalGraph << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
41 if(graph != nullptr) {
42 Add(graph, label);
44 }
45}
46
47TCalibrationGraphSet::TCalibrationGraphSet(std::string xAxisLabel, std::string yAxisLabel)
48 : fTotalGraph(new TGraphErrors), fTotalResidualGraph(new TGraphErrors), fXAxisLabel(std::move(xAxisLabel)), fYAxisLabel(std::move(yAxisLabel))
49{
50 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << __PRETTY_FUNCTION__ << " fTotalGraph " << fTotalGraph << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
51}
52
59
60int TCalibrationGraphSet::Add(TGraphErrors* graph, const std::string& label)
61{
63 std::cout << __PRETTY_FUNCTION__ << ", fTotalGraph " << fTotalGraph << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
64 Print();
65 }
66 if(graph->GetN() == 0) {
67 std::cout << __PRETTY_FUNCTION__ << ": graph \"" << graph->GetName() << "\", label \"" << label << "\" is empty, not adding it" << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
68 return -1;
69 }
70
71 graph->GetListOfFunctions()->Clear(); // we don't want to include any fits
72
73 // get points and error bars from our calibration
74 double* x = fTotalGraph->GetX();
75 double* y = fTotalGraph->GetY();
76 double* ex = fTotalGraph->GetEX();
77 double* ey = fTotalGraph->GetEY();
78
79 // get points and error bars from graph
80 double* rhsX = graph->GetX();
81 double* rhsY = graph->GetY();
82 double* rhsEX = graph->GetEX();
83 double* rhsEY = graph->GetEY();
84
85 // create one vector with x, y, ex, ey, index of graph, and index of point that we can use to sort the data
86 // TODO: create a (private?) enum to reference to x, y, ex, ey, graph index, and point index
87 std::vector<std::tuple<double, double, double, double, size_t, size_t>> data(fTotalGraph->GetN() + graph->GetN());
88 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "Filling vector of size " << data.size() << " with " << fTotalGraph->GetN() << " and " << graph->GetN() << " entries" << std::endl; }
89 for(int i = 0; i < fTotalGraph->GetN(); ++i) {
90 data[i] = std::make_tuple(x[i], y[i], ex[i], ey[i], fGraphIndex[i], fPointIndex[i]);
91 }
92 for(int i = 0; i < graph->GetN(); ++i) {
93 data[fTotalGraph->GetN() + i] = std::make_tuple(rhsX[i], rhsY[i], rhsEX[i], rhsEY[i], fGraphs.size(), i);
94 }
95
96 std::sort(data.begin(), data.end());
97
98 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "sorted vector, setting graph sizes" << std::endl; }
99
100 fTotalGraph->Set(static_cast<Int_t>(data.size()));
101 fTotalResidualGraph->Set(static_cast<Int_t>(data.size()));
102 fGraphIndex.resize(data.size());
103 fPointIndex.resize(data.size());
104
105 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "Filling fTotalGraph, fGraphIndex, and fPointIndex with " << data.size() << " points" << std::endl; }
106 for(size_t i = 0; i < data.size(); ++i) {
107 fTotalGraph->SetPoint(static_cast<Int_t>(i), std::get<0>(data[i]), std::get<1>(data[i]));
108 fTotalGraph->SetPointError(static_cast<Int_t>(i), std::get<2>(data[i]), std::get<3>(data[i]));
109 fGraphIndex[i] = std::get<4>(data[i]);
110 fPointIndex[i] = std::get<5>(data[i]);
111 }
112 fMinimumX = std::get<0>(data[0]);
113 fMinimumY = std::get<1>(data[0]);
114 fMaximumX = std::get<0>(data.back());
115 fMaximumY = std::get<1>(data.back());
116 // doesn't really make sense to calculate the residual here, as we don't have a fit of all the data yet
117 fResidualSet = false;
118
119 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "Adding new calibration graph and label to vectors" << std::endl; }
120 // add graph and label to our vectors
121 fGraphs.emplace_back(this, graph);
122 fResidualGraphs.emplace_back(this, 0, true);
123 fLabel.push_back(label);
124 // set initial color
125 fGraphs.back().SetLineColor(static_cast<Color_t>(fGraphs.size()));
126 fGraphs.back().SetMarkerColor(static_cast<Color_t>(fGraphs.size()));
127 fGraphs.back().SetMarkerStyle(static_cast<Style_t>(fGraphs.size()));
128 fResidualGraphs.back().SetLineColor(static_cast<Color_t>(fResidualGraphs.size()));
129 fResidualGraphs.back().SetMarkerColor(static_cast<Color_t>(fResidualGraphs.size()));
130 fResidualGraphs.back().SetMarkerStyle(static_cast<Style_t>(fResidualGraphs.size()));
131
133 std::cout << "done" << std::endl;
134 Print();
135 }
136 return static_cast<int>(fGraphs.size() - 1);
137}
138
140{
141 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << __PRETTY_FUNCTION__ << " gpad " << gPad->GetName() << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
142 TF1* calibration = FitFunction();
143 if(calibration != nullptr && (!fResidualSet || force)) {
145 std::cout << "calibration " << calibration << ", fResidualSet " << (fResidualSet ? "true" : "false") << ", force " << (force ? "true" : "false") << ", calibration: " << std::endl;
146 calibration->Print();
147 }
148 double* x = fTotalGraph->GetX();
149 double* y = fTotalGraph->GetY();
150 double* ex = fTotalGraph->GetEX();
151 double* ey = fTotalGraph->GetEY();
152 fTotalResidualGraph->Set(fTotalGraph->GetN());
153 for(int i = 0; i < fTotalGraph->GetN(); ++i) {
154 fTotalResidualGraph->SetPoint(i, y[i] - calibration->Eval(x[i]), y[i]);
155 fTotalResidualGraph->SetPointError(i, TMath::Sqrt(TMath::Power(ey[i], 2) + TMath::Power(ex[i] * calibration->Derivative(x[i]), 2)), ey[i]);
156 }
158 std::cout << "Done calculating total residual graph with " << fTotalResidualGraph->GetN() << " points" << std::endl;
159 }
160 for(size_t g = 0; g < fGraphs.size(); ++g) {
161 x = fGraphs[g].GetX();
162 y = fGraphs[g].GetY();
163 ex = fGraphs[g].GetEX();
164 ey = fGraphs[g].GetEY();
165 fResidualGraphs[g].Set(fGraphs[g].GetN());
166 for(int i = 0; i < fGraphs[g].GetN(); ++i) {
167 fResidualGraphs[g].SetPoint(i, y[i] - calibration->Eval(x[i]), y[i]);
168 fResidualGraphs[g].SetPointError(i, TMath::Sqrt(TMath::Power(ey[i], 2) + TMath::Power(ex[i] * calibration->Derivative(x[i]), 2)), ey[i]);
169 }
170 }
172 std::cout << "Done calculating all " << fGraphs.size() << " individual residual graphs" << std::endl;
173 }
174 fResidualSet = true;
175 auto* mother = gPad->GetMother();
176 int pad = 0;
177 while(mother->GetPad(pad) != nullptr) {
178 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "pad " << pad << " = " << std::flush << mother->GetPad(pad) << " = " << std::flush << mother->GetPad(pad)->GetName() << ": modifying, " << std::flush; }
179 mother->GetPad(pad)->Modified();
180 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "updating, " << std::flush; }
181 mother->GetPad(pad)->Update();
182 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "changing to pad, " << std::flush; }
183 mother->cd(pad);
184 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "done!" << std::endl; }
185 pad++;
186 }
187 } else {
188 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << __PRETTY_FUNCTION__ << ": didn't find calibration (" << calibration << "), or the residual was already set (" << (fResidualSet ? "true" : "false") << ") and we don't force it (" << (force ? "true" : "false") << ")" << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
189 if(calibration == nullptr) { fResidualSet = false; }
190 }
192 return fResidualSet;
193}
194
196{
197 if(fTotalGraph == nullptr) {
198 std::cerr << "Trying to set axis title to \"" << title << "\" failed because no total graph is present for this title to be applied to!" << std::endl;
199 return;
200 }
201 fTotalGraph->SetTitle(Form("%s;%s", fTotalGraph->GetTitle(), title));
202}
203
204void TCalibrationGraphSet::DrawCalibration(Option_t* opt, TLegend* legend)
205{
206 TString options = opt;
207 options.ToLower();
208 if(options.Contains("same")) {
209 options.Remove(options.Index("same"), 4);
210 opt = options.Data();
211 } else {
212 options.Append("a");
213 }
214 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << __PRETTY_FUNCTION__ << " drawing total graph with option \"" << options.Data() << "\" on gPad " << gPad->GetName() << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
215 fTotalGraph->Draw(options.Data());
216
217 if(fTotalGraph->GetHistogram() != nullptr) {
218 fTotalGraph->GetHistogram()->GetXaxis()->CenterTitle();
219 fTotalGraph->GetHistogram()->GetXaxis()->SetTitle(fXAxisLabel.c_str());
220 fTotalGraph->GetHistogram()->GetYaxis()->CenterTitle();
221 fTotalGraph->GetHistogram()->GetYaxis()->SetTitle(fYAxisLabel.c_str());
222 }
223
224 for(size_t i = 0; i < fGraphs.size(); ++i) {
225 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << __PRETTY_FUNCTION__ << " drawing " << i << ". graph with option \"" << opt << "\", marker color " << fGraphs[i].GetMarkerColor() << " on gPad " << gPad->GetName() << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
226 fGraphs[i].Draw(opt);
227 if(legend != nullptr) {
228 legend->AddEntry(&(fGraphs[i]), fLabel[i].c_str());
229 }
230 }
231 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << __PRETTY_FUNCTION__ << " drawing legend " << legend << " on gPad " << gPad->GetName() << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
232 if(legend != nullptr) {
233 legend->Draw();
234 }
235}
236
238{
239 if(fZeroResidual == nullptr) { return; }
240
241 auto* hist = fTotalResidualGraph->GetHistogram();
242 double minY = hist->GetYaxis()->GetXmin();
243 double maxY = hist->GetYaxis()->GetXmax();
244 fZeroResidual->SetY1(minY);
245 fZeroResidual->SetY2(maxY);
247 std::cout << __PRETTY_FUNCTION__ << " changing line at zero to go from " << minY << " to " << maxY << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
248 }
249}
250
251void TCalibrationGraphSet::DrawResidual(Option_t* opt, TLegend* legend)
252{
253 TString options = opt;
254 options.ToLower();
255 options.Append("a");
256 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << __PRETTY_FUNCTION__ << " drawing total residual graph with option \"" << options.Data() << "\" on gPad " << gPad->GetName() << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
257 fTotalResidualGraph->Draw(options.Data());
258 auto* hist = fTotalResidualGraph->GetHistogram();
259 if(hist != nullptr) {
260 hist->GetXaxis()->SetLabelSize(0.06);
261 // check if a reference line at zero has been requested
262 if(options.Contains("r0")) {
263 // find the y-range to plot the line on
264 double minY = hist->GetYaxis()->GetXmin();
265 double maxY = hist->GetYaxis()->GetXmax();
267 std::cout << __PRETTY_FUNCTION__ << " drawing line at zero from " << minY << " to " << maxY << " as requested (\"" << options.Data() << "\") " << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
268 }
269 if(fZeroResidual == nullptr) {
270 fZeroResidual = new TLine(0., minY, 0., maxY);
271 } else {
272 fZeroResidual->SetY1(minY);
273 fZeroResidual->SetY2(maxY);
274 }
275 fZeroResidual->Draw("same");
276 }
277 } else {
278 std::cout << "Failed to get histogram for graph:" << std::endl;
279 fTotalResidualGraph->Print();
280 }
281
282 options = opt;
283 options.ReplaceAll("r0", "");
284 for(size_t i = 0; i < fResidualGraphs.size(); ++i) {
285 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << __PRETTY_FUNCTION__ << " drawing " << i << ". residual graph with option \"" << opt << "\", marker color " << fResidualGraphs[i].GetMarkerColor() << " on gPad " << gPad->GetName() << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
286 fResidualGraphs[i].Draw(options.Data());
287 if(legend != nullptr) {
288 legend->AddEntry(&(fResidualGraphs[i]), fLabel[i].c_str());
289 }
290 }
291 if(legend != nullptr) {
292 legend->Draw();
293 }
294}
295
296Int_t TCalibrationGraphSet::RemovePoint(const Int_t& px, const Int_t& py)
297{
298 /// This function calls RemovePoint on the total graph.
299
301 std::cout << __PRETTY_FUNCTION__ << ": point " << px << ", " << py << "; gPad " << gPad->GetName() << ": " << gPad->AbsPixeltoX(px) << ", " << gPad->AbsPixeltoY(py) << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
302 }
303
304 return RemovePoint(fTotalGraph, px, py);
305}
306
307Int_t TCalibrationGraphSet::RemoveResidualPoint(const Int_t& px, const Int_t& py)
308{
309 /// This function calls RemovePoint on the total graph.
310
312 std::cout << __PRETTY_FUNCTION__ << ": point " << px << ", " << py << "; gPad " << gPad->GetName() << ": " << gPad->AbsPixeltoX(px) << ", " << gPad->AbsPixeltoY(py) << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
313 }
314
315 return RemovePoint(fTotalResidualGraph, px, py);
316}
317
318Int_t TCalibrationGraphSet::RemovePoint(TGraphErrors* graph, const Int_t& px, const Int_t& py)
319{
320 /// This function is primarily a copy of TGraph::RemovePoint with some added bits to remove a point that has been selected in the residual graph from it and the corresponding point from the calibration graph and the total graphs
321
323
324 // localize point to be deleted
325 Int_t point = -2;
326 Int_t i = 0;
327 // start with a small window (in case the mouse is very close to one point)
328 double* x = graph->GetX();
329 double* y = graph->GetY();
330 for(i = 0; i < graph->GetN(); i++) {
331 Int_t dpx = px - gPad->XtoAbsPixel(gPad->XtoPad(x[i]));
332 Int_t dpy = py - gPad->YtoAbsPixel(gPad->YtoPad(y[i]));
333 // TODO replace 100 with member variable?
334 if(dpx * dpx + dpy * dpy < 100) {
336 std::cout << i << ": dpx = " << dpx << " = " << px << " - " << gPad->XtoAbsPixel(gPad->XtoPad(x[i])) << ", dpy = " << dpy << " = " << py << " - " << gPad->YtoAbsPixel(gPad->YtoPad(y[i])) << " this is the point we're looking for at " << x[i] << ", " << y[i] << std::endl;
337 }
338 point = i;
339 break;
340 }
342 std::cout << i << ": dpx = " << dpx << " = " << px << " - " << gPad->XtoAbsPixel(gPad->XtoPad(x[i])) << ", dpy = " << dpy << " = " << py << " - " << gPad->YtoAbsPixel(gPad->YtoPad(y[i])) << " not the point we're looking for" << std::endl;
343 }
344 }
345
346 if(point < 0) {
347 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << "Failed to find point close to " << px << ", " << py << std::endl; }
349 std::cout << __PRETTY_FUNCTION__ << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
350 Print();
351 }
352 return point;
353 }
355 Print();
356 }
357
358 // now that we have the point we want to delete, remove it from the total graphs
359 fTotalGraph->RemovePoint(point);
360 if(fTotalResidualGraph->RemovePoint(point) < 0) {
361 // we failed to remove the point in the residual, so we assume it's out of whack
362 fResidualSet = false;
363 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << point << " didn't removed residual point" << std::endl; }
365 std::cout << point << " removed residual point" << std::endl;
366 }
367
368 // need to find which of the graphs we have to remove this point from -> use fGraphIndex[point]
369 // and also which point this is of the graph -> use fPointIndex[point]
370 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << point << " - " << fGraphIndex[point] << ", " << fPointIndex[point] << std::endl; }
371 if(fGraphs[fGraphIndex[point]].RemovePoint(fPointIndex[point]) == -1 || fResidualGraphs[fGraphIndex[point]].RemovePoint(fPointIndex[point]) == -1) {
372 std::cout << "point " << point << " out of range?" << std::endl;
373 }
374
375 // now we need to update the indices by removing this one, and updating all that come after it
376 auto oldGraph = fGraphIndex[point];
377 auto oldPoint = fPointIndex[point];
378 fGraphIndex.erase(fGraphIndex.begin() + point);
379 fPointIndex.erase(fPointIndex.begin() + point);
380 for(size_t p = 0; p < fGraphIndex.size(); ++p) {
381 if(fGraphIndex[p] == oldGraph && fPointIndex[p] >= oldPoint) {
382 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << fGraphIndex[p] << " == " << oldGraph << " && " << fPointIndex[p] << " >= " << oldPoint << ": decrementing index" << std::endl; }
383 --fPointIndex[p];
385 std::cout << fGraphIndex[p] << " != " << oldGraph << " || " << fPointIndex[p] << " < " << oldPoint << ": decrementing index" << std::endl;
386 }
387 }
388
389 // update the graphics
390 auto* mother = gPad->GetMother();
391 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "Got mother pad " << mother->GetName() << " from pad " << gPad->GetName() << std::endl; }
392 int pad = 0;
393 while(mother->GetPad(pad) != nullptr) {
394 mother->GetPad(pad)->Modified(); // one version also used Update and cd(pad)?
395 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "Modified pad " << pad << " = " << mother->GetPad(pad)->GetName() << std::endl; }
396 pad++;
397 }
399
400 // emit signal that this point has been removed
401 std::array<Longptr_t, 2> args = {static_cast<int64_t>(oldGraph), static_cast<int64_t>(oldPoint)};
402 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "Emitting RemovePoint(Int_t, Int_t) with " << args[0] << ", " << args[1] << " - " << args.data() << std::endl; }
403 Emit("RemovePoint(Int_t, Int_t)", args.data());
404
405 return point;
406}
407
409{
410 /// This function removes the indicated point from the total graph.
411
413
414 // now that we have the point we want to delete, remove it from the total graphs
415 fTotalGraph->RemovePoint(point);
416 if(fTotalResidualGraph->RemovePoint(point) < 0) {
417 // we failed to remove the point in the residual, so we assume it's out of whack
418 fResidualSet = false;
419 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << point << " didn't removed residual point" << std::endl; }
421 std::cout << point << " removed residual point" << std::endl;
422 }
423
424 // need to find which of the graphs we have to remove this point from -> use fGraphIndex[point]
425 // and also which point this is of the graph -> use fPointIndex[point]
426 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "Going to remove " << point << " - " << fGraphIndex[point] << ", " << fPointIndex[point] << std::endl; }
427 if(fGraphs[fGraphIndex[point]].RemovePoint(fPointIndex[point]) == -1 || fResidualGraphs[fGraphIndex[point]].RemovePoint(fPointIndex[point]) == -1) {
428 std::cout << "point " << point << " out of range?" << std::endl;
429 }
430
431 // now we need to update the indices by removing this one, and updating all that come after it
432 auto oldGraph = fGraphIndex[point];
433 auto oldPoint = fPointIndex[point];
434 fGraphIndex.erase(fGraphIndex.begin() + point);
435 fPointIndex.erase(fPointIndex.begin() + point);
436 for(size_t p = 0; p < fGraphIndex.size(); ++p) {
437 if(fGraphIndex[p] == oldGraph && fPointIndex[p] >= oldPoint) {
438 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << fGraphIndex[p] << " == " << oldGraph << " && " << fPointIndex[p] << " >= " << oldPoint << ": decrementing index" << std::endl; }
439 --fPointIndex[p];
441 std::cout << fGraphIndex[p] << " != " << oldGraph << " || " << fPointIndex[p] << " < " << oldPoint << ": decrementing index" << std::endl;
442 }
443 }
444
445 // update the graphics
446 auto* mother = gPad->GetMother();
447 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "Got mother pad " << mother->GetName() << " from pad " << gPad->GetName() << std::endl; }
448 int pad = 0;
449 while(mother->GetPad(pad) != nullptr) {
450 mother->GetPad(pad)->Modified(); // one version also used Update and cd(pad)?
451 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "Modified pad " << pad << " = " << mother->GetPad(pad)->GetName() << std::endl; }
452 pad++;
453 }
455
456 // emit signal that this point has been removed
457 std::array<Longptr_t, 2> args = {static_cast<int64_t>(oldGraph), static_cast<int64_t>(oldPoint)};
458 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "Emitting RemovePoint(Int_t, Int_t) with " << args[0] << ", " << args[1] << " - " << args.data() << std::endl; }
459 Emit("RemovePoint(Int_t, Int_t)", args.data());
460}
461
463{
464 for(auto& graph : fGraphs) {
465 graph.Sort();
466 }
467 for(auto& graph : fResidualGraphs) {
468 // residuals plot energy vs. residual, so we want to sort by y, not x!
469 graph.Sort(&TGraph::CompareY);
470 }
471}
472
473void TCalibrationGraphSet::Scale(bool useAllPrevious)
474{
475 /// Scale all graphs to fit each other (based on the first "previous" graph found or just the first graph).
476 /// If no overlap is being found between the graph that is being scaled and the first graph (or all graphs before this one),
477 /// the current graph isn't being scaled and we continue with the next graph.
479 std::cout << __PRETTY_FUNCTION__ << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
480 Print();
481 }
482 Sort();
483 for(size_t g = 1; g < fGraphs.size(); ++g) {
484 auto& graph = fGraphs[g];
485 if(graph.GetN() == 0) {
486 std::cout << "No entries in " << g << ". graph" << std::endl;
487 }
488 double* x = graph.GetX();
489 double* y = graph.GetY();
490 // loop over all other graphs before this one and use the first one with an overlap (should this be extended to use all possible ones before this one?)
491 size_t g2 = 0;
492 for(g2 = 0; (useAllPrevious ? g2 < g : g2 < 1); ++g2) {
493 double minRef = fGraphs[g2].GetX()[0];
494 double maxRef = fGraphs[g2].GetX()[fGraphs[g2].GetN() - 1];
496 std::cout << "Checking overlap between " << g2 << ". graph (" << fGraphs[g2].GetN() << ": " << minRef << " - " << maxRef << ") and " << g << ". graph (" << graph.GetN() << std::flush << ": " << x[0] << " - " << x[graph.GetN() - 1] << ")" << std::endl;
497 }
498 if(maxRef < x[0] || x[graph.GetN() - 1] < minRef) {
499 // no overlap between the two graphs, for now we just skip this one, but we could try and compare it to all the other ones?
500 std::cout << "No overlap between " << g2 << ". graph (" << fGraphs[g2].GetN() << ": " << minRef << " - " << maxRef << ") and " << g << ". graph (" << graph.GetN() << ": " << x[0] << " - " << x[graph.GetN() - 1] << ")" << std::endl;
501 continue;
502 }
503 // we have an overlap, so we calculate the scaling factor for each point and take the average (maybe should add some weight from the errors bars)
504 int count = 0;
505 double sum = 0.;
506 for(int p = 0; p < graph.GetN(); ++p) {
507 if(minRef < x[p] && x[p] < maxRef) {
508 sum += fGraphs[g2].Eval(x[p]) / y[p];
509 ++count;
510 if(fVerboseLevel > EVerbosity::kLoops) { std::cout << g << ", " << p << ": " << count << " - " << sum << ", " << fGraphs[g2].Eval(x[p]) / y[p] << std::endl; }
511 }
512 }
513 sum /= count;
514 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << g << ": scaling with " << sum << std::endl; }
515 graph.Scale(sum);
516 break;
517 }
518 if(g2 == g && fVerboseLevel > EVerbosity::kQuiet) {
519 std::cout << "No overlap(s) between 0. to " << g2 - 1 << ". graph and " << g << ". graph (" << graph.GetN() << ": " << x[0] << " - " << x[graph.GetN() - 1] << "), not scaling this one!" << std::endl;
520 }
521 }
524}
525
527{
528 if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << __PRETTY_FUNCTION__ << std::endl; } // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
529 if(fGraphs.empty()) {
530 std::cerr << "No graphs added yet, makes no sense to reset total graph?" << std::endl;
531 return;
532 }
533
534 size_t newSize = 0;
535 for(auto& graph : fGraphs) {
536 newSize += graph.GetN();
537 }
538 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "creating graph with " << newSize << " points" << std::endl; }
539 // create one vector with x, y, ex, ey, index of graph, and index of point that we can use to sort the data
540 std::vector<std::tuple<double, double, double, double, size_t, size_t>> data(newSize);
541 size_t counter = 0;
542 for(size_t i = 0; i < fGraphs.size(); ++i) {
543 // get points and error bars from graph
544 double* x = fGraphs[i].GetX();
545 double* y = fGraphs[i].GetY();
546 double* eX = fGraphs[i].GetEX();
547 double* eY = fGraphs[i].GetEY();
548 for(int p = 0; p < fGraphs[i].GetN(); ++p, ++counter) {
549 if(fVerboseLevel > EVerbosity::kLoops) { std::cout << "filling point " << counter << " of vector of size " << newSize << std::endl; }
550 data[counter] = std::make_tuple(x[p], y[p], eX[p], eY[p], i, p);
551 }
552 }
553
554 std::sort(data.begin(), data.end());
555
556 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "sorted vector, setting graph sizes" << std::endl; }
557
558 fTotalGraph->Set(data.size());
559 fTotalResidualGraph->Set(data.size());
560 fGraphIndex.resize(data.size());
561 fPointIndex.resize(data.size());
562
563 if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "Filling fTotalGraph, fGraphIndex, and fPointIndex with " << data.size() << " points" << std::endl; }
564 for(size_t i = 0; i < data.size(); ++i) {
565 fTotalGraph->SetPoint(i, std::get<0>(data[i]), std::get<1>(data[i]));
566 fTotalGraph->SetPointError(i, std::get<2>(data[i]), std::get<3>(data[i]));
567 fGraphIndex[i] = std::get<4>(data[i]);
568 fPointIndex[i] = std::get<5>(data[i]);
569 }
570 // doesn't really make sense to calculate the residual here, as we don't have a fit of all the data yet
571 fResidualSet = false;
572 // set the new minima and maxima (always assuming the first and last points are minima and maxima, respectively)
573 fMinimumX = std::get<0>(data[0]);
574 fMinimumY = std::get<1>(data[0]);
575 fMaximumX = std::get<0>(data.back());
576 fMaximumY = std::get<1>(data.back());
577}
578
579void TCalibrationGraphSet::Print(Option_t* opt) const
580{
582 std::cout << __PRETTY_FUNCTION__ << ", fTotalGraph " << fTotalGraph << std::endl; // NOLINT(cppcoreguidelines-pro-type-const-cast, cppcoreguidelines-pro-bounds-array-to-pointer-decay)
583 }
584
585 std::cout << "TCalibrationGraphSet " << this << " - " << GetName() << ": " << fGraphs.size() << " calibration graphs, " << fResidualGraphs.size() << " residual graphs, " << fLabel.size() << " labels, ";
586 if(fTotalGraph != nullptr) {
587 std::cout << fTotalGraph->GetN() << " calibration points, and ";
588 } else {
589 std::cout << " no calibration points, and ";
590 }
591 if(fTotalResidualGraph != nullptr) {
592 std::cout << fTotalResidualGraph->GetN() << " residual points" << std::endl;
593 } else {
594 std::cout << " no residual points" << std::endl;
595 }
596 TString options = opt;
597 bool errors = options.Contains("e", TString::ECaseCompare::kIgnoreCase);
598 for(const auto& g : fGraphs) {
599 double* x = g.GetX();
600 double* y = g.GetY();
601 double* ex = g.GetEX();
602 double* ey = g.GetEY();
603 for(int p = 0; p < g.GetN(); ++p) {
604 if(errors) {
605 std::cout << p << " - " << x[p] << "(" << ex[p] << "), " << y[p] << "(" << ey[p] << "); ";
606 } else {
607 std::cout << p << " - " << x[p] << ", " << y[p] << "; ";
608 }
609 }
610 std::cout << std::endl;
611 }
612 std::cout << fGraphIndex.size() << " graph indices: ";
613 for(const auto& i : fGraphIndex) { std::cout << std::setw(3) << i << " "; }
614 std::cout << std::endl;
615 std::cout << fPointIndex.size() << " point indices: ";
616 for(const auto& i : fPointIndex) { std::cout << std::setw(3) << i << " "; }
617 std::cout << std::endl;
618 std::cout << "---- total graph ----" << std::endl;
619 double* x = fTotalGraph->GetX();
620 double* y = fTotalGraph->GetY();
621 for(int p = 0; p < fTotalGraph->GetN(); ++p) {
622 std::cout << p << " - " << x[p] << ", " << y[p] << "; ";
623 }
624 std::cout << std::endl;
625 std::cout << "---------------------" << std::endl;
626}
627
628void TCalibrationGraphSet::Clear(Option_t* option)
629{
630 fGraphs.clear();
631 fResidualGraphs.clear();
632 fLabel.clear();
633 fGraphIndex.clear();
634 fPointIndex.clear();
635 fResidualSet = false;
636 fMinimumX = 0.;
637 fMaximumX = 0.;
638 fMinimumY = 0.;
639 fMaximumY = 0.;
640 TNamed::Clear(option);
641 fTotalGraph->Set(0);
642 fTotalResidualGraph->Set(0);
643}
EVerbosity
Definition Globals.h:143
TH1D * hist
Definition UserFillObj.h:3
TCalibrationGraphSet * fParent
pointer to the set this graph belongs to
bool fIsResidual
flag to indicate that this graph is for residuals
Int_t RemovePoint() override
void ResetTotalGraph()
reset the total graph and add the individual ones again (used e.g. after scaling of individual graphs...
void DrawResidual(Option_t *opt="", TLegend *legend=nullptr)
std::vector< size_t > fGraphIndex
Index of the graph this point belongs to.
int GetN()
Returns GetN(), i.e. number of points of the total graph.
std::vector< size_t > fPointIndex
Index of the point within the graph this point corresponds to.
bool SetResidual(const bool &force=false)
TCalibrationGraphSet(TGraphErrors *graph=nullptr, const std::string &label="")
std::string fYAxisLabel
The label of the y-axis.
static EVerbosity fVerboseLevel
Changes verbosity.
std::vector< TCalibrationGraph > fResidualGraphs
These are the graphs used for plotting the residuals per source.
TGraphErrors * fTotalGraph
The sum of the other graphs, used for fitting.
double fMaximumX
Maximum x-value.
double fMinimumY
Minimum y-value.
TLine * fZeroResidual
! The line denoting zero residual
void Print(Option_t *opt="") const override
double fMaximumY
Maximum y-value.
std::string fXAxisLabel
The label of the x-axis.
void Scale(bool useAllPrevious=true)
int Add(TGraphErrors *, const std::string &label)
Add new graph to set, using the label when creating legends during plotting.
void DrawCalibration(Option_t *opt="", TLegend *legend=nullptr)
TGraphErrors * fTotalResidualGraph
The sum of the residuals. Not really used apart from plotting (but overlayed with the individual grap...
bool fResidualSet
Flag to indicate if the residual has been set correctly.
void Clear(Option_t *option="") override
double fMinimumX
Minimum x-value.
Int_t RemovePoint(const Int_t &px, const Int_t &py)
void SetAxisTitle(const char *title)
Set axis title for the graph (form "x-axis title;y-axis title")
std::vector< TCalibrationGraph > fGraphs
These are the graphs used for plotting the calibration points per source.
TF1 * FitFunction()
Gets the calibration from the total graph (might be nullptr!).
Int_t RemoveResidualPoint(const Int_t &px, const Int_t &py)
static EVerbosity VerboseLevel()
std::vector< std::string > fLabel
The labels for the different graphs.