NAGASH 0.9.8
Next Generation Analysis System
Loading...
Searching...
No Matches
NGHist.h
Go to the documentation of this file.
1//***************************************************************************************
4//***************************************************************************************
5
6#pragma once
7
8#include "NAGASH/HistBase.h"
9
10namespace NAGASH
11{
16 template <typename HistType>
17 class NGHist : public HistBase
18 {
19 friend class PlotGroup;
20
21 private:
22 HistType *FillHist = nullptr;
23 HistType *Nominal = nullptr;
24
25 std::map<TString, int> SystematicIndexMap;
26 std::vector<TString> SystematicVariationNames;
27 std::vector<HistType *> SystematicVariations;
28 std::vector<double> CorrelationFactors;
29 std::vector<std::vector<size_t>> SystematicVariationIndices;
30 std::vector<SystPolicy> SystematicVariationPolicies;
31 std::map<TString, int> SystematicVariationBackup;
32
33 HistType *Final = nullptr;
34 TH2D *Covariance = nullptr;
35 TH2D *Correlation = nullptr;
36 bool Processed = false;
37
38 // hist link
39 std::vector<std::shared_ptr<HistBase>> LinkedPlots;
40 TString LinkType = "";
41 std::function<double(const std::vector<double> &)> LinkContentFunction;
42 std::function<double(const std::vector<double> &, const std::vector<double> &)> LinkErrorFunction;
43 std::function<void(const std::vector<TH1 *> &, TH1 *)> LinkHistFunction;
44
45 protected:
46 std::shared_ptr<HistBase> CloneVirtual(const TString &name) override;
47
48 public:
50 template <typename... Args>
51 void Fill(Args &&...args) { FillHist->Fill(std::forward<Args>(args)...); }
52
54 HistType *GetNominalHist() { return Nominal; }
55
57 HistType *GetFillHist() { return FillHist; }
58
60 size_t GetNVariations() override { return SystematicVariationNames.size(); }
61
63 HistType *GetFinalHist()
64 {
65 if (!Final)
66 {
67 auto titleguard = MSGUser()->StartTitleWithGuard("NGHist::GetFinalHist");
68 MSGUser()->MSG_WARNING("Plot ", GetResultName(), " has not been processed, please call NGHist::Process. Returning nullptr");
69 }
70 return Final;
71 }
72
74 TH2D *GetCovariance()
75 {
76 if (!Covariance)
77 {
78 auto titleguard = MSGUser()->StartTitleWithGuard("NGHist::GetCovariance");
79 MSGUser()->MSG_WARNING("Plot ", GetResultName(), " has not been processed, please call NGHist::Process. Returning nullptr");
80 }
81 return Covariance;
82 }
83
85 TH2D *GetCorrelation()
86 {
87 if (!Correlation)
88 {
89 auto titleguard = MSGUser()->StartTitleWithGuard("NGHist::GetCorrelation");
90 MSGUser()->MSG_WARNING("Plot ", GetResultName(), " has not been processed, please call NGHist::Process. Returning nullptr");
91 }
92 return Correlation;
93 }
94
96 template <typename... Args>
97 void SetBinContent(Args &&...args) { FillHist->SetBinContent(std::forward<Args>(args)...); }
99 template <typename... Args>
100 void SetBinError(Args &&...args) { FillHist->SetBinError(std::forward<Args>(args)...); }
102 template <typename... Args>
103 double GetBinContent(Args &&...args) { return FillHist->GetBinContent(std::forward<Args>(args)...); }
105 template <typename... Args>
106 double GetBinError(Args &&...args) { return FillHist->GetBinError(std::forward<Args>(args)...); }
108 template <typename... Args>
109 double GetBinNominalError(Args &&...args) { return Nominal->GetBinError(std::forward<Args>(args)...); }
110
112 template <typename... Args>
113 double GetBinTotalError(Args &&...args)
114 {
115 if (Final != nullptr)
116 return Final->GetBinError(std::forward<Args>(args)...);
117 else
118 return Nominal->GetBinError(std::forward<Args>(args)...);
119 }
120
122 template <typename... Args>
123 double GetBinFinalError(Args &&...args)
124 {
125 if (Final != nullptr)
126 return Final->GetBinError(std::forward<Args>(args)...);
127 else
128 return Nominal->GetBinError(std::forward<Args>(args)...);
129 }
130
132 template <typename... Args>
133 double GetBinSystematicError(Args &&...args)
134 {
135 if (Final != nullptr)
136 {
137 double finalerror = Final->GetBinError(std::forward<Args>(args)...);
138 double staterror = Nominal->GetBinError(std::forward<Args>(args)...);
139 return std::sqrt(finalerror * finalerror - staterror * staterror);
140 }
141 else
142 return 0;
143 }
144
146 TString GetVariationName(uint64_t index) override
147 {
148 if (index == 0)
149 return "Nominal";
150 if (index > SystematicVariationNames.size() || index < 0)
151 return "UNKNOWN";
152 else
153 return SystematicVariationNames[index - 1];
154 }
155
157 std::vector<TString> GetVariationNames()
158 {
159 return SystematicVariationNames;
160 }
161
164 HistType *GetVariation(const TString &name)
165 {
166 if (name == "Nominal")
167 return Nominal;
168 if (auto findresult = SystematicIndexMap.find(name); findresult != SystematicIndexMap.end())
169 {
170 return SystematicVariations[findresult->second];
171 }
172 else
173 {
174 auto guard = MSGUser()->StartTitleWithGuard("NGHist::GetVariation");
175 MSGUser()->MSG_WARNING("NGHist ", GetResultName(), " does not have variation ", name, ", returning nullptr");
176 return nullptr;
177 }
178 }
179
182 HistType *GetVariation(uint64_t index)
183 {
184 if (index == 0)
185 return Nominal;
186 else if (index >= 1 && index <= SystematicVariations.size())
187 return SystematicVariations[index - 1];
188 else
189 {
190 auto guard = MSGUser()->StartTitleWithGuard("NGHist::GetVariation");
191 MSGUser()->MSG_WARNING("NGHist ", GetResultName(), " does not have variation ", index, ", returning nullptr");
192 return nullptr;
193 }
194 }
195
198 double GetVariationCorrelationFactor(const TString &name) override
199 {
200 if (auto findresult = SystematicIndexMap.find(name); findresult != SystematicIndexMap.end())
201 {
202 return CorrelationFactors[findresult->second];
203 }
204 else
205 {
206 auto guard = MSGUser()->StartTitleWithGuard("NGHist::GetVariationCorrelationFactor");
207 MSGUser()->MSG_ERROR("NGHist ", GetResultName(), " does not have variation ", name, ", returning zero");
208 return 0;
209 }
210 }
211
213 virtual ~NGHist()
214 {
215 if (Nominal)
216 delete Nominal;
217
218 for (auto &h : SystematicVariations)
219 if (h)
220 delete h;
221
222 if (Processed)
223 {
224 if (Final)
225 delete Final;
226 if (Covariance)
227 delete Covariance;
228 if (Correlation)
229 delete Correlation;
230 }
231 }
232
233 // default constructor
234 template <typename... Args>
235 NGHist(std::shared_ptr<MSGTool> MSG, std::shared_ptr<ConfigTool> c, const TString &name, const TString &filename, const Args &...args);
236 // auxiliary constructor
237 NGHist(std::shared_ptr<MSGTool> MSG, std::shared_ptr<ConfigTool> c, const TString &name);
238
239 // helper function
240 std::vector<std::vector<double>> GetPartialCovariance(const std::vector<TString> &systlist);
241 std::shared_ptr<NGHist<HistType>> Rebin(const TString &name, int rebinnum);
242 std::shared_ptr<NGHist<HistType>> Rebin2D(const TString &name, int rebinnum_x, int rebinum_y);
243 std::shared_ptr<NGHist<HistType>> Rebin3D(const TString &name, int rebinnum_x, int rebinum_y, int rebinum_z);
244 std::shared_ptr<NGHist<HistType>> Rebin(const TString &name, std::vector<double> bindiv);
245 std::shared_ptr<NGHist<HistType>> Rebin2D(const TString &name, std::vector<double> bindiv_x, std::vector<double> bindiv_y);
246 std::shared_ptr<NGHist<HistType>> Rebin3D(const TString &name, std::vector<double> bindiv_x, std::vector<double> bindiv_y, std::vector<double> bindiv_z);
247
248 // function inherited from HistBase
249 bool BookSystematicVariation(const TString &name, double corr_factor = 0, SystPolicy policy = SystPolicy::Maximum) override;
250 bool BookSystematicVariation(const TString &name1, const TString &name2, double corr_factor = 0, SystPolicy policy = SystPolicy::Maximum) override;
251 bool BookSystematicVariation(const std::vector<TString> &names, double corr_factor = 0, SystPolicy policy = SystPolicy::Maximum) override;
252 void ClearLinkPlot() override;
253 std::shared_ptr<NGHist<HistType>> Clone(const TString &name);
254 void Combine(std::shared_ptr<Result> result) override;
255 TH1 *GetNominalHistVirtual() override;
256 std::vector<std::tuple<std::vector<TString>, double, SystPolicy>> GetVariationTypes() override;
257 TH1 *GetVariationVirtual(const TString &name) override;
258 TH1 *GetVariationVirtual(uint64_t index) override;
259 bool IsBackupVariation(uint64_t) override;
260 bool IsBackupVariation(const TString &) override;
261 bool Process() override;
262 bool Recover() override;
263 bool Recover(TFile *file) override;
264 bool Recover(const TString &filename) override;
265 bool Recover(std::shared_ptr<TFileHelper>) override;
266 bool RegroupSystematicVariation(const std::vector<TString> &names, SystPolicy policy) override;
267 bool RemoveSystematicVariation(const TString &name) override;
268 bool RenameSystematicVariation(const TString &name_old, const TString &name_new) override;
269 void Reset() override;
270 void Scale(double) override;
271 void SetLinkPlot(uint64_t index, std::shared_ptr<HistBase> subplot) override;
272 void SetLinkType(const TString &type) override;
273 void SetLinkType(std::function<double(const std::vector<double> &)>, std::function<double(const std::vector<double> &, const std::vector<double> &)>) override;
274 void SetLinkType(std::function<void(const std::vector<TH1 *> &, TH1 *)>) override;
275 bool SetSystematicVariation(const TString &name = "Nominal") override;
276 bool SetSystematicVariation(uint64_t index) override;
277 void SetUnprocessed() override;
278 void Write() override;
279 void WriteToFile() override;
280
282 static double Variance(const std::vector<double> &v)
283 {
284 if (v.size() == 0)
285 return 0;
286
287 double mean = 0;
288 double mean_square = 0;
289
290 for (auto &value : v)
291 {
292 mean += value;
293 mean_square += value * value;
294 }
295
296 mean /= v.size();
297 mean_square /= v.size();
298
299 return (mean_square - mean * mean);
300 }
301 };
302
303 // protected function
304 template <typename HistType>
305 inline std::shared_ptr<HistBase> NGHist<HistType>::CloneVirtual(const TString &name)
306 {
307 std::shared_ptr<NGHist<HistType>> newPlot = std::shared_ptr<NGHist<HistType>>(new NGHist<HistType>(MSGUser(), ConfigUser(), name));
308 newPlot->SetOutputFileName(GetOutputFileName());
309 delete newPlot->Nominal;
310 newPlot->Nominal = nullptr;
311 newPlot->Nominal = (HistType *)this->Nominal->Clone(name);
312 newPlot->FillHist = newPlot->Nominal;
313
314 // copy systs
315 newPlot->SystematicIndexMap = this->SystematicIndexMap;
316 newPlot->SystematicVariationNames = this->SystematicVariationNames;
317 newPlot->CorrelationFactors = this->CorrelationFactors;
318 newPlot->SystematicVariationIndices = this->SystematicVariationIndices;
319 newPlot->SystematicVariationPolicies = this->SystematicVariationPolicies;
320 newPlot->SystematicVariationBackup = this->SystematicVariationBackup;
321
322 for (int i = 0; i < this->SystematicVariationNames.size(); i++)
323 newPlot->SystematicVariations.emplace_back((HistType *)(this->SystematicVariations[i]->Clone(name + "_Syst_" + this->SystematicVariationNames[i])));
324
325 // copy links
326 newPlot->LinkedPlots = this->LinkedPlots;
327 newPlot->LinkType = this->LinkType;
328 newPlot->LinkContentFunction = this->LinkContentFunction;
329 newPlot->LinkErrorFunction = this->LinkErrorFunction;
330 newPlot->LinkHistFunction = this->LinkHistFunction;
331
332 // copy processed plots
333 if (Processed)
334 {
335 newPlot->Final = (HistType *)this->Final->Clone("Final_" + name);
336 if (this->Covariance)
337 newPlot->Covariance = (TH2D *)this->Covariance->Clone("Covariance_" + name);
338 if (this->Correlation)
339 newPlot->Correlation = (TH2D *)this->Correlation->Clone("Correlation_" + name);
340 newPlot->Processed = true;
341 }
342
343 return newPlot;
344 }
345
353 template <typename HistType>
354 template <typename... Args>
355 inline NGHist<HistType>::NGHist(std::shared_ptr<MSGTool> MSG, std::shared_ptr<ConfigTool> c, const TString &name, const TString &filename, const Args &...args)
356 : HistBase(MSG, c, name, filename)
357 {
358 static_assert(std::is_base_of<TH1, HistType>::value, "NGHist : the template type must be inherited from ROOT::TH1");
359 constexpr size_t argcount = sizeof...(Args);
360 if constexpr (argcount != 0)
361 Nominal = new HistType(GetResultName(), GetResultName(), args...);
362 else // if no more arguments provided, use the default constructor of TH1
363 Nominal = new HistType();
364 Nominal->SetDirectory(0);
365 Nominal->Sumw2();
366 FillHist = Nominal;
367 }
368
373 template <typename HistType>
374 inline NGHist<HistType>::NGHist(std::shared_ptr<MSGTool> MSG, std::shared_ptr<ConfigTool> c, const TString &name) : HistBase(MSG, c, name, "")
375 {
376 static_assert(std::is_base_of<TH1, HistType>::value, "NGHist : the template type must be inherited from ROOT::TH1");
377 Nominal = new HistType();
378 Nominal->SetDirectory(0);
379 FillHist = Nominal;
380 }
381
385 template <typename HistType>
386 inline std::vector<std::vector<double>> NGHist<HistType>::GetPartialCovariance(const std::vector<TString> &systlist)
387 {
388 auto titleguard = MSGUser()->StartTitleWithGuard("NGHist::GetPartialCovariance");
389
390 // find the indices
391 bool is_contain_nominal = false;
392 std::vector<size_t> inputsystindices;
393 for (size_t i = 0; i < systlist.size(); i++)
394 {
395 auto findresult = SystematicIndexMap.find(systlist[i]);
396 if (findresult != SystematicIndexMap.end())
397 inputsystindices.emplace_back(findresult->second);
398 else if (systlist[i] != "Nominal")
399 MSGUser()->MSG_WARNING("Syst ", systlist[i], " does not exist in plot ", GetResultName(), ", will skip that.");
400 else
401 is_contain_nominal = true;
402
403 // check for backup
404 findresult = SystematicVariationBackup.find(systlist[i]);
405 if (findresult != SystematicVariationBackup.end())
406 MSGUser()->MSG_WARNING("Syst ", systlist[i], " is in the backup list of plot ", GetResultName(), ", will skip that.");
407 }
408
409 // sort and remove duplicates
410 std::sort(inputsystindices.begin(), inputsystindices.end());
411 auto last = std::unique(inputsystindices.begin(), inputsystindices.end());
412 inputsystindices.erase(last, inputsystindices.end());
413
414 // check and get the correct syst pair
415 std::vector<std::vector<size_t>> input_syst_indices;
416 std::vector<SystPolicy> input_syst_policies;
417 int systcount = 0;
418 for (auto &indices : this->SystematicVariationIndices)
419 {
420 for (auto &index : indices)
421 if (std::binary_search(inputsystindices.begin(), inputsystindices.end(), index))
422 {
423 input_syst_indices.emplace_back(indices);
424 input_syst_policies.emplace_back(SystematicVariationPolicies[systcount]);
425 break;
426 }
427 systcount++;
428 }
429
430 // calculate cov
431 int Nbins = Nominal->GetNbinsX() * Nominal->GetNbinsY() * Nominal->GetNbinsY();
432 auto histtool = ToolkitUser()->template GetTool<HistTool>();
433 auto convert2vector = [=](TH1 *h)
434 {
435 std::vector<double> content;
436 content.reserve(Nbins);
437 for (int i = 1; i <= h->GetXaxis()->GetNbins(); i++)
438 for (int j = 1; j <= h->GetYaxis()->GetNbins(); j++)
439 for (int k = 1; k <= h->GetZaxis()->GetNbins(); k++)
440 content.emplace_back(h->GetBinContent(i, j, k));
441 return content;
442 };
443
444 std::vector<std::vector<double>> cov(Nbins, std::vector<double>(Nbins, 0));
445 std::vector<double> nominal_content = convert2vector(Nominal);
446
447 systcount = 0;
448 for (auto &indices : input_syst_indices)
449 {
450 double corr = CorrelationFactors[indices[0]];
451 SystPolicy policy = input_syst_policies[systcount];
452 std::vector<std::vector<double>> syst_contents;
453 for (auto &index : indices)
454 syst_contents.emplace_back(convert2vector(SystematicVariations[index]));
455
456 std::vector<double> vdeltai(indices.size());
457 std::vector<double> vdeltaj(indices.size());
458 std::vector<double> vsystcontenti(indices.size());
459 std::vector<double> vsystcontentj(indices.size());
460
461 for (int i = 0; i < Nbins; i++)
462 {
463 for (int j = 0; j < Nbins; j++)
464 {
465 double deltai, deltaj;
466 int signi = syst_contents[0][i] - nominal_content[i] > 0 ? 1 : -1;
467 int signj = syst_contents[0][j] - nominal_content[j] > 0 ? 1 : -1;
468 for (int k = 0; k < indices.size(); k++)
469 {
470 vdeltai[k] = fabs(syst_contents[k][i] - nominal_content[i]);
471 vdeltaj[k] = fabs(syst_contents[k][j] - nominal_content[j]);
472 vsystcontenti[k] = syst_contents[k][i];
473 vsystcontentj[k] = syst_contents[k][j];
474 }
475
476 if (policy == SystPolicy::Maximum)
477 {
478 deltai = *std::max_element(vdeltai.begin(), vdeltai.end());
479 deltaj = *std::max_element(vdeltaj.begin(), vdeltaj.end());
480 }
481 else if (policy == SystPolicy::Average)
482 {
483 deltai = std::accumulate(vdeltai.begin(), vdeltai.end(), (double)0) / indices.size();
484 deltaj = std::accumulate(vdeltaj.begin(), vdeltaj.end(), (double)0) / indices.size();
485 }
486 else if (policy == SystPolicy::Variance)
487 {
488 deltai = std::sqrt(std::fabs(Variance(vsystcontenti)));
489 deltaj = std::sqrt(std::fabs(Variance(vsystcontentj)));
490 }
491 else if (policy == SystPolicy::CTEQ)
492 {
493 deltai = (vsystcontenti[0] - vsystcontenti[1]) / 2;
494 deltaj = (vsystcontentj[0] - vsystcontentj[1]) / 2;
495 signi = 1;
496 signj = 1;
497 }
498
499 if (i != j)
500 cov[i][j] += signi * signj * deltai * deltaj * corr;
501 else
502 cov[i][j] += signi * signj * deltai * deltaj;
503 }
504 }
505
506 systcount++;
507 }
508
509 // add nominal stat unc. if it is given
510 if (is_contain_nominal)
511 {
512 std::vector<double> nominal_error;
513 nominal_error.reserve(Nbins);
514 for (int i = 1; i <= Nominal->GetXaxis()->GetNbins(); i++)
515 for (int j = 1; j <= Nominal->GetYaxis()->GetNbins(); j++)
516 for (int k = 1; k <= Nominal->GetZaxis()->GetNbins(); k++)
517 nominal_error.emplace_back(Nominal->GetBinError(i, j, k));
518
519 for (int i = 0; i < Nbins; i++)
520 cov[i][i] += nominal_error[i] * nominal_error[i];
521 }
522
523 return cov;
524 }
525
530 template <typename HistType>
531 inline std::shared_ptr<NGHist<HistType>> NGHist<HistType>::Rebin(const TString &name, int rebinnum)
532 {
533 static_assert(!std::is_base_of<TH2, HistType>::value, "NGHist : function Rebin is not supported for 2D histograms");
534 static_assert(!std::is_base_of<TH3, HistType>::value, "NGHist : function Rebin is not supported for 3D histograms");
535
536 auto titleguard = MSGUser()->StartTitleWithGuard("NGHist::Rebin");
537
538 if (rebinnum <= 0)
539 {
540 MSGUser()->MSG_ERROR("Rebin number is not legal(", rebinnum, "), will do nothing for plot ", GetResultName());
541 return nullptr;
542 }
543
544 if (!Processed)
545 Process();
546
547 auto newh = this->Clone(name);
548 newh->ClearLinkPlot();
549 newh->SetUnprocessed();
550
551 newh->Nominal->Rebin(rebinnum);
552 for (auto &h : newh->SystematicVariations)
553 h->Rebin(rebinnum);
554
555 newh->Process();
556
557 return newh;
558 }
559
565 template <typename HistType>
566 inline std::shared_ptr<NGHist<HistType>> NGHist<HistType>::Rebin2D(const TString &name, int rebinnum_x, int rebinnum_y)
567 {
568 static_assert(!std::is_base_of<TH1, HistType>::value || std::is_base_of<TH2, HistType>::value, "NGHist : function Rebin2D is not supported for 1D histograms");
569 static_assert(!std::is_base_of<TH3, HistType>::value, "NGHist : function Rebin2D is not supported for 3D histograms");
570
571 auto titleguard = MSGUser()->StartTitleWithGuard("NGHist::Rebin2D");
572
573 if (rebinnum_x <= 0 || rebinnum_y <= 0)
574 {
575 MSGUser()->MSG_ERROR("Rebin number is not legal, will do nothing for plot ", GetResultName());
576 return nullptr;
577 }
578
579 if (!Processed)
580 Process();
581
582 auto newh = this->Clone(name);
583 newh->ClearLinkPlot();
584 newh->SetUnprocessed();
585
586 newh->Nominal->Rebin(rebinnum_x, rebinnum_y);
587 for (auto &h : newh->SystematicVariations)
588 h->Rebin(rebinnum_x, rebinnum_y);
589
590 newh->Process();
591
592 return newh;
593 }
594
601 template <typename HistType>
602 inline std::shared_ptr<NGHist<HistType>> NGHist<HistType>::Rebin3D(const TString &name, int rebinnum_x, int rebinnum_y, int rebinnum_z)
603 {
604 static_assert(!std::is_base_of<TH1, HistType>::value || std::is_base_of<TH3, HistType>::value, "NGHist : function Rebin3D is not supported for 1D histograms");
605 static_assert(!std::is_base_of<TH2, HistType>::value || std::is_base_of<TH3, HistType>::value, "NGHist : function Rebin3D is not supported for 2D histograms");
606
607 auto titleguard = MSGUser()->StartTitleWithGuard("NGHist::Rebin3D");
608
609 if (rebinnum_x <= 0 || rebinnum_y <= 0 || rebinnum_z <= 0)
610 {
611 MSGUser()->MSG_ERROR("Rebin number is not legal, will do nothing for plot ", GetResultName());
612 return nullptr;
613 }
614
615 if (!Processed)
616 Process();
617
618 auto newh = this->Clone(name);
619 newh->ClearLinkPlot();
620 newh->SetUnprocessed();
621
622 newh->Nominal->Rebin(rebinnum_x, rebinnum_y, rebinnum_z);
623 for (auto &h : newh->SystematicVariations)
624 h->Rebin(rebinnum_x, rebinnum_y, rebinnum_z);
625
626 newh->Process();
627
628 return newh;
629 }
630
635 template <typename HistType>
636 inline std::shared_ptr<NGHist<HistType>> NGHist<HistType>::Rebin(const TString &name, std::vector<double> bindiv)
637 {
638 static_assert(!std::is_base_of<TH2, HistType>::value, "NGHist : function Rebin is not supported for 2D histograms");
639 static_assert(!std::is_base_of<TH3, HistType>::value, "NGHist : function Rebin is not supported for 3D histograms");
640
641 auto titleguard = MSGUser()->StartTitleWithGuard("NGHist::Rebin");
642
643 if (bindiv.size() == 0)
644 {
645 MSGUser()->MSG_ERROR("No bin division is given, will do nothing for plot ", GetResultName());
646 return nullptr;
647 }
648
649 if (!Processed)
650 Process();
651
652 std::shared_ptr<NGHist<HistType>> newPlot = std::shared_ptr<NGHist<HistType>>(new NGHist<HistType>(MSGUser(), ConfigUser(), name));
653 newPlot->SetOutputFileName(GetOutputFileName());
654 delete newPlot->Nominal;
655 newPlot->Nominal = nullptr;
656 newPlot->Nominal = (HistType *)this->Nominal->Rebin(bindiv.size() - 1, name, bindiv.data());
657 newPlot->FillHist = newPlot->Nominal;
658
659 // copy systs
660 newPlot->SystematicIndexMap = this->SystematicIndexMap;
661 newPlot->SystematicVariationNames = this->SystematicVariationNames;
662 newPlot->CorrelationFactors = this->CorrelationFactors;
663 newPlot->SystematicVariationIndices = this->SystematicVariationIndices;
664 newPlot->SystematicVariationPolicies = this->SystematicVariationPolicies;
665 newPlot->SystematicVariationBackup = this->SystematicVariationBackup;
666
667 for (int i = 0; i < this->SystematicVariationNames.size(); i++)
668 newPlot->SystematicVariations.emplace_back((HistType *)(this->SystematicVariations[i]->Rebin(bindiv.size() - 1, name + "_Syst_" + this->SystematicVariationNames[i], bindiv.data())));
669
670 newPlot->Process();
671
672 return newPlot;
673 }
674
680 template <typename HistType>
681 inline std::shared_ptr<NGHist<HistType>> NGHist<HistType>::Rebin2D(const TString &name, std::vector<double> bindiv_x, std::vector<double> bindiv_y)
682 {
683 static_assert(!std::is_base_of<TH1, HistType>::value || std::is_base_of<TH2, HistType>::value, "NGHist : function Rebin2D is not supported for 1D histograms");
684 static_assert(!std::is_base_of<TH3, HistType>::value, "NGHist : function Rebin2D is not supported for 3D histograms");
685
686 auto titleguard = MSGUser()->StartTitleWithGuard("NGHist::Rebin2D");
687
688 if (bindiv_x.size() == 0 && bindiv_y.size() == 0)
689 {
690 MSGUser()->MSG_ERROR("No bin division is given, will do nothing for plot ", GetResultName());
691 return nullptr;
692 }
693 int nbinX, nbinY;
694 bool fixX = false, fixY = false;
695 double minX, maxX, minY, maxY;
696 if (bindiv_x.size() == 0)
697 {
698 auto bindivarr = Nominal->GetXaxis()->GetXbins()->GetArray();
699 nbinX = Nominal->GetXaxis()->GetNbins();
700 if (bindivarr)
701 bindiv_x = std::vector<double>(bindivarr, bindivarr + nbinX + 1);
702 else
703 {
704 minX = Nominal->GetXaxis()->GetBinLowEdge(1);
705 maxX = Nominal->GetXaxis()->GetBinUpEdge(nbinX);
706 fixX = true;
707 }
708 }
709 if (bindiv_y.size() == 0)
710 {
711 auto bindivarr = Nominal->GetYaxis()->GetXbins()->GetArray();
712 nbinY = Nominal->GetYaxis()->GetNbins();
713 if (bindivarr)
714 bindiv_y = std::vector<double>(bindivarr, bindivarr + nbinY + 1);
715 else
716 {
717 minY = Nominal->GetYaxis()->GetBinLowEdge(1);
718 maxY = Nominal->GetYaxis()->GetBinUpEdge(nbinY);
719 fixY = true;
720 }
721 }
722
723 // check bin edge
724 if (!fixX)
725 for (int newX = 0, oldX = 1; newX < bindiv_x.size(); newX++)
726 {
727 if (bindiv_x.at(newX) < Nominal->GetXaxis()->GetBinLowEdge(1))
728 continue;
729 if (bindiv_x.at(newX) > Nominal->GetXaxis()->GetBinUpEdge(Nominal->GetXaxis()->GetNbins()))
730 break;
731 if (!TMath::AreEqualAbs(Nominal->GetXaxis()->GetBinLowEdge(oldX), bindiv_x.at(newX), TMath::Max(1.E-8 * Nominal->GetXaxis()->GetBinWidth(oldX), 1.E-16)) && Nominal->GetXaxis()->GetBinLowEdge(oldX) > bindiv_x.at(newX))
732 MSGUser()->MSG_WARNING("Bin X edge ", newX + 1, " of rebinned histogram does not match any bin edges of the old histogram. Result can be inconsistent");
733 else
734 oldX++;
735 }
736 if (!fixY)
737 for (int newY = 0, oldY = 1; newY < bindiv_y.size(); newY++)
738 {
739 if (bindiv_y.at(newY) < Nominal->GetYaxis()->GetBinLowEdge(1))
740 continue;
741 if (bindiv_y.at(newY) > Nominal->GetYaxis()->GetBinUpEdge(Nominal->GetYaxis()->GetNbins()))
742 break;
743 if (!TMath::AreEqualAbs(Nominal->GetYaxis()->GetBinLowEdge(oldY), bindiv_y.at(newY), TMath::Max(1.E-8 * Nominal->GetYaxis()->GetBinWidth(oldY), 1.E-16)) && Nominal->GetYaxis()->GetBinLowEdge(oldY) > bindiv_y.at(newY))
744 MSGUser()->MSG_WARNING("Bin Y edge ", newY + 1, " of rebinned histogram does not match any bin edges of the old histogram. Result can be inconsistent");
745 else
746 oldY++;
747 }
748
749 if (!Processed)
750 Process();
751
752 // rebin algorithm
753 auto Rebin = [](HistType *in, HistType *out)
754 {
755 auto xAxis = in->GetXaxis();
756 auto yAxis = in->GetYaxis();
757 for (int i = 1; i <= xAxis->GetNbins(); i++)
758 for (int j = 1; j <= yAxis->GetNbins(); j++)
759 {
760 double binCenterX = xAxis->GetBinCenter(i);
761 double binCenterY = yAxis->GetBinCenter(j);
762 int binX = out->GetXaxis()->FindBin(binCenterX);
763 int binY = out->GetYaxis()->FindBin(binCenterY);
764 out->SetBinContent(binX, binY, out->GetBinContent(binX, binY) + in->GetBinContent(i, j));
765 out->SetBinError(binX, binY, Uncertainty::AplusB(out->GetBinError(binX, binY), in->GetBinError(i, j)));
766 }
767 };
768
769 std::shared_ptr<NGHist<HistType>> newPlot;
770 if (fixX)
771 newPlot = std::make_shared<NGHist<HistType>>(MSGUser(), ConfigUser(), name, GetOutputFileName(), nbinX, minX, maxX, bindiv_y.size() - 1, bindiv_y.data());
772 else if (fixY)
773 newPlot = std::make_shared<NGHist<HistType>>(MSGUser(), ConfigUser(), name, GetOutputFileName(), bindiv_x.size() - 1, bindiv_x.data(), nbinY, minY, maxY);
774 else
775 newPlot = std::make_shared<NGHist<HistType>>(MSGUser(), ConfigUser(), name, GetOutputFileName(), bindiv_x.size() - 1, bindiv_x.data(), bindiv_y.size() - 1, bindiv_y.data());
776
777 Rebin(this->Nominal, newPlot->Nominal);
778 newPlot->FillHist = newPlot->Nominal;
779
780 // copy systs
781 newPlot->SystematicIndexMap = this->SystematicIndexMap;
782 newPlot->SystematicVariationNames = this->SystematicVariationNames;
783 newPlot->CorrelationFactors = this->CorrelationFactors;
784 newPlot->SystematicVariationIndices = this->SystematicVariationIndices;
785 newPlot->SystematicVariationPolicies = this->SystematicVariationPolicies;
786 newPlot->SystematicVariationBackup = this->SystematicVariationBackup;
787
788 newPlot->SystematicVariations.resize(this->SystematicVariations.size());
789 for (int i = 0; i < this->SystematicVariationNames.size(); i++)
790 Rebin(this->SystematicVariations.at(i), newPlot->SystematicVariations.at(i));
791
792 newPlot->Process();
793
794 return newPlot;
795 }
796
803 template <typename HistType>
804 inline std::shared_ptr<NGHist<HistType>> NGHist<HistType>::Rebin3D(const TString &name, std::vector<double> bindiv_x, std::vector<double> bindiv_y, std::vector<double> bindiv_z)
805 {
806 static_assert(!std::is_base_of<TH1, HistType>::value || std::is_base_of<TH3, HistType>::value, "NGHist : function Rebin3D is not supported for 1D histograms");
807 static_assert(!std::is_base_of<TH2, HistType>::value || std::is_base_of<TH3, HistType>::value, "NGHist : function Rebin3D is not supported for 2D histograms");
808
809 auto titleguard = MSGUser()->StartTitleWithGuard("NGHist::Rebin3D");
810
811 if (bindiv_x.size() == 0 && bindiv_y.size() == 0 && bindiv_z.size() == 0)
812 {
813 MSGUser()->MSG_ERROR("No bin division is given, will do nothing for plot ", GetResultName());
814 return nullptr;
815 }
816 if (bindiv_x.size() == 0)
817 {
818 auto bindivarr = Nominal->GetXaxis()->GetXbins()->GetArray();
819 int size = Nominal->GetXaxis()->GetNbins() + 1;
820 bindiv_x = std::vector<double>(bindivarr, bindivarr + size);
821 }
822 if (bindiv_y.size() == 0)
823 {
824 auto bindivarr = Nominal->GetYaxis()->GetXbins()->GetArray();
825 int size = Nominal->GetYaxis()->GetNbins() + 1;
826 bindiv_y = std::vector<double>(bindivarr, bindivarr + size);
827 }
828 if (bindiv_z.size() == 0)
829 {
830 auto bindivarr = Nominal->GetZaxis()->GetXbins()->GetArray();
831 int size = Nominal->GetZaxis()->GetNbins() + 1;
832 bindiv_z = std::vector<double>(bindivarr, bindivarr + size);
833 }
834
835 int oldX = 1, oldY = 1, oldZ = 1;
836 for (int newX = 0; newX < bindiv_x.size(); newX++)
837 {
838 if (bindiv_x.at(newX) < Nominal->GetXaxis()->GetBinLowEdge(1))
839 continue;
840 if (bindiv_x.at(newX) > Nominal->GetXaxis()->GetBinUpEdge(Nominal->GetXaxis()->GetNbins()))
841 break;
842 if (!TMath::AreEqualAbs(Nominal->GetXaxis()->GetBinLowEdge(oldX), bindiv_x.at(newX), TMath::Max(1.E-8 * Nominal->GetXaxis()->GetBinWidth(oldX), 1.E-16)) && Nominal->GetXaxis()->GetBinLowEdge(oldX) > bindiv_x.at(newX))
843 MSGUser()->MSG_WARNING("Bin X edge ", newX + 1, " of rebinned histogram does not match any bin edges of the old histogram. Result can be inconsistent");
844 else
845 oldX++;
846 }
847 for (int newY = 0; newY < bindiv_y.size(); newY++)
848 {
849 if (bindiv_y.at(newY) < Nominal->GetYaxis()->GetBinLowEdge(1))
850 continue;
851 if (bindiv_y.at(newY) > Nominal->GetYaxis()->GetBinUpEdge(Nominal->GetYaxis()->GetNbins()))
852 break;
853 if (!TMath::AreEqualAbs(Nominal->GetYaxis()->GetBinLowEdge(oldY), bindiv_y.at(newY), TMath::Max(1.E-8 * Nominal->GetYaxis()->GetBinWidth(oldY), 1.E-16)) && Nominal->GetYaxis()->GetBinLowEdge(oldY) > bindiv_y.at(newY))
854 MSGUser()->MSG_WARNING("Bin Y edge ", newY + 1, " of rebinned histogram does not match any bin edges of the old histogram. Result can be inconsistent");
855 else
856 oldY++;
857 }
858 for (int newZ = 0; newZ < bindiv_z.size(); newZ++)
859 {
860 if (bindiv_z.at(newZ) < Nominal->GetZaxis()->GetBinLowEdge(1))
861 continue;
862 if (bindiv_z.at(newZ) > Nominal->GetZaxis()->GetBinUpEdge(Nominal->GetZaxis()->GetNbins()))
863 break;
864 if (!TMath::AreEqualAbs(Nominal->GetZaxis()->GetBinLowEdge(oldZ), bindiv_z.at(newZ), TMath::Max(1.E-8 * Nominal->GetZaxis()->GetBinWidth(oldZ), 1.E-16)) && Nominal->GetZaxis()->GetBinLowEdge(oldZ) > bindiv_z.at(newZ))
865 MSGUser()->MSG_WARNING("Bin Z edge ", newZ + 1, " of rebinned histogram does not match any bin edges of the old histogram. Result can be inconsistent");
866 else
867 oldZ++;
868 }
869
870 if (!Processed)
871 Process();
872
873 // rebin algorithm
874 auto Rebin = [](HistType *in, HistType *out)
875 {
876 auto xAxis = in->GetXaxis();
877 auto yAxis = in->GetYaxis();
878 auto zAxis = in->GetZaxis();
879 for (int i = 1; i <= xAxis->GetNbins(); i++)
880 for (int j = 1; j <= yAxis->GetNbins(); j++)
881 for (int k = 1; k <= zAxis->GetNbins(); k++)
882 {
883 double binCenterX = xAxis->GetBinCenter(i);
884 double binCenterY = yAxis->GetBinCenter(j);
885 double binCenterZ = zAxis->GetBinCenter(k);
886 int binX = out->GetXaxis()->FindBin(binCenterX);
887 int binY = out->GetYaxis()->FindBin(binCenterY);
888 int binZ = out->GetZaxis()->FindBin(binCenterZ);
889 out->SetBinContent(binX, binY, binZ, out->GetBinContent(binX, binY, binZ) + in->GetBinContent(i, j, k));
890 out->SetBinError(binX, binY, binZ, Uncertainty::AplusB(out->GetBinError(binX, binY, binZ), in->GetBinError(i, j, k)));
891 }
892 };
893
894 std::shared_ptr<NGHist<HistType>> newPlot = std::make_shared<NGHist<HistType>>(MSGUser(), ConfigUser(), name, GetOutputFileName(), bindiv_x.size() - 1, bindiv_x.data(), bindiv_y.size() - 1, bindiv_y.data(), bindiv_z.size() - 1, bindiv_z.data());
895 Rebin(this->Nominal, newPlot->Nominal);
896 newPlot->FillHist = newPlot->Nominal;
897
898 // copy systs
899 newPlot->SystematicIndexMap = this->SystematicIndexMap;
900 newPlot->SystematicVariationNames = this->SystematicVariationNames;
901 newPlot->CorrelationFactors = this->CorrelationFactors;
902 newPlot->SystematicVariationIndices = this->SystematicVariationIndices;
903 newPlot->SystematicVariationPolicies = this->SystematicVariationPolicies;
904 newPlot->SystematicVariationBackup = this->SystematicVariationBackup;
905
906 newPlot->SystematicVariations.resize(this->SystematicVariations.size());
907 for (int i = 0; i < this->SystematicVariationNames.size(); i++)
908 Rebin(this->SystematicVariations.at(i), newPlot->SystematicVariations.at(i));
909
910 newPlot->Process();
911
912 return newPlot;
913 }
914
915 // virtual functions
916
923 template <typename HistType>
924 inline bool NGHist<HistType>::BookSystematicVariation(const TString &name, double corr_factor, SystPolicy policy)
925 {
926 return BookSystematicVariation(std::vector<TString>{name}, corr_factor, policy);
927 }
928
936 template <typename HistType>
937 inline bool NGHist<HistType>::BookSystematicVariation(const TString &name1, const TString &name2, double corr_factor, SystPolicy policy)
938 {
939 return BookSystematicVariation(std::vector<TString>{name1, name2}, corr_factor, policy);
940 }
941
948 template <typename HistType>
949 inline bool NGHist<HistType>::BookSystematicVariation(const std::vector<TString> &names, double corr_factor, SystPolicy policy)
950 {
951 auto titleguard = MSGUser()->StartTitleWithGuard("NGHist::BookSystematicVariation");
952 if (names.size() == 0)
953 {
954 MSGUser()->MSG_ERROR("No systematic name is given, will do nothing for plot ", GetResultName());
955 return false;
956 }
957 // check not duplicate
958 for (auto name : names)
959 {
960 if (SystematicIndexMap.find(name) != SystematicIndexMap.end())
961 {
962 MSGUser()->MSG_ERROR("Systematic variation ", name, " has already existed in plot ", GetResultName());
963 return false;
964 }
965
966 if (name == "Nominal")
967 {
968 MSGUser()->MSG_ERROR("can not book variation named Nominal for plot ", GetResultName());
969 return false;
970 }
971 }
972
973 if (names.size() != 2 && policy == SystPolicy::CTEQ)
974 {
975 MSGUser()->MSG_ERROR("The CTEQ policy only applies to two systematic variations at the same time for plot ", GetResultName());
976 return false;
977 }
978
979 SetUnprocessed();
980
981 if (policy != SystPolicy::Backup)
982 {
983 std::vector<size_t> indices(names.size());
984 std::iota(indices.begin(), indices.end(), CorrelationFactors.size());
985 SystematicVariationIndices.emplace_back(indices);
986 SystematicVariationPolicies.emplace_back(policy);
987 }
988
989 for (auto name : names)
990 {
991 // for backup policy, store all of them into special map backup
992 if (policy == SystPolicy::Backup)
993 SystematicVariationBackup.emplace(std::pair<TString, int>(name, CorrelationFactors.size()));
994
995 SystematicIndexMap.emplace(std::pair<TString, int>(name, CorrelationFactors.size()));
996 SystematicVariationNames.emplace_back(name);
997
998 TString SystHistName = GetResultName() + "_Syst_" + name;
999 HistType *SystHist = (HistType *)Nominal->Clone(SystHistName.TString::Data());
1000 SystHist->SetTitle(SystHistName.TString::Data());
1001 // SystHist->Sumw2();
1002 SystHist->SetDirectory(0);
1003 SystHist->Reset("ICESM");
1004 SystematicVariations.emplace_back(SystHist);
1005 CorrelationFactors.emplace_back(corr_factor);
1006 }
1007
1008 return true;
1009 }
1010
1012 template <typename HistType>
1013 inline void NGHist<HistType>::ClearLinkPlot()
1014 {
1015 LinkedPlots.clear();
1016 LinkType = "";
1017 LinkContentFunction = std::function<double(const std::vector<double> &)>();
1018 LinkErrorFunction = std::function<double(const std::vector<double> &, const std::vector<double> &)>();
1019 LinkHistFunction = std::function<void(const std::vector<TH1 *> &, TH1 *)>();
1020 }
1021
1023 template <typename HistType>
1024 inline std::shared_ptr<NGHist<HistType>> NGHist<HistType>::Clone(const TString &name)
1025 {
1026 return std::static_pointer_cast<NGHist<HistType>>(this->CloneVirtual(name));
1027 }
1028
1030 template <typename HistType>
1031 inline void NGHist<HistType>::Combine(std::shared_ptr<Result> result)
1032 {
1033 auto titleguard = MSGUser()->StartTitleWithGuard("NGHist::Combine");
1034 auto subPlot = std::dynamic_pointer_cast<typename std::remove_pointer<decltype(this)>::type>(result);
1035 if (subPlot != nullptr)
1036 {
1037 Nominal->Add(subPlot->GetNominalHist());
1038
1039 bool systCompatible = true;
1040 for (auto &[varnames, corr_factor, policy] : subPlot->GetVariationTypes())
1041 if (SystematicIndexMap.find(varnames[0]) == SystematicIndexMap.end() && SystematicVariationBackup.find(varnames[0]) == SystematicVariationBackup.end())
1042 if (systCompatible)
1043 systCompatible = this->BookSystematicVariation(varnames, corr_factor, policy);
1044
1045 if (!systCompatible)
1046 {
1047 MSGUser()->MSG_ERROR("the systematic variations of plot ", this->GetResultName(), " are not compatible with input plot ", subPlot->GetResultName());
1048 return;
1049 }
1050
1051 for (size_t i = 1; i <= subPlot->GetNVariations(); i++)
1052 {
1053 TString varname = subPlot->GetVariationName(i);
1054 this->GetVariation(varname)->Add(subPlot->GetVariation(i));
1055 }
1056
1057 SetUnprocessed();
1058 }
1059 }
1060
1062 template <typename HistType>
1063 inline TH1 *NGHist<HistType>::GetVariationVirtual(const TString &name)
1064 {
1065 return (TH1 *)GetVariation(name);
1066 }
1067
1069 template <typename HistType>
1070 inline TH1 *NGHist<HistType>::GetVariationVirtual(uint64_t index)
1071 {
1072 return (TH1 *)GetVariation(index);
1073 }
1074
1076 template <typename HistType>
1077 inline TH1 *NGHist<HistType>::GetNominalHistVirtual()
1078 {
1079 return (TH1 *)Nominal;
1080 }
1081
1084 template <typename HistType>
1085 inline std::vector<std::tuple<std::vector<TString>, double, HistBase::SystPolicy>> NGHist<HistType>::GetVariationTypes()
1086 {
1087 std::vector<std::tuple<std::vector<TString>, double, SystPolicy>> types;
1088
1089 int systcount = 0;
1090 for (auto &indices : this->SystematicVariationIndices)
1091 {
1092 double corr = CorrelationFactors[indices[0]];
1093 SystPolicy policy = SystematicVariationPolicies[systcount];
1094 std::vector<TString> varnames;
1095
1096 for (auto &index : indices)
1097 varnames.emplace_back(SystematicVariationNames[index]);
1098
1099 types.emplace_back(std::tuple<std::vector<TString>, double, SystPolicy>(varnames, corr, policy));
1100
1101 ++systcount;
1102 }
1103
1104 // emplace backup
1105 for (auto &iter : SystematicVariationBackup)
1106 types.emplace_back(std::tuple<std::vector<TString>, double, SystPolicy>({iter.first}, CorrelationFactors[iter.second], HistBase::SystPolicy::Backup));
1107
1108 return types;
1109 }
1110
1112 template <typename HistType>
1113 inline bool NGHist<HistType>::IsBackupVariation(uint64_t index)
1114 {
1115 TString name;
1116 if (index == 0)
1117 return false;
1118 else if (index > 0 && index <= SystematicVariations.size())
1119 name = SystematicVariationNames[index - 1];
1120 else
1121 {
1122 auto guard = MSGUser()->StartTitleWithGuard("NGHist::IsBackupVariation");
1123 MSGUser()->MSG_WARNING("NGHist ", GetResultName(), " does not have variation ", index);
1124 FillHist = Nominal;
1125 return false;
1126 }
1127
1128 if (auto findresult = SystematicVariationBackup.find(name); findresult != SystematicVariationBackup.end())
1129 return true;
1130 else
1131 return false;
1132 }
1133
1135 template <typename HistType>
1136 inline bool NGHist<HistType>::IsBackupVariation(const TString &name)
1137 {
1138 if (auto findresult = SystematicVariationBackup.find(name); findresult != SystematicVariationBackup.end())
1139 return true;
1140 else
1141 return false;
1142 }
1143
1154 template <typename HistType>
1155 inline bool NGHist<HistType>::Process()
1156 {
1157 if (Processed)
1158 return true;
1159
1160 auto titleguard = MSGUser()->StartTitleWithGuard("NGHist::Process");
1161 auto histtool = ToolkitUser()->template GetTool<HistTool>();
1162 if (LinkedPlots.size() > 0)
1163 {
1164 for (size_t i = 0; i < LinkedPlots.size(); i++)
1165 LinkedPlots[i]->Process();
1166
1167 // check for variations
1168 bool systCompatible = true;
1169 for (size_t i = 0; i < LinkedPlots.size(); i++)
1170 for (auto &[varnames, corr_factor, policy] : LinkedPlots[i]->GetVariationTypes())
1171 if (SystematicIndexMap.find(varnames[0]) == SystematicIndexMap.end() && SystematicVariationBackup.find(varnames[0]) == SystematicVariationBackup.end())
1172 if (systCompatible)
1173 systCompatible = this->BookSystematicVariation(varnames, corr_factor, policy);
1174
1175 if (!systCompatible)
1176 {
1177 MSGUser()->MSG_ERROR("the systematic variations of plot ", GetResultName(), " are not compatible with those of the linked plots");
1178 return false;
1179 }
1180
1181 // Process as LinkedType with HistTool
1182 std::vector<TH1 *> LinkedHists;
1183 for (size_t i = 0; i < LinkedPlots.size(); i++)
1184 {
1185 if (LinkedPlots[i] == nullptr)
1186 LinkedHists.emplace_back(nullptr);
1187 else
1188 LinkedHists.emplace_back(LinkedPlots[i]->GetNominalHistVirtual());
1189 }
1190 if (LinkType != "")
1191 histtool->ProcessHistLink(LinkType, LinkedHists, GetNominalHist());
1192 else if (LinkContentFunction && LinkErrorFunction)
1193 histtool->ProcessHistLink(LinkContentFunction, LinkErrorFunction, LinkedHists, GetNominalHist());
1194 else if (LinkHistFunction)
1195 histtool->ProcessHistLink(LinkHistFunction, LinkedHists, GetNominalHist());
1196 else
1197 {
1198 MSGUser()->MSG_ERROR("Link type is not set for plot ", GetResultName(), " the hist link will not be processed");
1199 return false;
1200 }
1201 LinkedHists.clear();
1202
1203 for (size_t i = 1; i <= GetNVariations(); i++)
1204 {
1205 TString vname = GetVariationName(i);
1206 auto level = MSGUser()->Level();
1207 MSGUser()->SetOutputLevel(MSGLevel::SILENT);
1208 size_t linkSystCount = 0;
1209 for (size_t j = 0; j < LinkedPlots.size(); j++)
1210 {
1211 if (LinkedPlots[j] == nullptr)
1212 LinkedHists.emplace_back(nullptr);
1213 else
1214 {
1215 auto varh = LinkedPlots[j]->GetVariationVirtual(vname);
1216 if (!varh)
1217 LinkedHists.emplace_back(LinkedPlots[j]->GetNominalHistVirtual());
1218 else
1219 {
1220 LinkedHists.emplace_back(varh);
1221 ++linkSystCount;
1222 }
1223 }
1224 }
1225 MSGUser()->SetOutputLevel(level);
1226
1227 if (linkSystCount == 0)
1228 MSGUser()->MSG_WARNING("No variation ", vname, " is found in all linked plots for plot ", GetResultName(), ", will use nominal instead.");
1229
1230 if (LinkType != "")
1231 histtool->ProcessHistLink(LinkType, LinkedHists, GetVariation(vname));
1232 else if (LinkContentFunction && LinkErrorFunction)
1233 histtool->ProcessHistLink(LinkContentFunction, LinkErrorFunction, LinkedHists, GetVariation(vname));
1234 else if (LinkHistFunction)
1235 histtool->ProcessHistLink(LinkHistFunction, LinkedHists, GetVariation(vname));
1236
1237 LinkedHists.clear();
1238 }
1239 }
1240
1241 // Process final errors
1242 TString fname = "Final_" + GetResultName();
1243 Final = (HistType *)Nominal->Clone(fname.TString::Data());
1244 Final->SetDirectory(0);
1245 Final->SetTitle(fname.TString::Data());
1246
1247 int Nbins = Nominal->GetXaxis()->GetNbins() * Nominal->GetYaxis()->GetNbins() * Nominal->GetZaxis()->GetNbins();
1248
1249 // Process Total Uncertainty
1250 for (int i = 1; i <= Nominal->GetXaxis()->GetNbins(); i++)
1251 for (int j = 1; j <= Nominal->GetYaxis()->GetNbins(); j++)
1252 for (int k = 1; k <= Nominal->GetZaxis()->GetNbins(); k++)
1253 {
1254 double val = Nominal->GetBinContent(i, j, k);
1255 double unc = Nominal->GetBinError(i, j, k);
1256
1257 for (size_t ivar = 0; ivar < SystematicVariationIndices.size(); ivar++)
1258 {
1259 std::vector<double> deltas(SystematicVariationIndices[ivar].size());
1260 std::vector<double> systcontents(SystematicVariationIndices[ivar].size());
1261 // Add Syst. Unc.
1262 std::transform(SystematicVariationIndices[ivar].begin(), SystematicVariationIndices[ivar].end(), deltas.begin(), [=](size_t index)
1263 { return std::fabs(SystematicVariations[index]->GetBinContent(i, j, k) - val); });
1264 std::transform(SystematicVariationIndices[ivar].begin(), SystematicVariationIndices[ivar].end(), systcontents.begin(), [=](size_t index)
1265 { return SystematicVariations[index]->GetBinContent(i, j, k); });
1266
1267 if (SystematicVariationPolicies[ivar] == SystPolicy::Average)
1268 unc = Uncertainty::AplusB(unc, std::accumulate(deltas.begin(), deltas.end(), (double)0) / (double)(deltas.size()));
1269 if (SystematicVariationPolicies[ivar] == SystPolicy::Maximum)
1270 unc = Uncertainty::AplusB(unc, *std::max_element(deltas.begin(), deltas.end()));
1271 if (SystematicVariationPolicies[ivar] == SystPolicy::Variance)
1272 unc = Uncertainty::AplusB(unc, std::sqrt(std::fabs(Variance(systcontents))));
1273 if (SystematicVariationPolicies[ivar] == SystPolicy::CTEQ)
1274 unc = Uncertainty::AplusB(unc, (systcontents[0] - systcontents[1]) / 2);
1275 }
1276
1277 Final->SetBinContent(i, j, k, val);
1278 Final->SetBinError(i, j, k, unc);
1279 }
1280
1281 // Process correlation and covariance if syst booked and not too many bins
1282 // since it consumes so much time when bin number is great
1283 if (SystematicVariations.size() > 0 && Nbins <= 1e3)
1284 {
1285 fname = "Correlation_" + GetResultName();
1286 Correlation = new TH2D(fname.TString::Data(), fname.TString::Data(),
1287 Nbins, 0, Nbins,
1288 Nbins, 0, Nbins);
1289 Correlation->SetDirectory(0);
1290 Correlation->Sumw2();
1291
1292 fname = "Covariance_" + GetResultName();
1293 Covariance = new TH2D(fname.TString::Data(), fname.TString::Data(),
1294 Nbins, 0, Nbins,
1295 Nbins, 0, Nbins);
1296 Covariance->SetDirectory(0);
1297 Covariance->Sumw2();
1298
1299 histtool->ConvertCovToCorr(Covariance, Correlation);
1300
1301 std::vector<std::vector<double>> cov(Nbins, std::vector<double>(Nbins, 0));
1302
1303 auto convert2vector = [=](TH1 *h)
1304 {
1305 std::vector<double> content;
1306 content.reserve(Nbins);
1307 for (int i = 1; i <= h->GetXaxis()->GetNbins(); i++)
1308 for (int j = 1; j <= h->GetYaxis()->GetNbins(); j++)
1309 for (int k = 1; k <= h->GetZaxis()->GetNbins(); k++)
1310 content.emplace_back(h->GetBinContent(i, j, k));
1311 return content;
1312 };
1313
1314 auto converterror2vector = [=](TH1 *h)
1315 {
1316 std::vector<double> error;
1317 error.reserve(Nbins);
1318 for (int i = 1; i <= h->GetXaxis()->GetNbins(); i++)
1319 for (int j = 1; j <= h->GetYaxis()->GetNbins(); j++)
1320 for (int k = 1; k <= h->GetZaxis()->GetNbins(); k++)
1321 error.emplace_back(h->GetBinError(i, j, k));
1322 return error;
1323 };
1324
1325 std::vector<double> nominal_content = convert2vector(Nominal);
1326 std::vector<double> nominal_error = converterror2vector(Nominal);
1327
1328 int systcount = 0;
1329 for (auto &indices : this->SystematicVariationIndices)
1330 {
1331 double corr = CorrelationFactors[indices[0]];
1332 SystPolicy policy = SystematicVariationPolicies[systcount];
1333 std::vector<std::vector<double>> syst_contents;
1334 for (auto &index : indices)
1335 syst_contents.emplace_back(convert2vector(SystematicVariations[index]));
1336
1337 std::vector<double> vdeltai(indices.size());
1338 std::vector<double> vdeltaj(indices.size());
1339 std::vector<double> vsystcontenti(indices.size());
1340 std::vector<double> vsystcontentj(indices.size());
1341
1342 for (int i = 0; i < Nbins; i++)
1343 {
1344 for (int j = 0; j < Nbins; j++)
1345 {
1346 double deltai, deltaj;
1347 int signi = syst_contents[0][i] - nominal_content[i] > 0 ? 1 : -1;
1348 int signj = syst_contents[0][j] - nominal_content[j] > 0 ? 1 : -1;
1349 for (int k = 0; k < indices.size(); k++)
1350 {
1351 vdeltai[k] = fabs(syst_contents[k][i] - nominal_content[i]);
1352 vdeltaj[k] = fabs(syst_contents[k][j] - nominal_content[j]);
1353 vsystcontenti[k] = syst_contents[k][i];
1354 vsystcontentj[k] = syst_contents[k][j];
1355 }
1356
1357 if (policy == SystPolicy::Maximum)
1358 {
1359 deltai = *std::max_element(vdeltai.begin(), vdeltai.end());
1360 deltaj = *std::max_element(vdeltaj.begin(), vdeltaj.end());
1361 }
1362 else if (policy == SystPolicy::Average)
1363 {
1364 deltai = std::accumulate(vdeltai.begin(), vdeltai.end(), (double)0) / indices.size();
1365 deltaj = std::accumulate(vdeltaj.begin(), vdeltaj.end(), (double)0) / indices.size();
1366 }
1367 else if (policy == SystPolicy::Variance)
1368 {
1369 deltai = std::sqrt(Variance(vsystcontenti));
1370 deltaj = std::sqrt(Variance(vsystcontentj));
1371 }
1372 else if (policy == SystPolicy::CTEQ)
1373 {
1374 deltai = (vsystcontenti[0] - vsystcontenti[1]) / 2;
1375 deltaj = (vsystcontentj[0] - vsystcontentj[1]) / 2;
1376 signi = 1;
1377 signj = 1;
1378 }
1379
1380 if (i != j)
1381 cov[i][j] += signi * signj * deltai * deltaj * corr;
1382 else
1383 cov[i][j] += signi * signj * deltai * deltaj;
1384 }
1385 }
1386
1387 systcount++;
1388 }
1389
1390 // add nominal stat unc.
1391 for (int i = 0; i < Nbins; i++)
1392 cov[i][i] += nominal_error[i] * nominal_error[i];
1393
1394 // here only syst unc are considered
1395 histtool->ConvertVectorToTH2D(Covariance, cov);
1396 histtool->ConvertCovToCorr(Covariance, Correlation);
1397 }
1398 else if (SystematicVariations.size() > 0 && Nbins > 1e3)
1399 {
1400 MSGUser()->MSG_WARNING("NGHist ", GetResultName(), " has too many bins, can not effectively generate a covariance matrix");
1401 }
1402
1403 Processed = true;
1404
1405 return true;
1406 }
1407
1411 template <typename HistType>
1412 inline bool NGHist<HistType>::Recover(TFile *file)
1413 {
1414 auto titleguard = MSGUser()->StartTitleWithGuard("NGHist::Recover");
1415 if (!file)
1416 {
1417 MSGUser()->MSG_ERROR("input file does not exist for plot ", GetResultName());
1418 return false;
1419 }
1420
1421 HistType *nominal_backup = nullptr;
1422 if (Nominal)
1423 {
1424 nominal_backup = Nominal;
1425 Nominal = nullptr;
1426 }
1427
1428 Nominal = (HistType *)file->Get(GetResultName());
1429
1430 if (!Nominal)
1431 {
1432 Nominal = nominal_backup;
1433 MSGUser()->MSG_ERROR("Hist ", GetResultName(), " does not exist in file ", file->GetName(), " will keep the original one.");
1434 return false;
1435 }
1436 else
1437 {
1438 Nominal->SetDirectory(0);
1439 delete nominal_backup;
1440 }
1441
1442 FillHist = Nominal;
1443
1444 TTree *infotree = (TTree *)file->Get(GetResultName() + "_SystInfoTree");
1445 if (infotree != nullptr)
1446 {
1447 std::vector<std::string> *systnames = nullptr;
1448 int policy;
1449 double corrfactor;
1450 int index;
1451 bool doindexcheck = false;
1452
1453 infotree->SetBranchAddress("SystNames", &systnames);
1454 infotree->SetBranchAddress("Policy", &policy);
1455 infotree->SetBranchAddress("CorrFactor", &corrfactor);
1456 // for old version
1457 if (infotree->SetBranchAddress("Index", &index) != TTree::ESetBranchAddressStatus::kMissingBranch)
1458 doindexcheck = true;
1459
1460 for (int i = 0; i < infotree->GetEntries(); i++)
1461 {
1462 infotree->GetEntry(i);
1463 // for avoid the case that hadd auto adds syst-into tree
1464 if (doindexcheck && index != i)
1465 break;
1466 std::vector<TString> tstring_systnames;
1467 for (int j = 0; j < systnames->size(); j++)
1468 tstring_systnames.emplace_back(systnames->at(j));
1469
1470 this->BookSystematicVariation(tstring_systnames, corrfactor, static_cast<SystPolicy>(policy));
1471 }
1472 }
1473
1474 for (size_t i = 0; i < SystematicVariationNames.size(); i++)
1475 {
1476 TString SystHistName = GetResultName() + "_Syst_" + SystematicVariationNames[i];
1477 if (SystematicVariations[i])
1478 {
1479 delete SystematicVariations[i];
1480 SystematicVariations[i] = nullptr;
1481 }
1482 SystematicVariations[i] = (HistType *)file->Get(SystHistName.TString::Data());
1483 if (!SystematicVariations[i])
1484 MSGUser()->MSG_ERROR("Hist ", SystHistName, " does not exist");
1485 else
1486 SystematicVariations[i]->SetDirectory(0);
1487 }
1488
1489 // get final, correlation and covariance if possible
1490 if (Final)
1491 {
1492 delete Final;
1493 Final = nullptr;
1494 }
1495 if (Covariance)
1496 {
1497 delete Covariance;
1498 Covariance = nullptr;
1499 }
1500 if (Correlation)
1501 {
1502 delete Correlation;
1503 Correlation = nullptr;
1504 }
1505 Final = (HistType *)file->Get("Final_" + GetResultName());
1506 Covariance = (TH2D *)file->Get("Covariance_" + GetResultName());
1507 Correlation = (TH2D *)file->Get("Correlation_" + GetResultName());
1508
1509 if (Final)
1510 Final->SetDirectory(0);
1511 if (Covariance)
1512 Covariance->SetDirectory(0);
1513 if (Correlation)
1514 Correlation->SetDirectory(0);
1515
1516 FillHist = Nominal;
1517 Processed = true;
1518
1519 return true;
1520 }
1521
1525 template <typename HistType>
1526 inline bool NGHist<HistType>::Recover(std::shared_ptr<TFileHelper> infile)
1527 {
1528 auto titleguard = MSGUser()->StartTitleWithGuard("NGHist::Recover");
1529
1530 if (!infile)
1531 {
1532 MSGUser()->MSG_ERROR("input file does not exist for plot ", GetResultName());
1533 return false;
1534 }
1535
1536 HistType *nominal_backup = nullptr;
1537 if (Nominal)
1538 {
1539 nominal_backup = Nominal;
1540 Nominal = nullptr;
1541 }
1542
1543 Nominal = (HistType *)infile->Get(GetResultName().Data());
1544
1545 if (!Nominal)
1546 {
1547 Nominal = nominal_backup;
1548 return false;
1549 }
1550 else
1551 {
1552 Nominal->SetDirectory(0);
1553 delete nominal_backup;
1554 }
1555
1556 FillHist = Nominal;
1557
1558 TTree *infotree = infile->GetTree((GetResultName() + "_SystInfoTree").Data());
1559 if (infotree != nullptr)
1560 {
1561 std::vector<std::string> *systnames = nullptr;
1562 int policy;
1563 double corrfactor;
1564 int index;
1565 bool doindexcheck = false;
1566
1567 infotree->SetBranchAddress("SystNames", &systnames);
1568 infotree->SetBranchAddress("Policy", &policy);
1569 infotree->SetBranchAddress("CorrFactor", &corrfactor);
1570 // for old version
1571 if (infotree->SetBranchAddress("Index", &index) != TTree::ESetBranchAddressStatus::kMissingBranch)
1572 doindexcheck = true;
1573
1574 for (int i = 0; i < infotree->GetEntries(); i++)
1575 {
1576 infotree->GetEntry(i);
1577 // for avoid the case that hadd auto adds syst-into tree
1578 if (doindexcheck && index != i)
1579 break;
1580
1581 std::vector<TString> tstring_systnames;
1582 for (int j = 0; j < systnames->size(); j++)
1583 tstring_systnames.emplace_back(systnames->at(j));
1584
1585 this->BookSystematicVariation(tstring_systnames, corrfactor, static_cast<SystPolicy>(policy));
1586 }
1587
1588 delete systnames;
1589 delete infotree;
1590 }
1591
1592 for (size_t i = 0; i < SystematicVariationNames.size(); i++)
1593 {
1594 TString SystHistName = GetResultName() + "_Syst_" + SystematicVariationNames[i];
1595 if (SystematicVariations[i])
1596 {
1597 delete SystematicVariations[i];
1598 SystematicVariations[i] = nullptr;
1599 }
1600
1601 SystematicVariations[i] = (HistType *)infile->Get(SystHistName.Data());
1602 }
1603
1604 // get final, correlation and covariance if possible
1605 if (Final)
1606 {
1607 delete Final;
1608 Final = nullptr;
1609 }
1610 if (Covariance)
1611 {
1612 delete Covariance;
1613 Covariance = nullptr;
1614 }
1615 if (Correlation)
1616 {
1617 delete Correlation;
1618 Correlation = nullptr;
1619 }
1620
1621 TString finalname = "Final_" + GetResultName();
1622 Final = (HistType *)infile->Get(finalname.Data());
1623
1624 if (Nominal)
1625 {
1626 if (Nominal->GetNbinsX() * Nominal->GetNbinsY() * Nominal->GetNbinsZ() <= 1e3 && SystematicVariations.size() > 0)
1627 {
1628 TString covname = "Covariance_" + GetResultName();
1629 TString corrname = "Correlation_" + GetResultName();
1630 Covariance = (TH2D *)infile->Get(covname.Data());
1631 Correlation = (TH2D *)infile->Get(corrname.Data());
1632 }
1633 }
1634
1635 FillHist = Nominal;
1636 Processed = true;
1637
1638 return true;
1639 }
1640
1644 template <typename HistType>
1645 inline bool NGHist<HistType>::Recover(const TString &filename)
1646 {
1647 TFile *file = new TFile(filename.TString::Data());
1648 if (file == nullptr)
1649 {
1650 auto guard = MSGUser()->StartTitleWithGuard("NGHist::Recover");
1651 MSGUser()->MSG_ERROR("ROOT file ", filename, " does not exist for plot ", GetResultName());
1652 return false;
1653 }
1654
1655 auto returnvalue = Recover(file);
1656 file->Close();
1657
1658 return returnvalue;
1659 }
1660
1663 template <typename HistType>
1664 inline bool NGHist<HistType>::Recover()
1665 {
1666 return this->Recover(this->GetOutputFileName());
1667 }
1668
1672 template <typename HistType>
1673 inline bool NGHist<HistType>::RegroupSystematicVariation(const std::vector<TString> &names, SystPolicy policy)
1674 {
1675 auto titleguard = MSGUser()->StartTitleWithGuard("NGHist::RegroupSystematicVariation");
1676 if (names.size() != 2 && policy == SystPolicy::CTEQ)
1677 {
1678 MSGUser()->MSG_ERROR("The CTEQ policy only applies to two systematic variations at the same time for plot ", GetResultName());
1679 return false;
1680 }
1681
1682 this->SetUnprocessed();
1683
1684 std::vector<size_t> indices;
1685 for (auto &name : names)
1686 {
1687 auto findresult = SystematicVariationBackup.find(name);
1688 if (findresult != SystematicVariationBackup.end())
1689 {
1690 indices.emplace_back(findresult->second);
1691 SystematicVariationBackup.erase(findresult);
1692 }
1693 else
1694 MSGUser()->MSG_WARNING("Syst ", name, " is not in the backup list of plot ", GetResultName());
1695 }
1696
1697 if (indices.size() > 0)
1698 {
1699 SystematicVariationIndices.emplace_back(indices);
1700 SystematicVariationPolicies.emplace_back(policy);
1701 }
1702 else
1703 {
1704 MSGUser()->MSG_ERROR("No suitable variation found in plot ", GetResultName());
1705 return false;
1706 }
1707
1708 // check corr value same or not
1709 if (indices.size() > 1)
1710 {
1711 auto corr = CorrelationFactors[indices[0]];
1712 for (auto &index : indices)
1713 if (CorrelationFactors[index] != corr)
1714 MSGUser()->MSG_WARNING("Correlation factors are not equal for input variations, will take the first one for plot ", GetResultName());
1715 }
1716
1717 return true;
1718 }
1719
1722 template <typename HistType>
1723 inline bool NGHist<HistType>::RemoveSystematicVariation(const TString &name)
1724 {
1725 // move the input variation to backup
1726 auto guard = MSGUser()->StartTitleWithGuard("NGHist::RemoveSystematicVariation");
1727
1728 if (auto findresult_backup = SystematicVariationBackup.find(name); findresult_backup != SystematicVariationBackup.end())
1729 {
1730 MSGUser()->MSG_WARNING("variation ", name, " has already been in the backup list of plot ", GetResultName());
1731 return false;
1732 }
1733
1734 if (auto findresult = this->SystematicIndexMap.find(name); findresult != SystematicIndexMap.end())
1735 {
1736 for (int i = 0; i < SystematicVariationIndices.size(); i++)
1737 {
1738 if (std::find(SystematicVariationIndices[i].begin(), SystematicVariationIndices[i].end(), findresult->second) != SystematicVariationIndices[i].end())
1739 {
1740 for (auto &index : SystematicVariationIndices[i])
1741 SystematicVariationBackup[SystematicVariationNames[index]] = index;
1742
1743 SystematicVariationIndices.erase(SystematicVariationIndices.begin() + i);
1744 SystematicVariationPolicies.erase(SystematicVariationPolicies.begin() + i);
1745
1746 break;
1747 }
1748 }
1749
1750 return true;
1751 }
1752 else
1753 {
1754 MSGUser()->MSG_WARNING("variation ", name, " does not exist in plot ", GetResultName());
1755 return false;
1756 }
1757 }
1758
1762 template <typename HistType>
1763 inline bool NGHist<HistType>::RenameSystematicVariation(const TString &name_old, const TString &name_new)
1764 {
1765 auto guard = MSGUser()->StartTitleWithGuard("NGHist::RenameSystematicVariation");
1766
1767 if (auto findresult = this->SystematicIndexMap.find(name_old); findresult != SystematicIndexMap.end())
1768 {
1769 SystematicIndexMap[name_new] = findresult->second;
1770 SystematicVariationNames[findresult->second] = name_new;
1771 TString SystHistName = GetResultName() + "_Syst_" + name_new;
1772 SystematicVariations[findresult->second]->SetName(SystHistName.Data());
1773 SystematicVariations[findresult->second]->SetTitle(SystHistName.Data());
1774 SystematicIndexMap.erase(findresult);
1775
1776 if (auto findresult2 = this->SystematicVariationBackup.find(name_old); findresult2 != SystematicVariationBackup.end())
1777 {
1778 SystematicVariationBackup[name_new] = findresult2->second;
1779 SystematicVariationBackup.erase(findresult2);
1780 }
1781
1782 return true;
1783 }
1784 else
1785 {
1786 MSGUser()->MSG_WARNING("variation ", name_old, " does not exist in plot ", GetResultName());
1787 return false;
1788 }
1789 }
1790
1792 template <typename HistType>
1793 inline void NGHist<HistType>::Reset()
1794 {
1795 Nominal->Reset("ICESM");
1796 for (auto &sh : SystematicVariations)
1797 sh->Reset("ICESM");
1798
1799 if (Final)
1800 Final->Reset("ICESM");
1801 if (Covariance)
1802 Covariance->Reset("ICESM");
1803 if (Correlation)
1804 Correlation->Reset("ICESM");
1805
1806 ClearLinkPlot();
1807 SetUnprocessed();
1808 }
1809
1812 template <typename HistType>
1813 inline void NGHist<HistType>::Scale(double sf)
1814 {
1815 if (Nominal)
1816 Nominal->Scale(sf);
1817 if (Final)
1818 Final->Scale(sf);
1819 if (Covariance)
1820 Covariance->Scale(sf);
1821
1822 for (size_t i = 0; i < SystematicVariationNames.size(); i++)
1823 if (SystematicVariations[i])
1824 SystematicVariations[i]->Scale(sf);
1825 }
1826
1832 template <typename HistType>
1833 inline void NGHist<HistType>::SetLinkPlot(uint64_t index, std::shared_ptr<HistBase> subplot)
1834 {
1835 if (subplot != nullptr)
1836 {
1837 SetUnprocessed();
1838 if (LinkedPlots.size() < (index + 1))
1839 {
1840 int num_needed = index + 1 - LinkedPlots.size();
1841 for (int i = 0; i < num_needed; i++)
1842 LinkedPlots.emplace_back(nullptr);
1843 }
1844
1845 LinkedPlots[index] = subplot;
1846 }
1847 else
1848 {
1849 auto guard = MSGUser()->StartTitleWithGuard("NGHist::SetLinkPlot");
1850 MSGUser()->MSG_WARNING("for plot ", GetResultName(), ", given link plot ", index, " is null");
1851 return;
1852 }
1853 }
1854
1867 template <typename HistType>
1868 inline void NGHist<HistType>::SetLinkType(const TString &type)
1869 {
1870 SetUnprocessed();
1871 LinkType = type;
1872 }
1873
1909 template <typename HistType>
1910 inline void NGHist<HistType>::SetLinkType(std::function<double(const std::vector<double> &)> contentfunction, std::function<double(const std::vector<double> &, const std::vector<double> &)> errorfunction)
1911 {
1912 SetUnprocessed();
1913 LinkContentFunction = contentfunction;
1914 LinkErrorFunction = errorfunction;
1915 }
1916
1936 template <typename HistType>
1937 inline void NGHist<HistType>::SetLinkType(std::function<void(const std::vector<TH1 *> &, TH1 *)> f)
1938 {
1939 SetUnprocessed();
1940 LinkHistFunction = f;
1941 }
1942
1944 template <typename HistType>
1945 inline bool NGHist<HistType>::SetSystematicVariation(const TString &name)
1946 {
1947 if (name == "Nominal")
1948 {
1949 FillHist = Nominal;
1950 return true;
1951 }
1952 auto findresult = SystematicIndexMap.find(name);
1953 if (findresult != SystematicIndexMap.end())
1954 {
1955 FillHist = SystematicVariations[findresult->second];
1956 return true;
1957 }
1958 else
1959 {
1960 auto guard = MSGUser()->StartTitleWithGuard("NGHist::SetSystematicVariation");
1961 MSGUser()->MSG_WARNING("NGHist ", GetResultName(), " does not have variation ", name, ", will set to nominal");
1962 FillHist = Nominal;
1963 return false;
1964 }
1965 }
1966
1968 template <typename HistType>
1969 inline bool NGHist<HistType>::SetSystematicVariation(uint64_t index)
1970 {
1971 if (index == 0)
1972 {
1973 FillHist = Nominal;
1974 return true;
1975 }
1976 else if (index > 0 && index <= SystematicVariations.size())
1977 {
1978 FillHist = SystematicVariations[index - 1];
1979 return true;
1980 }
1981 else
1982 {
1983 auto guard = MSGUser()->StartTitleWithGuard("NGHist::SetSystematicVariation");
1984 MSGUser()->MSG_WARNING("NGHist ", GetResultName(), " does not have variation ", index, ", will set to nominal");
1985 FillHist = Nominal;
1986 return false;
1987 }
1988 }
1989
1991 template <typename HistType>
1992 inline void NGHist<HistType>::SetUnprocessed()
1993 {
1994 Processed = false;
1995 if (Final)
1996 {
1997 delete Final;
1998 Final = nullptr;
1999 }
2000
2001 if (Covariance)
2002 {
2003 delete Covariance;
2004 Covariance = nullptr;
2005 }
2006
2007 if (Correlation)
2008 {
2009 delete Correlation;
2010 Correlation = nullptr;
2011 }
2012 }
2013
2021 template <typename HistType>
2022 inline void NGHist<HistType>::Write()
2023 {
2024 Nominal->Write();
2025 for (size_t i = 0; i < SystematicVariations.size(); i++)
2026 SystematicVariations[i]->Write();
2027
2028 if (Final != nullptr)
2029 Final->Write();
2030 if (Covariance != nullptr)
2031 Covariance->Write();
2032 if (Correlation != nullptr)
2033 Correlation->Write();
2034
2035 // write syst info for fast recover
2036 if (SystematicVariations.size() > 0)
2037 {
2038 TTree *infotree = new TTree(GetResultName() + "_SystInfoTree", GetResultName() + "_SystInfoTree");
2039 std::vector<std::string> systnames;
2040 int policy;
2041 double corrfactor;
2042 int index = 0;
2043 infotree->Branch("SystNames", &systnames);
2044 infotree->Branch("Policy", &policy);
2045 infotree->Branch("CorrFactor", &corrfactor);
2046 infotree->Branch("Index", &index);
2047
2048 for (size_t i = 0; i < SystematicVariationIndices.size(); i++)
2049 {
2050 systnames.clear();
2051 policy = static_cast<std::underlying_type<SystPolicy>::type>(SystematicVariationPolicies[i]);
2052 corrfactor = CorrelationFactors[SystematicVariationIndices[i][0]];
2053 for (size_t j = 0; j < SystematicVariationIndices[i].size(); j++)
2054 systnames.emplace_back(SystematicVariationNames[SystematicVariationIndices[i][j]]);
2055
2056 infotree->Fill();
2057 ++index;
2058 }
2059
2060 // write backup
2061 for (auto &iter : SystematicVariationBackup)
2062 {
2063 systnames.clear();
2064 policy = static_cast<std::underlying_type<SystPolicy>::type>(SystPolicy::Backup);
2065 corrfactor = CorrelationFactors[iter.second];
2066 systnames.emplace_back(iter.first.Data());
2067 infotree->Fill();
2068 ++index;
2069 }
2070
2071 infotree->Write();
2072 delete infotree;
2073 }
2074 }
2075
2077 template <typename HistType>
2078 inline void NGHist<HistType>::WriteToFile()
2079 {
2080 auto titleguard = MSGUser()->StartTitleWithGuard("NGHist::WriteToFile");
2081 Process();
2082 if (GetOutputFileName() == "")
2083 {
2084 MSGUser()->MSG_WARNING("input an empty output file name for plot ", GetResultName(), " will be ignored");
2085 return;
2086 }
2087 TFile *file = new TFile(GetOutputFileName(), "UPDATE");
2088 file->cd();
2089 Write();
2090 file->Close();
2091 }
2092}
NAGASH::HistBase
Virtual base class for histograms.
Definition HistBase.h:27
NAGASH::HistBase::SystPolicy
SystPolicy
Policies of dealing with systematic uncertainties.
Definition HistBase.h:33
NAGASH::HistBase::SystPolicy::Backup
@ Backup
not add to total systematic uncertainties, just for future use.
NAGASH::HistBase::SystPolicy::Maximum
@ Maximum
take the maximum difference in the same group.
NAGASH::NGHist
NAGASH interface for using ROOT histograms.
Definition NGHist.h:18
NAGASH::NGHist::RenameSystematicVariation
bool RenameSystematicVariation(const TString &name_old, const TString &name_new) override
Rename a systematic variation.
Definition NGHist.h:1763
NAGASH::NGHist::GetBinError
double GetBinError(Args &&...args)
Get the bin error of the current variation using TH1::GetBinError().
Definition NGHist.h:106
NAGASH::NGHist::Rebin
std::shared_ptr< NGHist< HistType > > Rebin(const TString &name, int rebinnum)
Rebin 1D histogram with the given rebin number, store the result in a new NGHist. Note that the link ...
Definition NGHist.h:531
NAGASH::NGHist::RemoveSystematicVariation
bool RemoveSystematicVariation(const TString &name) override
Move a systematic variation to backup the variations in the same group will also be moved to the back...
Definition NGHist.h:1723
NAGASH::NGHist::SetUnprocessed
void SetUnprocessed() override
Set this NGHist to a unprocessed state.
Definition NGHist.h:1992
NAGASH::NGHist::LinkType
TString LinkType
Definition NGHist.h:40
NAGASH::NGHist::SystematicVariationNames
std::vector< TString > SystematicVariationNames
Definition NGHist.h:26
NAGASH::NGHist::GetFillHist
HistType * GetFillHist()
Get the pointer to the current variation histogram.
Definition NGHist.h:57
NAGASH::NGHist::SystematicVariations
std::vector< HistType * > SystematicVariations
Definition NGHist.h:27
NAGASH::NGHist::GetVariationVirtual
TH1 * GetVariationVirtual(const TString &name) override
Virtual function to get the variation histogram.
Definition NGHist.h:1063
NAGASH::NGHist::SetLinkType
void SetLinkType(const TString &type) override
Set the type of the link.
Definition NGHist.h:1868
NAGASH::NGHist::GetBinContent
double GetBinContent(Args &&...args)
Get the bin content of the current variation using TH1::GetBinContent().
Definition NGHist.h:103
NAGASH::NGHist::GetVariation
HistType * GetVariation(const TString &name)
Get the pointer to the variation histogram based on the name.
Definition NGHist.h:164
NAGASH::NGHist::GetCorrelation
TH2D * GetCorrelation()
Get the pointer to the correlation histogram.
Definition NGHist.h:85
NAGASH::NGHist::LinkedPlots
std::vector< std::shared_ptr< HistBase > > LinkedPlots
Definition NGHist.h:39
NAGASH::NGHist::SystematicIndexMap
std::map< TString, int > SystematicIndexMap
Definition NGHist.h:25
NAGASH::NGHist::Nominal
HistType * Nominal
Definition NGHist.h:23
NAGASH::NGHist::LinkErrorFunction
std::function< double(const std::vector< double > &, const std::vector< double > &)> LinkErrorFunction
Definition NGHist.h:42
NAGASH::NGHist::Rebin2D
std::shared_ptr< NGHist< HistType > > Rebin2D(const TString &name, int rebinnum_x, int rebinum_y)
Rebin 2D histogram with the given rebin numbers, store the result in a new NGHist....
Definition NGHist.h:566
NAGASH::NGHist::Combine
void Combine(std::shared_ptr< Result > result) override
Combine the current NGHist with another NGHist.
Definition NGHist.h:1031
NAGASH::NGHist::Scale
void Scale(double) override
Scale all histograms using TH1::Scale().
Definition NGHist.h:1813
NAGASH::NGHist::SystematicVariationIndices
std::vector< std::vector< size_t > > SystematicVariationIndices
Definition NGHist.h:29
NAGASH::NGHist::GetBinNominalError
double GetBinNominalError(Args &&...args)
Get the bin error of the nominal histogram using TH1::GetBinError().
Definition NGHist.h:109
NAGASH::NGHist::LinkContentFunction
std::function< double(const std::vector< double > &)> LinkContentFunction
Definition NGHist.h:41
NAGASH::NGHist::~NGHist
virtual ~NGHist()
Destructor.
Definition NGHist.h:213
NAGASH::NGHist::GetVariation
HistType * GetVariation(uint64_t index)
Get the pointer to the variation histogram based on the index.
Definition NGHist.h:182
NAGASH::NGHist::GetBinTotalError
double GetBinTotalError(Args &&...args)
Get the total bin error of the final histogram using TH1::GetBinError().
Definition NGHist.h:113
NAGASH::NGHist::GetVariationTypes
std::vector< std::tuple< std::vector< TString >, double, SystPolicy > > GetVariationTypes() override
Get the structure of the systematic variations.
Definition NGHist.h:1085
NAGASH::NGHist::Recover
bool Recover() override
Recover the NGHist from a ROOT file with the name passed to the constructor.
Definition NGHist.h:1664
NAGASH::NGHist::GetPartialCovariance
std::vector< std::vector< double > > GetPartialCovariance(const std::vector< TString > &systlist)
Get the covariance matrix for the given list of systematic variations.
Definition NGHist.h:386
NAGASH::NGHist::Clone
std::shared_ptr< NGHist< HistType > > Clone(const TString &name)
Clone the current NGHist with the given name.
Definition NGHist.h:1024
NAGASH::NGHist::FillHist
HistType * FillHist
Definition NGHist.h:22
NAGASH::NGHist::Final
HistType * Final
Definition NGHist.h:33
NAGASH::NGHist::GetVariationCorrelationFactor
double GetVariationCorrelationFactor(const TString &name) override
Get the correlation factor for the systematic variation with the given name.
Definition NGHist.h:198
NAGASH::NGHist::NGHist
NGHist(std::shared_ptr< MSGTool > MSG, std::shared_ptr< ConfigTool > c, const TString &name, const TString &filename, const Args &...args)
Constructor.
Definition NGHist.h:355
NAGASH::NGHist::GetBinSystematicError
double GetBinSystematicError(Args &&...args)
Get the total systematic error of ith bin.
Definition NGHist.h:133
NAGASH::NGHist::Variance
static double Variance(const std::vector< double > &v)
Calculate the variance of the given vector.
Definition NGHist.h:282
NAGASH::NGHist::Reset
void Reset() override
Reset all histograms.
Definition NGHist.h:1793
NAGASH::NGHist::LinkHistFunction
std::function< void(const std::vector< TH1 * > &, TH1 *)> LinkHistFunction
Definition NGHist.h:43
NAGASH::NGHist::Processed
bool Processed
Definition NGHist.h:36
NAGASH::NGHist::RegroupSystematicVariation
bool RegroupSystematicVariation(const std::vector< TString > &names, SystPolicy policy) override
Regroup the backup variations with the give name of policy.
Definition NGHist.h:1673
NAGASH::NGHist::GetCovariance
TH2D * GetCovariance()
Get the pointer to the covariance histogram.
Definition NGHist.h:74
NAGASH::NGHist::CorrelationFactors
std::vector< double > CorrelationFactors
Definition NGHist.h:28
NAGASH::NGHist::GetBinFinalError
double GetBinFinalError(Args &&...args)
Get the total bin error of the final histogram using TH1::GetBinError().
Definition NGHist.h:123
NAGASH::NGHist::SetBinError
void SetBinError(Args &&...args)
Set bin error of the current variation using TH1::SetBinError().
Definition NGHist.h:100
NAGASH::NGHist::WriteToFile
void WriteToFile() override
Write the NGHist to a ROOT file.
Definition NGHist.h:2078
NAGASH::NGHist::GetNVariations
size_t GetNVariations() override
Get the number of systematic variations.
Definition NGHist.h:60
NAGASH::NGHist::Correlation
TH2D * Correlation
Definition NGHist.h:35
NAGASH::NGHist::Covariance
TH2D * Covariance
Definition NGHist.h:34
NAGASH::NGHist::SetLinkPlot
void SetLinkPlot(uint64_t index, std::shared_ptr< HistBase > subplot) override
Link this NGHist to anothe one.
Definition NGHist.h:1833
NAGASH::NGHist::GetVariationName
TString GetVariationName(uint64_t index) override
Get the name of the ith systematic variation.
Definition NGHist.h:146
NAGASH::NGHist::Process
bool Process() override
process the NGHist.
Definition NGHist.h:1155
NAGASH::NGHist::Write
void Write() override
Write the NGHist to a ROOT file.
Definition NGHist.h:2022
NAGASH::NGHist::CloneVirtual
std::shared_ptr< HistBase > CloneVirtual(const TString &name) override
Definition NGHist.h:305
NAGASH::NGHist::Fill
void Fill(Args &&...args)
Fill the current variation using TH1::Fill().
Definition NGHist.h:51
NAGASH::NGHist::GetFinalHist
HistType * GetFinalHist()
Get the pointer to the final histogram, whose contents are nominal contents and errors are the total ...
Definition NGHist.h:63
NAGASH::NGHist::SetSystematicVariation
bool SetSystematicVariation(const TString &name="Nominal") override
Set the NGHist to the give variation name.
Definition NGHist.h:1945
NAGASH::NGHist::SystematicVariationBackup
std::map< TString, int > SystematicVariationBackup
Definition NGHist.h:31
NAGASH::NGHist::GetNominalHist
HistType * GetNominalHist()
Get the pointer to the nominal variation histogram.
Definition NGHist.h:54
NAGASH::NGHist::Rebin3D
std::shared_ptr< NGHist< HistType > > Rebin3D(const TString &name, int rebinnum_x, int rebinum_y, int rebinum_z)
Rebin 3D histogram with the given rebin numbers, store the result in a new NGHist....
Definition NGHist.h:602
NAGASH::NGHist::ClearLinkPlot
void ClearLinkPlot() override
Clear all linked plots.
Definition NGHist.h:1013
NAGASH::NGHist::IsBackupVariation
bool IsBackupVariation(uint64_t) override
If ith variation is a backup variation.
Definition NGHist.h:1113
NAGASH::NGHist::GetNominalHistVirtual
TH1 * GetNominalHistVirtual() override
Virtual function to get the nominal histogram.
Definition NGHist.h:1077
NAGASH::NGHist::SetBinContent
void SetBinContent(Args &&...args)
Set bin content of the current variation using TH1::SetBinContent().
Definition NGHist.h:97
NAGASH::NGHist::SystematicVariationPolicies
std::vector< SystPolicy > SystematicVariationPolicies
Definition NGHist.h:30
NAGASH::NGHist::GetVariationNames
std::vector< TString > GetVariationNames()
Get the vector of the name of all systematic variations.
Definition NGHist.h:157
NAGASH::NGHist::BookSystematicVariation
bool BookSystematicVariation(const TString &name, double corr_factor=0, SystPolicy policy=SystPolicy::Maximum) override
Book a systematic variation.
Definition NGHist.h:924
NAGASH::PlotGroup
Provide a base class for manipulating a group of histograms at the same time.
Definition PlotGroup.h:18
NAGASH::Result::GetResultName
const TString & GetResultName()
Definition Result.h:70
NAGASH::Result::MSGUser
std::shared_ptr< MSGTool > MSGUser()
Return the internal MSGTool.
Definition Result.h:106
NAGASH::Uncertainty::AplusB
static double AplusB(double A_err, double B_err)
Calculate the uncertainty of .
Definition Uncertainty.cxx:13
Generated on Mon Sep 2 2024 14:24:32 for NAGASH by doxygen 1.9.8