TEfficiencyCalibrator.cxx

Go to the documentation of this file.

#include "TEfficiencyCalibrator.h"
 
#if __cplusplus >= 201703L
 
#include <cstdarg>
#include <cstdio>
#include <fcntl.h>
 
#include "TGTableLayout.h"
#include "TTimer.h"
#include "TCanvas.h"
 
#include "TGauss.h"
#include "TABPeak.h"
#include "TAB3Peak.h"
#include "TRWPeak.h"
#include "TRedirect.h"
 
//////////////////////////////////////// TEfficiencyTab ////////////////////////////////////////
TEfficiencyTab::TEfficiencyTab(TEfficiencyDatatypeTab* parent, TNucleus* nucleus, std::tuple<TH1*, TH1*, TH2*> hists, TGCompositeFrame* frame)
   : fFrame(frame), fNucleus(nucleus), fParent(parent),
     fSingles(std::get<0>(hists)), fSummingOut(std::get<1>(hists)), fSummingIn(std::get<2>(hists))
// name *180Corr                 // name *Sum180Corr
{
   if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kSubroutines) {
      std::cout << "Assigned singles, summing out, and summing in histograms (" << fSingles->GetName() << ", " << fSummingOut->GetName() << ", " << fSummingIn->GetName() << ")" << std::endl;
   }
   BuildInterface();
   {
      TRedirect redirect("fitOutput.txt");
      std::cout << "======================================== Finding peaks in " << fSingles->GetName() << ", " << fSummingOut->GetName() << ", and " << fSummingIn->GetName() << std::endl;
      FindPeaks();
   }
}
 
void TEfficiencyTab::BuildInterface()
{
   // top frame with one canvas, status bar, and controls
   fTopFrame         = new TGHorizontalFrame(fFrame, TEfficiencyCalibrator::WindowWidth() / 2, 450);
   fProjectionCanvas = new TRootEmbeddedCanvas("ProjectionCanvas", fTopFrame, TEfficiencyCalibrator::WindowWidth() / 2, 400);
 
   fTopFrame->AddFrame(fProjectionCanvas, new TGLayoutHints(kLHintsLeft | kLHintsTop | kLHintsExpandY | kLHintsExpandX, 2, 2, 2, 2));
 
   fFrame->AddFrame(fTopFrame, new TGLayoutHints(kLHintsTop | kLHintsExpandX, 2, 2, 2, 2));
 
   fStatusBar               = new TGStatusBar(fFrame, TEfficiencyCalibrator::WindowWidth() / 2, 50);   // NOLINT(cppcoreguidelines-prefer-member-initializer)
   std::array<int, 5> parts = {35, 15, 20, 15, 15};
   fStatusBar->SetParts(parts.data(), parts.size());
 
   fFrame->AddFrame(fStatusBar, new TGLayoutHints(kLHintsBottom | kLHintsExpandX, 2, 2, 2, 2));
   if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kLoops) { std::cout << "Done with " << __PRETTY_FUNCTION__ << std::endl; }   // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay)
}
 
void TEfficiencyTab::FindPeaks()
{
   /// Find and fit the peaks in the singles histogram, then checks for each peak how much summing in and summing out happens.
 
   if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kLoops) {
      std::cout << __PRETTY_FUNCTION__ << std::endl;   // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay)
      std::cout << "Using parent " << fParent << " and transition list " << fNucleus->GetTransitionListByEnergy() << std::endl;
      fNucleus->GetTransitionListByEnergy()->Print();
      std::cout << "Using parent " << fParent << std::flush << ", trying to get peak type " << TEfficiencyCalibrator::PeakType() << std::endl;
   }
   fProjectionCanvas->GetCanvas()->cd();
   // loop through all gamma-rays of the source
   // TODO: combine ranges if peak ranges overlap
   // probably need to: change loop to using size and ->At(index), plus extra loop that finds all peaks within the range
   // maybe also change the creation of the peak to either a function or using TClass::New()
   for(int t = 0; t < fNucleus->GetTransitionListByEnergy()->GetSize(); ++t) {
      auto* transition = static_cast<TTransition*>(fNucleus->GetTransitionListByEnergy()->At(t));
      auto  energy     = transition->GetEnergy();
      if(energy > fSingles->GetXaxis()->GetXmax()) {
         // need to check if we also want to reject peaks whose fit range is partially outside the range of the histogram
         std::cout << "Skipping peak at " << energy << " keV of " << fNucleus->GetA() << fNucleus->GetSymbol() << ", maximum range of histogram is " << fSingles->GetXaxis()->GetXmax() << std::endl;
         continue;
      }
      if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kLoops) {
         std::cout << "Fitting peak at " << energy << " keV, using range " << energy - TEfficiencyCalibrator::Range() << ", " << energy + TEfficiencyCalibrator::Range() << " peak type " << TEfficiencyCalibrator::PeakType() << " singles histogram " << fSingles->GetName() << std::endl;
      }
      fPeakFitter.RemoveAllPeaks();
      fPeakFitter.ResetInitFlag();
      fSummingOutPeakFitter.RemoveAllPeaks();
      fSummingOutPeakFitter.ResetInitFlag();
 
      std::vector<TSinglePeak*> peaks;
      peaks.push_back(NewPeak(energy));
      double lastEnergy = energy;
      for(int t2 = t + 1; t2 < fNucleus->GetTransitionListByEnergy()->GetSize(); ++t2) {
         auto* transition2 = static_cast<TTransition*>(fNucleus->GetTransitionListByEnergy()->At(t2));
         auto  energy2     = transition2->GetEnergy();
         if(lastEnergy + TEfficiencyCalibrator::Range() > energy2 - TEfficiencyCalibrator::Range()) {
            if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kSubroutines) {
               std::cout << t << ". peak: Found " << peaks.size() << ". overlapping peak for " << lastEnergy << ", and " << energy2 << " with range " << TEfficiencyCalibrator::Range() << std::endl;
            }
            // peak ranges overlap, so we add this peak to the fit, skip it in the main loop, and update lastEnergy to refer to it instead
            peaks.push_back(NewPeak(energy2));
            ++t;
            lastEnergy = energy2;
         } else {
            // list of transitions is sorted by energy, so all others will be outside the range as well
            if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kSubroutines) {
               std::cout << t << ". peak: Found " << peaks.size() << ". overlapping peaks, stopping here with " << lastEnergy << ", and " << energy2 << " with range " << TEfficiencyCalibrator::Range() << std::endl;
            }
            break;
         }
      }
 
      fPeakFitter.SetRange(energy - TEfficiencyCalibrator::Range(), lastEnergy + TEfficiencyCalibrator::Range());
      for(auto& peak : peaks) {
         fPeakFitter.AddPeak(peak);
      }
      std::cout << "---------------------------------------- Fitting peak(s) at " << energy << "/" << lastEnergy << " from " << energy - TEfficiencyCalibrator::Range() << " to " << energy + TEfficiencyCalibrator::Range() << std::endl;
      fPeakFitter.Print();
      fPeakFitter.PrintParameters();
      fPeakFitter.Fit(fSingles, "retryfit");
      // we also need to fit the summing out peaks, so we create a copy of the vector of peaks to be added to the summing out peak fitter
      // we could also simply create them in the same loop as the other peaks?
      std::vector<TSinglePeak*> summingOutPeaks = peaks;
      fSummingOutPeakFitter.SetRange(energy - TEfficiencyCalibrator::Range(), lastEnergy + TEfficiencyCalibrator::Range());
      // before we add the peaks to the peak fitter, we need to copy them
      for(auto& peak : summingOutPeaks) {
         peak = NewPeak(peak);
         fSummingOutPeakFitter.AddPeak(peak);
      }
      fSummingOutPeakFitter.Fit(fSummingOut, "retryfit");
      for(size_t p = 0; p < peaks.size(); ++p) {
         auto* peak = peaks[p];
         if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kSubroutines) {
            std::cout << "================================================================================" << std::endl;
            std::cout << "Got centroid " << peak->Centroid() << " +- " << peak->CentroidErr() << " (" << 100 * peak->CentroidErr() / peak->Centroid() << " %), area " << peak->Area() << " +- " << peak->AreaErr() << " (" << 100 * peak->AreaErr() / peak->Area() << " %), FWHM " << peak->FWHM() << ", red. chi sq. " << peak->GetReducedChi2() << std::endl;
         }
         if(peak->Area() < TEfficiencyCalibrator::Threshold() || peak->CentroidErr() > peak->Centroid() || peak->AreaErr() > peak->Area() || std::isnan(peak->Centroid()) || std::isnan(peak->Area()) || std::isnan(peak->CentroidErr()) || std::isnan(peak->AreaErr())) {
            if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kSubroutines) {
               std::cout << "Skipping peak at " << energy << " keV with centroid " << peak->Centroid() << " +- " << peak->CentroidErr() << " (" << 100 * peak->CentroidErr() / peak->Centroid() << " %), FWHM " << peak->FWHM() << ", and area " << peak->Area() << " +- " << peak->AreaErr() << " (" << 100 * peak->AreaErr() / peak->Area() << " %)" << std::endl;
               std::cout << "================================================================================" << std::endl;
            }
            fRemovedFitFunctions.push_back(fPeakFitter.GetFitFunction()->Clone(Form("removed_fit_%.1f", peak->Centroid())));
            continue;
         }
         if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kSubroutines) {
            std::cout << "================================================================================" << std::endl;
         }
         // increase number of points of fit function
         fPeakFitter.GetFitFunction()->SetNpx(1000);
         peak->GetPeakFunction()->SetNpx(1000);
         double fwhm        = peak->FWHM();
         double centroid    = peak->Centroid();
         double centroidErr = peak->CentroidErr();
         if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kLoops) {
            std::cout << "Got centroid " << centroid << " keV, FWHM " << fwhm << ":" << std::endl;
         }
         if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kSubroutines) {
            std::cout << "Writing '" << Form("%s_%.0fkeV", fSingles->GetName(), centroid) << "' to '" << gDirectory->GetName() << "'" << std::endl;
         }
         fSingles->Write(Form("%s_%.0fkeV", fSingles->GetName(), centroid), TObject::kOverwrite);
         if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kSubroutines) {
            std::cout << "Writing '" << Form("%s_%.0fkeV", fSummingOut->GetName(), centroid) << "' to '" << gDirectory->GetName() << "'" << std::endl;
         }
         fSummingOut->Write(Form("%s_%.0fkeV", fSummingOut->GetName(), centroid), TObject::kOverwrite);
         // for summing in we need to estimate the background and subtract that from the integral
         fSummingInProj.push_back(fSummingIn->ProjectionY(Form("%s_%.0f_to_%.0f", fSummingIn->GetName(), centroid - fwhm / 2., centroid + fwhm / 2.), fSummingIn->GetXaxis()->FindBin(centroid - fwhm / 2.), fSummingIn->GetXaxis()->FindBin(centroid + fwhm / 2.)));
         fSummingInProjBg.push_back(fSummingInProj.back()->ShowBackground(TEfficiencyCalibrator::BgParam()));
         if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kSubroutines) {
            std::cout << "Writing '" << Form("%s_%.0fkeV", fSummingInProj.back()->GetName(), centroid) << "' and '" << Form("%s_%.0fkeV", fSummingInProjBg.back()->GetName(), centroid) << "' to '" << gDirectory->GetName() << "'" << std::endl;
         }
         double summingIn = fSummingInProj.back()->Integral() - fSummingInProjBg.back()->Integral();
         fSummingInProj.back()->Write(Form("%s_%.0fkeV", fSummingInProj.back()->GetName(), centroid), TObject::kOverwrite);
         fSummingInProjBg.back()->Write(Form("%s_%.0fkeV", fSummingInProjBg.back()->GetName(), centroid), TObject::kOverwrite);
 
         double summingOut    = summingOutPeaks[p]->Area();
         double correctedArea = peak->Area() - summingIn + summingOutPeaks[p]->AreaErr();
         // uncertainties for summing in and summing out is sqrt(N) (?)
         double correctedAreaErr = TMath::Sqrt(TMath::Power(peak->AreaErr(), 2) + summingIn + summingOut);
         if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kLoops) {
            std::cout << "Got summingIn " << summingIn << ", summingOut " << summingOut << ", correctedArea " << correctedArea << " +- " << correctedAreaErr << std::endl;
         }
         fPeaks.emplace_back(transition, centroid, centroidErr, correctedArea, correctedAreaErr, peak->Area(), peak->AreaErr(), summingIn, summingOut);
         fFitFunctions.push_back(fPeakFitter.GetFitFunction()->Clone(Form("fit_%.1f", centroid)));
      }
   }
 
   Redraw();
   if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kLoops) {
      std::cout << "Unsorted peak vector:";
      for(const auto& peak : fPeaks) { std::cout << " " << std::get<0>(peak)->GetEnergy(); }
      std::cout << std::endl;
   }
   std::sort(fPeaks.begin(), fPeaks.end(), [](const std::tuple<TTransition*, double, double, double, double, double, double, double, double>& lhs, const std::tuple<TTransition*, double, double, double, double, double, double, double, double>& rhs) -> bool { return std::get<0>(lhs)->GetEnergy() < std::get<0>(rhs)->GetEnergy(); });
   if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kLoops) {
      std::cout << "Sorted peak vector:";
      for(auto& peak : fPeaks) { std::cout << " " << std::get<0>(peak)->GetEnergy(); }
      std::cout << std::endl;
   }
}
 
TSinglePeak* TEfficiencyTab::NewPeak(const double& energy)
{
   switch(TEfficiencyCalibrator::PeakType()) {
   case kRWPeak:
      return new TRWPeak(energy);
      break;
   case kABPeak:
      return new TABPeak(energy);
      break;
   case kAB3Peak:
      return new TAB3Peak(energy);
      break;
   case kGauss:
      return new TGauss(energy);
      break;
   default:
      std::cerr << "Unknow peak type " << TEfficiencyCalibrator::PeakType() << ", defaulting to TRWPeak!" << std::endl;
      return new TRWPeak(energy);
   }
}
 
TSinglePeak* TEfficiencyTab::NewPeak(TSinglePeak* peak)
{
   switch(TEfficiencyCalibrator::PeakType()) {
   case kRWPeak:
      return new TRWPeak(*static_cast<TRWPeak*>(peak));
      break;
   case kABPeak:
      return new TABPeak(*static_cast<TABPeak*>(peak));
      break;
   case kAB3Peak:
      return new TAB3Peak(*static_cast<TAB3Peak*>(peak));
      break;
   case kGauss:
      return new TGauss(*static_cast<TGauss*>(peak));
      break;
   default:
      std::cerr << "Unknow peak type " << TEfficiencyCalibrator::PeakType() << ", defaulting to TRWPeak!" << std::endl;
      return new TRWPeak(*static_cast<TRWPeak*>(peak));
   }
}
 
void TEfficiencyTab::Redraw()
{
   fProjectionCanvas->GetCanvas()->cd();
   fSingles->GetListOfFunctions()->Clear();
   for(auto& fit : fFitFunctions) {
      static_cast<TF1*>(fit)->SetLineColor(2);   // red line color for fits
      fSingles->GetListOfFunctions()->Add(fit);
   }
   if(TEfficiencyCalibrator::ShowRemovedFits()) {
      for(auto& fit : fRemovedFitFunctions) {
         static_cast<TF1*>(fit)->SetLineColor(15);   // grey line color for removed fits
         fSingles->GetListOfFunctions()->Add(fit);
      }
   }
   fSingles->Draw();
   fProjectionCanvas->GetCanvas()->Modified();
}
 
void TEfficiencyTab::MakeConnections()
{
   fProjectionCanvas->GetCanvas()->Connect("ProcessedEvent(Int_t, Int_t, Int_t, TObject*)", "TEfficiencyTab", this, "Status(Int_t, Int_t, Int_t, TObject*)");
}
 
void TEfficiencyTab::Disconnect()
{
   fProjectionCanvas->GetCanvas()->Disconnect("ProcessedEvent(Int_t, Int_t, Int_t, TObject*)", this, "Status(Int_t, Int_t, Int_t, TObject*)");
}
 
void TEfficiencyTab::Status(Int_t, Int_t px, Int_t py, TObject* selected)
{
   fStatusBar->SetText(selected->GetName(), 0);
   fStatusBar->SetText(selected->GetObjectInfo(px, py), 1);
}
 
//////////////////////////////////////// TEfficiencyDatatypeTab ////////////////////////////////////////
TEfficiencyDatatypeTab::TEfficiencyDatatypeTab(TEfficiencyCalibrator* parent, std::vector<TNucleus*> nucleus, std::vector<std::tuple<TH1*, TH1*, TH2*>> hists, TGCompositeFrame* frame, const std::string& dataType, TGHProgressBar* progressBar)
   : fFrame(frame), fNucleus(std::move(nucleus)), fParent(parent), fDataType(dataType)
{
   if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kBasicFlow) {
      std::cout << "======================================== creating tab for data type " << dataType << std::endl;
   }
   // create new subdirectory for this tab (and its channel tabs), but store the old directory
   fMainDirectory = gDirectory;
   fSubDirectory  = gDirectory->mkdir(dataType.c_str());
   if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kSubroutines) {
      std::cout << "switching from gDirectory " << gDirectory->GetName() << " to sub directory " << fSubDirectory;
   }
   fSubDirectory->cd();
   if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kSubroutines) {
      std::cout << " = " << fSubDirectory->GetName() << ", main directory is " << fMainDirectory << " = " << (fMainDirectory == nullptr ? "" : fMainDirectory->GetName()) << std::endl;
   }
 
   // left frame with histograms tabs
   fLeftFrame = new TGVerticalFrame(fFrame, TEfficiencyCalibrator::WindowWidth() / 2, TEfficiencyCalibrator::WindowWidth() / 2);
 
   fDataTab = new TGTab(fLeftFrame, TEfficiencyCalibrator::WindowWidth() / 2, 500);
   fEfficiencyTab.resize(hists.size(), nullptr);
   for(size_t i = 0; i < hists.size(); ++i) {
      if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kSubroutines) { std::cout << i << ": Creating efficiency tab using " << fNucleus[i] << " = " << fNucleus[i]->GetName() << ", " << std::get<0>(hists[i])->GetName() << ", " << std::get<1>(hists[i])->GetName() << ", " << std::get<2>(hists[i])->GetName() << ", " << TEfficiencyCalibrator::Range() << ", " << TEfficiencyCalibrator::Threshold() << ", " << TEfficiencyCalibrator::BgParam() << std::endl; }
      fEfficiencyTab[i] = new TEfficiencyTab(this, fNucleus[i], hists[i], fDataTab->AddTab(fNucleus[i]->GetName()));
      progressBar->Increment(1);
   }
 
   fFittingParameterFrame = new TGGroupFrame(fLeftFrame, "Fitting Parameters", kHorizontalFrame);
   fFittingParameterFrame->SetLayoutManager(new TGTableLayout(fFittingParameterFrame, 2, 4, false, 5));   // rows, columns, not homogenous cell sizes, 5 = interval between frames in pixel
   fRangeLabel     = new TGLabel(fFittingParameterFrame, "Range");
   fRangeEntry     = new TGNumberEntry(fFittingParameterFrame, TEfficiencyCalibrator::Range(), 5, kRangeEntry, TGNumberFormat::EStyle::kNESRealTwo, TGNumberFormat::EAttribute::kNEAPositive);
   fBgParamLabel   = new TGLabel(fFittingParameterFrame, "BG Parameter");
   fBgParamEntry   = new TGNumberEntry(fFittingParameterFrame, TEfficiencyCalibrator::BgParam(), 5, kBgParamEntry, TGNumberFormat::EStyle::kNESRealTwo, TGNumberFormat::EAttribute::kNEAPositive);
   fThresholdLabel = new TGLabel(fFittingParameterFrame, "Threshold");
   fThresholdEntry = new TGNumberEntry(fFittingParameterFrame, TEfficiencyCalibrator::Threshold(), 5, kThresholdEntry, TGNumberFormat::EStyle::kNESRealThree, TGNumberFormat::EAttribute::kNEAPositive);
   fPeakTypeBox    = new TGComboBox(fFittingParameterFrame, kPeakTypeBox);
   fPeakTypeBox->AddEntry("Radware", TEfficiencyTab::EPeakType::kRWPeak);
   fPeakTypeBox->AddEntry("Addback", TEfficiencyTab::EPeakType::kABPeak);
   fPeakTypeBox->AddEntry("Addback3", TEfficiencyTab::EPeakType::kAB3Peak);
   fPeakTypeBox->AddEntry("Gaussian", TEfficiencyTab::EPeakType::kGauss);
   fPeakTypeBox->SetMinHeight(200);
   fPeakTypeBox->Resize(100, TEfficiencyCalibrator::LineHeight());
   fPeakTypeBox->Select(TEfficiencyCalibrator::PeakType());
   fFittingControlGroup = new TGHButtonGroup(fFittingParameterFrame, "");
   fRefitButton         = new TGTextButton(fFittingControlGroup, "Refit this source");
   fRefitAllButton      = new TGTextButton(fFittingControlGroup, "Refit all sources");
 
   fFittingParameterFrame->AddFrame(fRangeLabel, new TGTableLayoutHints(0, 1, 0, 1));
   fFittingParameterFrame->AddFrame(fRangeEntry, new TGTableLayoutHints(1, 2, 0, 1));
   fFittingParameterFrame->AddFrame(fBgParamLabel, new TGTableLayoutHints(2, 3, 0, 1));
   fFittingParameterFrame->AddFrame(fBgParamEntry, new TGTableLayoutHints(3, 4, 0, 1));
   fFittingParameterFrame->AddFrame(fThresholdLabel, new TGTableLayoutHints(0, 1, 1, 2));
   fFittingParameterFrame->AddFrame(fThresholdEntry, new TGTableLayoutHints(1, 2, 1, 2));
   fFittingParameterFrame->AddFrame(fPeakTypeBox, new TGTableLayoutHints(2, 3, 1, 2));
   fFittingParameterFrame->AddFrame(fFittingControlGroup, new TGTableLayoutHints(3, 4, 1, 2));
 
   fLeftFrame->AddFrame(fDataTab, new TGLayoutHints(kLHintsTop | kLHintsExpandX | kLHintsExpandY, 2, 2, 2, 2));
   fLeftFrame->AddFrame(fFittingParameterFrame, new TGLayoutHints(kLHintsBottom | kLHintsExpandX, 2, 2, 2, 2));
 
   // right frame with canvas, status bar, and the degree entry
   fRightFrame              = new TGVerticalFrame(fFrame, TEfficiencyCalibrator::WindowWidth() / 2, TEfficiencyCalibrator::WindowWidth() / 2);
   fEfficiencyCanvas        = new TRootEmbeddedCanvas("EfficiencyCanvas", fRightFrame, TEfficiencyCalibrator::WindowWidth() / 2, 450);
   fStatusBar               = new TGStatusBar(fRightFrame, TEfficiencyCalibrator::WindowWidth() / 2, 50);
   std::array<int, 2> parts = {35, 65};
   fStatusBar->SetParts(parts.data(), parts.size());
   fGraphParameterFrame = new TGGroupFrame(fRightFrame, "Graph Parameters", kHorizontalFrame);
   fGraphParameterFrame->SetLayoutManager(new TGTableLayout(fGraphParameterFrame, 3, 4, false, 5));   // rows, columns, not homogenous cell sizes, 5 = interval between frames in pixel
   fDegreeEntry                 = new TGNumberEntry(fGraphParameterFrame, TEfficiencyCalibrator::Degree(), 2, kDegreeEntry, TGNumberFormat::EStyle::kNESInteger);
   fDegreeLabel                 = new TGLabel(fGraphParameterFrame, "Type of efficiency curve");
   fCalibrationUncertaintyLabel = new TGLabel(fGraphParameterFrame, "Calibration Uncertainty");
   fCalibrationUncertaintyEntry = new TGNumberEntry(fGraphParameterFrame, TEfficiencyCalibrator::CalibrationUncertainty(), 5, kCalibrationUncertaintyEntry, TGNumberFormat::EStyle::kNESRealOne, TGNumberFormat::EAttribute::kNEAPositive);
   fPlotOptionFrame             = new TGGroupFrame(fGraphParameterFrame, "Plot Selection", kHorizontalFrame);
   fPlotOptionFrame->SetLayoutManager(new TGTableLayout(fPlotOptionFrame, 2, 3, false, 5));   // rows, columns, not homogenous cell sizes, 5 = interval between frames in pixel
   fPlotEfficiencyCheck = new TGCheckButton(fPlotOptionFrame, "Efficiency", kPlotEfficiencyCheck);
   fPlotEfficiencyCheck->SetState(kButtonDown);
   fPlotUncorrEfficiencyCheck = new TGCheckButton(fPlotOptionFrame, "Uncorrected efficiency", kPlotUncorrEfficiencyCheck);
   fPlotUncorrEfficiencyCheck->SetState(kButtonUp);
   fPlotPeakAreaCheck = new TGCheckButton(fPlotOptionFrame, "Peak area", kPlotPeakAreaCheck);
   fPlotPeakAreaCheck->SetState(kButtonUp);
   fPlotSummingInCheck = new TGCheckButton(fPlotOptionFrame, "Summing in", kPlotSummingInCheck);
   fPlotSummingInCheck->SetState(kButtonUp);
   fPlotSummingOutCheck = new TGCheckButton(fPlotOptionFrame, "Summing out", kPlotSummingOutCheck);
   fPlotSummingOutCheck->SetState(kButtonUp);
   fPlotOptionFrame->AddFrame(fPlotEfficiencyCheck, new TGTableLayoutHints(0, 1, 0, 1));
   fPlotOptionFrame->AddFrame(fPlotUncorrEfficiencyCheck, new TGTableLayoutHints(1, 3, 0, 1));
   fPlotOptionFrame->AddFrame(fPlotPeakAreaCheck, new TGTableLayoutHints(0, 1, 1, 2));
   fPlotOptionFrame->AddFrame(fPlotSummingInCheck, new TGTableLayoutHints(1, 2, 1, 2));
   fPlotOptionFrame->AddFrame(fPlotSummingOutCheck, new TGTableLayoutHints(2, 3, 1, 2));
   fRecalculateButton = new TGTextButton(fGraphParameterFrame, "Recalculate graphs");
   fGraphParameterFrame->AddFrame(fDegreeLabel, new TGTableLayoutHints(0, 1, 0, 1));
   fGraphParameterFrame->AddFrame(fDegreeEntry, new TGTableLayoutHints(1, 2, 0, 1));
   fGraphParameterFrame->AddFrame(fCalibrationUncertaintyLabel, new TGTableLayoutHints(0, 1, 1, 2));
   fGraphParameterFrame->AddFrame(fCalibrationUncertaintyEntry, new TGTableLayoutHints(1, 2, 1, 2));
   fGraphParameterFrame->AddFrame(fPlotOptionFrame, new TGTableLayoutHints(2, 4, 0, 2));
   fGraphParameterFrame->AddFrame(fRecalculateButton, new TGTableLayoutHints(0, 4, 2, 3));
   fRightFrame->AddFrame(fEfficiencyCanvas, new TGLayoutHints(kLHintsExpandX));
   fRightFrame->AddFrame(fStatusBar, new TGLayoutHints(kLHintsExpandX));
   fRightFrame->AddFrame(fGraphParameterFrame, new TGLayoutHints(kLHintsExpandX));
 
   fFrame->SetLayoutManager(new TGHorizontalLayout(fFrame));
   fFrame->AddFrame(fLeftFrame, new TGLayoutHints(kLHintsLeft | kLHintsTop | kLHintsExpandY, 2, 2, 2, 2));
   fFrame->AddFrame(fRightFrame, new TGLayoutHints(kLHintsRight | kLHintsTop | kLHintsExpandY, 2, 2, 2, 2));
 
   DrawGraph();
   if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kSubroutines) {
      std::cout << "switching from sub directory " << gDirectory->GetName() << " to main directory " << fMainDirectory;
   }
   fMainDirectory->cd();
   if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kSubroutines) {
      std::cout << " = " << fMainDirectory->GetName() << ", sub directory is " << fSubDirectory << " = " << (fSubDirectory == nullptr ? "" : fSubDirectory->GetName()) << std::endl;
   }
   if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kBasicFlow) {
      std::cout << "======================================== done creating tab for data type " << dataType << std::endl;
   }
}
 
TEfficiencyDatatypeTab::~TEfficiencyDatatypeTab()
{
   fMainDirectory->cd();
}
 
void TEfficiencyDatatypeTab::MakeConnections()
{
   fFittingControlGroup->Connect("Clicked(Int_t)", "TEfficiencyDatatypeTab", this, "FittingControl(Int_t)");
   fPlotEfficiencyCheck->Connect("Clicked()", "TEfficiencyDatatypeTab", this, "DrawGraph()");
   fPlotUncorrEfficiencyCheck->Connect("Clicked()", "TEfficiencyDatatypeTab", this, "DrawGraph()");
   fPlotPeakAreaCheck->Connect("Clicked()", "TEfficiencyDatatypeTab", this, "DrawGraph()");
   fPlotSummingInCheck->Connect("Clicked()", "TEfficiencyDatatypeTab", this, "DrawGraph()");
   fPlotSummingOutCheck->Connect("Clicked()", "TEfficiencyDatatypeTab", this, "DrawGraph()");
   fRecalculateButton->Connect("Clicked()", "TEfficiencyDatatypeTab", this, "DrawGraph()");
   fEfficiencyCanvas->GetCanvas()->Connect("ProcessedEvent(Int_t, Int_t, Int_t, TObject*)", "TEfficiencyDatatypeTab", this, "Status(Int_t, Int_t, Int_t, TObject*)");
   for(auto& tab : fEfficiencyTab) {
      tab->MakeConnections();
   }
   fPeakTypeBox->Connect("Selected(Int_t)", "TEfficiencyDatatypeTab", this, "UpdatePeakType()");
}
 
void TEfficiencyDatatypeTab::Disconnect()
{
   fFittingControlGroup->Disconnect("Clicked(Int_t)", this, "FittingControl(Int_t)");
   fPlotEfficiencyCheck->Disconnect("Clicked()", this, "DrawGraph()");
   fPlotUncorrEfficiencyCheck->Disconnect("Clicked()", this, "DrawGraph()");
   fPlotPeakAreaCheck->Disconnect("Clicked()", this, "DrawGraph()");
   fPlotSummingInCheck->Disconnect("Clicked()", this, "DrawGraph()");
   fPlotSummingOutCheck->Disconnect("Clicked()", this, "DrawGraph()");
   fRecalculateButton->Disconnect("Clicked()", this, "DrawGraph()");
   fEfficiencyCanvas->GetCanvas()->Disconnect("ProcessedEvent(Int_t, Int_t, Int_t, TObject*)", this, "Status(Int_t, Int_t, Int_t, TObject*)");
   for(auto& tab : fEfficiencyTab) {
      tab->Disconnect();
   }
}
 
void TEfficiencyDatatypeTab::Status(Int_t, Int_t px, Int_t py, TObject* selected)
{
   fStatusBar->SetText(selected->GetName(), 0);
   fStatusBar->SetText(selected->GetObjectInfo(px, py), 1);
}
 
void TEfficiencyDatatypeTab::DrawGraph()
{
   ReadValues();
   UpdateEfficiencyGraph();
   FitEfficiency();
   if(fLegend == nullptr) {
      fLegend = new TLegend(0.8, 0.95 - static_cast<double>(fNucleus.size()) * 0.05, 0.95, 0.95);
   } else {
      fLegend->Clear();
   }
 
   fEfficiencyCanvas->GetCanvas()->cd();
   bool firstPlot = true;
   int  color     = 1;
   if(fPlotEfficiencyCheck->IsDown()) {
      if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kSubroutines) {
         std::cout << "Drawing efficiency graph " << fEfficiencyGraph << ":" << std::endl;
         fEfficiencyGraph->Print("e");
      }
      for(size_t g = 0; g < fEfficiencyGraph->NumberOfGraphs(); ++g) {
         fEfficiencyGraph->SetColorStyle(g, color++);
      }
      fEfficiencyGraph->DrawCalibration("*", fLegend);
      firstPlot = false;
   }
   if(fPlotUncorrEfficiencyCheck->IsDown()) {
      if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kSubroutines) {
         std::cout << "Drawing uncorrected efficiency graph " << fUncorrEfficiencyGraph << ":" << std::endl;
         fUncorrEfficiencyGraph->Print("e");
      }
      for(size_t g = 0; g < fUncorrEfficiencyGraph->NumberOfGraphs(); ++g) {
         fUncorrEfficiencyGraph->SetColorStyle(g, color++);
      }
      if(firstPlot) {
         fUncorrEfficiencyGraph->DrawCalibration("*", fLegend);
      } else {
         fUncorrEfficiencyGraph->DrawCalibration("same*", fLegend);
      }
      firstPlot = false;
   }
   if(fPlotPeakAreaCheck->IsDown()) {
      if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kSubroutines) {
         fPeakAreaGraph->Print("");
      }
      for(size_t g = 0; g < fPeakAreaGraph->NumberOfGraphs(); ++g) {
         fPeakAreaGraph->SetColorStyle(g, color++);
      }
      if(firstPlot) {
         fPeakAreaGraph->DrawCalibration("*", fLegend);
      } else {
         fPeakAreaGraph->DrawCalibration("same*", fLegend);
      }
      firstPlot = false;
   }
   if(fPlotSummingInCheck->IsDown()) {
      if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kSubroutines) {
         fSummingInGraph->Print("");
      }
      for(size_t g = 0; g < fSummingInGraph->NumberOfGraphs(); ++g) {
         fSummingInGraph->SetColorStyle(g, color++);
      }
      if(firstPlot) {
         fSummingInGraph->DrawCalibration("*", fLegend);
      } else {
         fSummingInGraph->DrawCalibration("same*", fLegend);
      }
      firstPlot = false;
   }
   if(fPlotSummingInCheck->IsDown()) {
      if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kSubroutines) {
         fSummingOutGraph->Print("");
      }
      for(size_t g = 0; g < fSummingOutGraph->NumberOfGraphs(); ++g) {
         fSummingOutGraph->SetColorStyle(g, color++);
      }
      if(firstPlot) {
         fSummingOutGraph->DrawCalibration("*", fLegend);
      } else {
         fSummingOutGraph->DrawCalibration("same*", fLegend);
      }
      firstPlot = false;
   }
   fLegend->Draw();
   fEfficiencyCanvas->GetCanvas()->Modified();
}
 
void TEfficiencyDatatypeTab::UpdateEfficiencyGraph()
{
   if(fMainDirectory == nullptr) {
      std::cout << "No main directory open (" << fMainDirectory << "), sub directory is " << fSubDirectory << ", gDirectory is " << gDirectory->GetName() << std::endl;
      exit(1);
   }
   if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kSubroutines) {
      std::cout << __PRETTY_FUNCTION__ << ": switching from gDirectory " << gDirectory->GetName() << " to main directory " << fMainDirectory << " = " << fMainDirectory->GetName() << std::endl;   // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay)
   }
   fMainDirectory->cd();
   //if(fSubDirectory == nullptr) {
   //   std::cout << "No subdirectory open (" << fSubDirectory << "), main directory is " << fMainDirectory << ", gDirectory is " << gDirectory->GetName() << std::endl;
   //} else {
   //   if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kSubroutines) {
   //      std::cout << __PRETTY_FUNCTION__ << ": switching from gDirectory " << gDirectory->GetName() << " to sub directory " << fSubDirectory;   // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay)
   //   }
   //   fSubDirectory->cd();
   //   if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kSubroutines) {
   //      std::cout << " = " << fSubDirectory->GetName() << ", main directory is " << fMainDirectory << " = " << (fMainDirectory == nullptr ? "" : fMainDirectory->GetName()) << std::endl;
   //   }
   //}
   delete fEfficiencyGraph;
   fEfficiencyGraph = new TCalibrationGraphSet();
   delete fUncorrEfficiencyGraph;
   fUncorrEfficiencyGraph = new TCalibrationGraphSet();
   delete fPeakAreaGraph;
   fPeakAreaGraph = new TCalibrationGraphSet();
   delete fSummingInGraph;
   fSummingInGraph = new TCalibrationGraphSet();
   delete fSummingOutGraph;
   fSummingOutGraph = new TCalibrationGraphSet();
 
   for(auto& tab : fEfficiencyTab) {
      // vector of tuple with transition and 8 doubles:
      // centroid, centroidErr, correctedArea, correctedAreaErr, peak->Area(), peak->AreaErr(), summingIn, summingOut
      auto                peaks            = tab->Peaks();
      TTransition*        transition       = nullptr;
      double              centroid         = 0.;
      double              centroidErr      = 0.;
      double              correctedArea    = 0.;
      double              correctedAreaErr = 0.;
      double              peakArea         = 0.;
      double              peakAreaErr      = 0.;
      double              summingIn        = 0.;
      double              summingOut       = 0.;
      std::vector<double> energy;
      std::vector<double> energyErr;
      std::vector<double> efficiency;
      std::vector<double> efficiencyErr;
      std::vector<double> uncorrEfficiency;
      std::vector<double> uncorrEfficiencyErr;
      std::vector<double> peakAreaVec;
      std::vector<double> peakAreaErrVec;
      std::vector<double> summingInVec;
      std::vector<double> summingOutVec;
      int                 goodPeaks = 0;
      for(auto& peak : peaks) {
         std::tie(transition, centroid, centroidErr, correctedArea, correctedAreaErr, peakArea, peakAreaErr, summingIn, summingOut) = peak;
         if(std::abs(transition->GetEnergy() - centroid) < transition->GetEnergyUncertainty() + centroidErr + TEfficiencyCalibrator::CalibrationUncertainty()) {
            energy.push_back(transition->GetEnergy());
            energyErr.push_back(transition->GetEnergyUncertainty());
            efficiency.push_back(correctedArea / transition->GetIntensity());
            efficiencyErr.push_back(TMath::Sqrt(TMath::Power(correctedAreaErr / transition->GetIntensity(), 2) + TMath::Power(correctedArea / transition->GetIntensity() * transition->GetIntensityUncertainty() / transition->GetIntensity(), 2)));
            uncorrEfficiency.push_back(peakArea / transition->GetIntensity());
            uncorrEfficiencyErr.push_back(TMath::Sqrt(TMath::Power(peakAreaErr / transition->GetIntensity(), 2) + TMath::Power(peakArea / transition->GetIntensity() * transition->GetIntensityUncertainty() / transition->GetIntensity(), 2)));
            peakAreaVec.push_back(peakArea);
            peakAreaErrVec.push_back(peakAreaErr);
            summingInVec.push_back(summingIn);
            summingOutVec.push_back(summingOut);
            ++goodPeaks;
         } else if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kSubroutines) {
            std::cout << "Rejecting centroid " << centroid << " +- " << centroidErr << " with area " << correctedArea << " +- " << correctedAreaErr << " for transition at " << transition->GetEnergy() << " +- " << transition->GetEnergyUncertainty() << " using calibration uncertainty " << TEfficiencyCalibrator::CalibrationUncertainty() << " (" << std::abs(transition->GetEnergy() - centroid) << " >= " << transition->GetEnergyUncertainty() + centroidErr + TEfficiencyCalibrator::CalibrationUncertainty() << ")" << std::endl;
         }
      }
      auto* effGraph = new TGraphErrors(goodPeaks, energy.data(), efficiency.data(), energyErr.data(), efficiencyErr.data());
      fEfficiencyGraph->Add(effGraph, tab->GetName());
      auto* uncorrEffGraph = new TGraphErrors(goodPeaks, energy.data(), uncorrEfficiency.data());
      fUncorrEfficiencyGraph->Add(uncorrEffGraph, Form("%s uncorr. Eff.", tab->GetName()));
      auto* peakAreaGraph = new TGraphErrors(goodPeaks, energy.data(), peakAreaVec.data(), energyErr.data(), peakAreaErrVec.data());
      fPeakAreaGraph->Add(peakAreaGraph, tab->GetName());
      auto* summingInGraph = new TGraphErrors(goodPeaks, energy.data(), summingInVec.data());
      fSummingInGraph->Add(summingInGraph, Form("%s summing in", tab->GetName()));
      auto* summingOutGraph = new TGraphErrors(goodPeaks, energy.data(), summingOutVec.data());
      fSummingOutGraph->Add(summingOutGraph, Form("%s summing out", tab->GetName()));
   }
   fEfficiencyGraph->Scale();
   fUncorrEfficiencyGraph->Scale();
   fEfficiencyGraph->SetAxisTitle("energy [keV];corrected efficiency [a.u.]");
   fUncorrEfficiencyGraph->SetAxisTitle("energy [keV];uncorrected efficiency [a.u.]");
   fPeakAreaGraph->SetAxisTitle("energy [keV];peak areas [a.u.]");
   fSummingInGraph->SetAxisTitle("energy [keV];summing in corrections [a.u.]");
   fSummingOutGraph->SetAxisTitle("energy [keV];summing out corrections [a.u.]");
   fEfficiencyGraph->Write(Form("%s_efficiency", fDataType.c_str()), TObject::kOverwrite);
   fUncorrEfficiencyGraph->Write(Form("%s_uncorr_efficiency", fDataType.c_str()), TObject::kOverwrite);
   fPeakAreaGraph->Write(Form("%s_peak_area", fDataType.c_str()), TObject::kOverwrite);
   fSummingInGraph->Write(Form("%s_summing_in", fDataType.c_str()), TObject::kOverwrite);
   fSummingOutGraph->Write(Form("%s_summing_out", fDataType.c_str()), TObject::kOverwrite);
   fMainDirectory->cd();
}
 
void TEfficiencyDatatypeTab::UpdatePeakType()
{
   if(fPeakTypeBox == nullptr) {
      throw std::runtime_error("Failed to find peak type box, but got a signal it was selected?");
   }
   if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kLoops) {
      std::cout << "peak type box " << fPeakTypeBox << std::flush << ", getting selected " << fPeakTypeBox->GetSelected() << std::endl;
   }
 
   TEfficiencyCalibrator::PeakType(static_cast<TEfficiencyTab::EPeakType>(fPeakTypeBox->GetSelected()));
}
 
int TEfficiencyDatatypeTab::Degree()
{
   if(fDegreeEntry != nullptr) {
      TEfficiencyCalibrator::Degree(static_cast<int>(fDegreeEntry->GetNumber()));
   }
   return TEfficiencyCalibrator::Degree();
}
 
/// The efficiency curves can be:
/// e(E) = ln E + (ln E)/E + (ln E)^2/E + (ln E)^3/E + (ln E)^4/E,
/// or the radware form
/// ln(e(E)) = ((a1 + a2 x + a3 x^2)^-a7 + (a4 + a5 y + a6 y^2)^-a7)^-1/a7
/// with x = ln(E/100), y = ln(E/1000)
/// or a polynomial ln(e(E)) = sum i 0->8 a_i (ln(E))^i (Ryan's & Andrew's PhD theses)
double TEfficiencyDatatypeTab::EfficiencyDebertin(double* x, double* par)   // NOLINT(readability-non-const-parameter)
{
   double eff = par[0] * TMath::Log(x[0]) + par[1] * TMath::Log(x[0]) / x[0] + par[2] * TMath::Power(TMath::Log(x[0]) / x[0], 2) + par[3] * TMath::Power(TMath::Log(x[0]) / x[0], 4) + par[4] * TMath::Power(TMath::Log(x[0]) / x[0], 5);
   return eff;
}
 
double TEfficiencyDatatypeTab::EfficiencyRadware(double* x, double* par)   // NOLINT(readability-non-const-parameter)
{
   double logX   = TMath::Log(x[0] / 100.);
   double logY   = TMath::Log(x[0] / 1000.);
   double logEff = TMath::Power(par[0] + par[1] * logX + par[2] * logX * logX, -par[6]) + TMath::Power(par[3] + par[4] * logY + par[5] * logY * logY, -1. / par[6]);
   return TMath::Exp(logEff);
}
 
double TEfficiencyDatatypeTab::EfficiencyPolynomial(double* x, double* par)   // NOLINT(readability-non-const-parameter)
{
   // first parameter: degree of polynomial
   double logEff = 0;
   for(int i = 0; i < par[0]; ++i) {
      logEff += par[i + 1] * TMath::Power(TMath::Log(x[0]), i);
   }
   return TMath::Exp(logEff);
}
 
void TEfficiencyDatatypeTab::FitEfficiency()
{
   ReadValues();
   delete fEfficiency;
   // degree of fit function: 0 = ln(E)/E form, 1 = radware form, everything else polynomial ln(e(E)) = sum i 0->8 a_i (ln(E))^i (Ryan's & Andrew's PhD theses)
   switch(TEfficiencyCalibrator::Degree()) {
   case 0:
      fEfficiency = new TF1("EfficiencyDebertin", TEfficiencyDatatypeTab::EfficiencyDebertin, 0., 10000., 5);
      break;
   case 1:
      fEfficiency = new TF1("EfficiencyRadware", TEfficiencyDatatypeTab::EfficiencyRadware, 0., 10000., 7);
      break;
   default:
      fEfficiency = new TF1("EfficiencyPolynomial", TEfficiencyDatatypeTab::EfficiencyPolynomial, 0., 10000., TEfficiencyCalibrator::Degree() + 1);
      fEfficiency->FixParameter(0, TEfficiencyCalibrator::Degree());
      break;
   }
 
   fEfficiencyGraph->Fit(fEfficiency);
}
 
void TEfficiencyDatatypeTab::FittingControl(Int_t id)
{
   int currentTabId = fDataTab->GetCurrent();
   if(TEfficiencyCalibrator::VerboseLevel() > EVerbosity::kBasicFlow) { std::cout << __PRETTY_FUNCTION__ << ": id " << id << ", current tab id " << currentTabId << std::endl; }   // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay)
   ReadValues();
   switch(id) {
   case 1:   // re-fit
      fEfficiencyTab[currentTabId]->FindPeaks();
      UpdateEfficiencyGraph();
      break;
   case 2:   // re-fit all
      for(auto& tab : fEfficiencyTab) {
         tab->FindPeaks();
      }
      UpdateEfficiencyGraph();
      break;
   default:
      break;
   }
}
 
void TEfficiencyDatatypeTab::ReadValues()
{
   TEfficiencyCalibrator::Range(fRangeEntry->GetNumber());
   TEfficiencyCalibrator::Threshold(fThresholdEntry->GetNumber());
   TEfficiencyCalibrator::BgParam(fBgParamEntry->GetNumberEntry()->GetIntNumber());
   TEfficiencyCalibrator::Degree(fDegreeEntry->GetNumberEntry()->GetIntNumber());
   TEfficiencyCalibrator::CalibrationUncertainty(fCalibrationUncertaintyEntry->GetNumber());
}
 
//////////////////////////////////////// TEfficiencyCalibrator ////////////////////////////////////////
EVerbosity TEfficiencyCalibrator::fVerboseLevel = EVerbosity::kQuiet;
 
double                    TEfficiencyCalibrator::fRange                  = 20.;
double                    TEfficiencyCalibrator::fThreshold              = 100.;
int                       TEfficiencyCalibrator::fBgParam                = 20;
TEfficiencyTab::EPeakType TEfficiencyCalibrator::fPeakType               = TEfficiencyTab::EPeakType::kRWPeak;
int                       TEfficiencyCalibrator::fDegree                 = 0;
double                    TEfficiencyCalibrator::fCalibrationUncertainty = 1.;
bool                      TEfficiencyCalibrator::fShowRemovedFits        = false;
 
unsigned int TEfficiencyCalibrator::fLineHeight  = 20;
unsigned int TEfficiencyCalibrator::fWindowWidth = 1200;
 
TEfficiencyCalibrator::TEfficiencyCalibrator(int n...)
   : TGMainFrame(nullptr, TEfficiencyCalibrator::WindowWidth() + 10, TEfficiencyCalibrator::WindowWidth() / 2)   // +10 for padding between frames
{
   // remove old file with redirect from fitting
   std::remove("fitOutput.txt");
 
   // we want to store the histograms in a vector of types (suppressed/addback), but we don't know yet if we have all types present
   // so we create a vector of 4 now (unsuppressed singles, suppressed singles, unsuppressed addback, and suppressed addback)
   // and then remove the missing ones at the end
   fHistograms.resize(5);
   fDataType = {"Unsuppressed Singles", "Suppressed Singles", "Unsuppressed Addback", "Suppressed Addback", "Single Crystal"};
 
   va_list args;
   va_start(args, n);   // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay)
   bool allSourcesFound = true;
   for(int i = 0; i < n; ++i) {
      const char* name = va_arg(args, const char*);   // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay)
      fFiles.push_back(new TFile(name));
      if(!fFiles.back()->IsOpen()) {
         std::cerr << "Failed to open " << i + 1 << ". file \"" << name << "\"" << std::endl;
         fFiles.pop_back();
         continue;
      }
      bool goodFile = false;
      if(fVerboseLevel > EVerbosity::kQuiet) {
         std::cout << "Checking file \"" << name << "\":";
      }
      if(fFiles.back()->FindKey("griffinE") != nullptr && fFiles.back()->FindKey("griffinGriffinEProj") != nullptr && fFiles.back()->FindKey("griffinGriffinESum") != nullptr) {
         // unsuppressed singles data
         fHistograms[0].emplace_back(
            static_cast<TH1*>(fFiles.back()->Get("griffinE")),
            static_cast<TH1*>(fFiles.back()->Get("griffinGriffinEProj")),
            static_cast<TH2*>(fFiles.back()->Get("griffinGriffinESum")));
         goodFile = true;
         if(fVerboseLevel > EVerbosity::kQuiet) {
            std::cout << " found unsuppressed singles";
         }
      }
      if(fFiles.back()->FindKey("griffinESupp") != nullptr && fFiles.back()->FindKey("griffinGriffinESuppProj") != nullptr && fFiles.back()->FindKey("griffinGriffinESuppSum") != nullptr) {
         // suppressed singles data
         fHistograms[1].emplace_back(
            static_cast<TH1*>(fFiles.back()->Get("griffinESupp")),
            static_cast<TH1*>(fFiles.back()->Get("griffinGriffinESuppProj")),
            static_cast<TH2*>(fFiles.back()->Get("griffinGriffinESuppSum")));
         goodFile = true;
         if(fVerboseLevel > EVerbosity::kQuiet) {
            std::cout << " found suppressed singles";
         }
      }
      if(fFiles.back()->FindKey("griffinEAddback") != nullptr && fFiles.back()->FindKey("griffinGriffinEAddbackProj") != nullptr && fFiles.back()->FindKey("griffinGriffinEAddbackSum") != nullptr) {
         // unsuppressed addback data
         fHistograms[2].emplace_back(
            static_cast<TH1*>(fFiles.back()->Get("griffinEAddback")),
            static_cast<TH1*>(fFiles.back()->Get("griffinGriffinEAddbackProj")),
            static_cast<TH2*>(fFiles.back()->Get("griffinGriffinEAddbackSum")));
         goodFile = true;
         if(fVerboseLevel > EVerbosity::kQuiet) {
            std::cout << " found unsuppressed addback";
         }
      }
      if(fFiles.back()->FindKey("griffinESuppAddback") != nullptr && fFiles.back()->FindKey("griffinGriffinESuppAddbackProj") != nullptr && fFiles.back()->FindKey("griffinGriffinESuppAddbackSum") != nullptr) {
         // suppressed addback data
         fHistograms[3].emplace_back(
            static_cast<TH1*>(fFiles.back()->Get("griffinESuppAddback")),
            static_cast<TH1*>(fFiles.back()->Get("griffinGriffinESuppAddbackProj")),
            static_cast<TH2*>(fFiles.back()->Get("griffinGriffinESuppAddbackSum")));
         goodFile = true;
         if(fVerboseLevel > EVerbosity::kQuiet) {
            std::cout << " found suppressed addback";
         }
      }
      if(fFiles.back()->FindKey("griffinESingleCrystal") != nullptr && fFiles.back()->FindKey("griffinGriffinESingleCrystalProj") != nullptr && fFiles.back()->FindKey("griffinGriffinESingleCrystalSum") != nullptr) {
         // suppressed addback data
         fHistograms[4].emplace_back(
            static_cast<TH1*>(fFiles.back()->Get("griffinESingleCrystal")),
            static_cast<TH1*>(fFiles.back()->Get("griffinGriffinESingleCrystalProj")),
            static_cast<TH2*>(fFiles.back()->Get("griffinGriffinESingleCrystalSum")));
         goodFile = true;
         if(fVerboseLevel > EVerbosity::kQuiet) {
            std::cout << " found single crystal";
         }
      }
      if(fVerboseLevel > EVerbosity::kQuiet) {
         std::cout << ": got " << fFiles.size() << " files, " << fHistograms[0].size() << " unsuppressed singles, " << fHistograms[1].size() << " suppressed singles, " << fHistograms[2].size() << " unsuppressed addback, " << fHistograms[3].size() << " suppressed addback" << std::endl;
      }
      if(!goodFile) {
         std::cerr << "Failed to find any histogram in " << i + 1 << ". file \"" << name << "\"" << std::endl;
         fFiles.pop_back();
         continue;
      }
      fSources.push_back(nullptr);
      // check if the file name matches any source
      if(std::strstr(name, "22Na") != nullptr) {
         fSources.back() = new TNucleus("22Na");
      } else if(std::strstr(name, "56Co") != nullptr) {
         fSources.back() = new TNucleus("56Co");
      } else if(std::strstr(name, "60Co") != nullptr) {
         fSources.back() = new TNucleus("60Co");
      } else if(std::strstr(name, "133Ba") != nullptr) {
         fSources.back() = new TNucleus("133Ba");
      } else if(std::strstr(name, "152Eu") != nullptr) {
         fSources.back() = new TNucleus("152Eu");
      } else if(std::strstr(name, "241Am") != nullptr) {
         fSources.back() = new TNucleus("241Am");
      } else {
         allSourcesFound = false;
      }
   }
   va_end(args);   // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay)
 
   // now we remove any type from the histogram vector that has no entries
   auto type = fDataType.begin();
   for(auto it = fHistograms.begin(); it != fHistograms.end();) {
      if(it->size() != fFiles.size()) {
         if(fVerboseLevel > EVerbosity::kBasicFlow) {
            std::cout << std::distance(fHistograms.begin(), it) << ": found " << it->size() << " histograms for " << fFiles.size() << " files, removing this type" << std::endl;
         }
         it   = fHistograms.erase(it);
         type = fDataType.erase(type);
      } else {
         ++it;
         ++type;
      }
   }
 
   if(fVerboseLevel > EVerbosity::kQuiet) {
      std::cout << "Read " << fFiles.size() << " files, got";
      for(auto& vec : fHistograms) {
         std::cout << " " << vec.size();
      }
      std::cout << " histograms" << std::endl;
   }
 
   if(fHistograms.empty()) {
      throw std::runtime_error("Unable to find any suitable histograms in the provided file(s)!");
   }
 
   // quick sanity check, should have at least one vector the size of that vector should equal the source size and file size
   if(fSources.size() != fFiles.size() || fHistograms[0].size() != fFiles.size() || fHistograms.size() != fDataType.size()) {
      std::ostringstream error;
      error << "Wrong sizes, from " << fFiles.size() << " file(s) we got " << fSources.size() << " source(s), and " << fHistograms[0].size() << " histogram(s), " << fDataType.size() << " data types, and " << fHistograms.size() << " data type histogram(s)!" << std::endl;
      throw std::runtime_error(error.str());
   }
 
   fOutput = new TFile("TEfficiencyCalibrator.root", "recreate");
   if(!fOutput->IsOpen()) {
      throw std::runtime_error("Unable to open output file \"TEfficiencyCalibrator.root\"!");
   }
 
   // build the first screen
   if(!allSourcesFound) {
      BuildFirstInterface();
      MakeFirstConnections();
   } else {
      SecondWindow();
   }
 
   // Set a name to the main frame
   SetWindowName("EffiencyCalibrator");
 
   // Map all subwindows of main frame
   MapSubwindows();
 
   // Initialize the layout algorithm
   Resize(GetDefaultSize());
 
   // Map main frame
   MapWindow();
 
   fOutput->Close();
}
 
TEfficiencyCalibrator::~TEfficiencyCalibrator()
{
   for(auto& file : fFiles) {
      file->Close();
   }
   if(fOutput != nullptr) { fOutput->Close(); }
}
 
void TEfficiencyCalibrator::DeleteElement(TGFrame* element)
{
   HideFrame(element);
   RemoveFrame(element);
   // delete element;
   // element = nullptr;
}
 
void TEfficiencyCalibrator::BuildFirstInterface()
{
   /// Build initial interface with histogram <-> source assignment
 
   auto* layoutManager = new TGTableLayout(this, fFiles.size() + 1, 2, true, 5);
   if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << "created table layout manager with 2 columns, " << fFiles.size() + 1 << " rows" << std::endl; }
   SetLayoutManager(layoutManager);
 
   // The matrices and source selection boxes
   size_t i = 0;
   for(i = 0; i < fFiles.size(); ++i) {
      fSourceLabel.push_back(new TGLabel(this, fFiles[i]->GetName()));
      if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "Text height " << fSourceLabel.back()->GetFont()->TextHeight() << std::endl; }
      fSourceBox.push_back(new TGComboBox(this, "Select source", kSourceBox + fSourceBox.size()));
      if(fSources[i] != nullptr) {
         fSourceBox.back()->AddEntry(fSources[i]->GetName(), i);
         fSourceBox.back()->Select(0);
         if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "Found source, setting entry to " << fSources[i]->GetName() << " and selecting " << fSourceBox.back()->GetSelected() << std::endl; }
      } else {
         fSourceBox.back()->AddEntry("22Na", k22Na);
         fSourceBox.back()->AddEntry("56Co", k56Co);
         fSourceBox.back()->AddEntry("60Co", k60Co);
         fSourceBox.back()->AddEntry("133Ba", k133Ba);
         fSourceBox.back()->AddEntry("152Eu", k152Eu);
         fSourceBox.back()->AddEntry("241Am", k241Am);
         if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "Didn't find source, created source box with " << fSourceBox.back()->GetNumberOfEntries() << std::endl; }
      }
      fSourceBox.back()->SetMinHeight(200);
      fSourceBox.back()->Resize(100, LineHeight());
      if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "Attaching " << i << ". label to 0, 1, " << i << ", " << i + 1 << ", and box to 1, 2, " << i << ", " << i + 1 << std::endl; }
      AddFrame(fSourceLabel.back(), new TGTableLayoutHints(0, 1, i, i + 1, kLHintsRight | kLHintsCenterY));
      AddFrame(fSourceBox.back(), new TGTableLayoutHints(1, 2, i, i + 1, kLHintsLeft | kLHintsCenterY));
   }
 
   // The buttons
   if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << "Attaching start button to 0, 2, " << i << ", " << i + 1 << std::endl; }
   fStartButton = new TGTextButton(this, "Accept && Continue", kStartButton);
   if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << "start button " << fStartButton << std::endl; }
   AddFrame(fStartButton, new TGTableLayoutHints(0, 2, i, i + 1, kLHintsCenterX | kLHintsCenterY));
   Layout();
}
 
void TEfficiencyCalibrator::MakeFirstConnections()
{
   /// Create connections for the file <-> source assignment interface
 
   // Connect the selection of the source
   for(auto* box : fSourceBox) {
      box->Connect("Selected(Int_t, Int_t)", "TEfficiencyCalibrator", this, "SetSource(Int_t, Int_t)");
   }
 
   // Connect the clicking of buttons
   fStartButton->Connect("Clicked()", "TEfficiencyCalibrator", this, "Start()");
}
 
void TEfficiencyCalibrator::DisconnectFirst()
{
   /// Disconnect all signals from the file <-> source assignment interface
   for(auto* box : fSourceBox) {
      box->Disconnect("Selected(Int_t, Int_t)", this, "SetSource(Int_t, Int_t)");
   }
   fStartButton->Disconnect("Clicked()", this, "Start()");
}
 
void TEfficiencyCalibrator::DeleteFirst()
{
 
   fSourceBox.clear();
   DeleteElement(fStartButton);
   if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "Deleted start button " << fStartButton << std::endl; }
}
 
void TEfficiencyCalibrator::SetSource(Int_t windowId, Int_t entryId)
{
   int index = windowId - kSourceBox;
   if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << __PRETTY_FUNCTION__ << ": windowId " << windowId << ", entryId " << entryId << " => " << index << std::endl; }   // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay)
   TNucleus* nucleus = fSources[index];
   delete nucleus;
   nucleus         = new TNucleus(fSourceBox[index]->GetListBox()->GetEntry(entryId)->GetTitle());
   fSources[index] = nucleus;
}
 
void TEfficiencyCalibrator::Start()
{
   if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << __PRETTY_FUNCTION__ << ": fEmitter " << fEmitter << ", fStartButton " << fStartButton << std::endl; }   // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay)
   if(fEmitter == nullptr) {                                                                                                                                         // we only want to do this once at the beginning (after fEmitter was initialized to nullptr)
      fEmitter = fStartButton;
      TTimer::SingleShot(100, "TEfficiencyCalibrator", this, "HandleTimer()");
   }
}
 
void TEfficiencyCalibrator::HandleTimer()
{
   if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << __PRETTY_FUNCTION__ << ": fEmitter " << fEmitter << ", fStartButton " << fStartButton << std::endl; }   // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay)
   if(fEmitter == fStartButton) {
      // disconnect signals of first screen and remove all elements
      DisconnectFirst();
      RemoveAll();
      DeleteFirst();
 
      SecondWindow();
   }
}
 
void TEfficiencyCalibrator::SecondWindow()
{
   if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << __PRETTY_FUNCTION__ << std::endl; }   // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay)
   // check that all sources have been set
   for(size_t i = 0; i < fSources.size(); ++i) {
      if(fSources[i] == nullptr) {
         std::cerr << "Source " << i << " not set!" << std::endl;
         return;
      }
      if(fVerboseLevel > EVerbosity::kQuiet) {
         std::cout << i << " - source " << fSources[i] << " = " << fSources[i]->GetName() << std::endl;
      }
   }
   // now check that we don't have the same source twice (which wouldn't make sense)
   for(size_t i = 0; i < fSources.size(); ++i) {
      for(size_t j = i + 1; j < fSources.size(); ++j) {
         if(*(fSources[i]) == *(fSources[j])) {
            std::cerr << "Duplicate sources: " << i << " - " << fSources[i]->GetName() << " and " << j << " - " << fSources[j]->GetName() << std::endl;
            return;
         }
      }
   }
 
   SetLayoutManager(new TGHorizontalLayout(this));
 
   // create intermediate progress bar
   fProgressBar      = new TGHProgressBar(this, TGProgressBar::kFancy, WindowWidth() / 2);
   int nofHistograms = 0;
   for(auto& type : fHistograms) {   // loop over all data types and add the number of source to the counter
      nofHistograms += type.size();
   }
   nofHistograms *= 3;   // there are three histograms for each data type/source combo
   fProgressBar->SetRange(0., static_cast<Float_t>(nofHistograms));
   fProgressBar->Percent(true);
   if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "Set range of progress bar to 0. - " << fProgressBar->GetMax() << " = " << fFiles.size() << std::endl; }
   AddFrame(fProgressBar, new TGLayoutHints(kLHintsCenterX | kLHintsCenterY | kLHintsExpandX | kLHintsExpandY, 0, 0, 0, 0));
 
   // Map all subwindows of main frame
   MapSubwindows();
 
   // Initialize the layout algorithm
   Resize(TGDimension(WindowWidth(), LineHeight()));
 
   // Map main frame
   MapWindow();
 
   // create second screen and its connections
   if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "Starting to build second interface:" << std::endl; }
   BuildSecondInterface();
   MakeSecondConnections();
 
   // remove progress bar
   DeleteElement(fProgressBar);
 
   // Map all subwindows of main frame
   MapSubwindows();
 
   // Initialize the layout algorithm
   Resize(TGDimension(WindowWidth(), WindowWidth() / 2));
 
   // Map main frame
   MapWindow();
   if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << __PRETTY_FUNCTION__ << " done" << std::endl; }   // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay)
}
 
void TEfficiencyCalibrator::BuildSecondInterface()
{
   SetLayoutManager(new TGHorizontalLayout(this));
 
   // left frame with tabs for each source
   fDatatypeTab = new TGTab(this, WindowWidth(), WindowWidth() / 2);
   fEfficiencyDatatypeTab.resize(fHistograms.size());
   fEfficiencyGraph.resize(fHistograms.size());
   if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << __PRETTY_FUNCTION__ << " creating " << fHistograms.size() << " tabs" << std::endl; }   // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay)
   for(size_t i = 0; i < fHistograms.size(); ++i) {
      if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << i << ": Creating efficiency source tab using " << fHistograms[i].size() << " histograms, and " << fSources.size() << " sources, " << fRange << ", " << fThreshold << ", " << fProgressBar << std::endl; }
      auto* frame = fDatatypeTab->AddTab(fDataType[i].c_str());
      std::replace(fDataType[i].begin(), fDataType[i].end(), ' ', '_');
      fEfficiencyDatatypeTab[i] = new TEfficiencyDatatypeTab(this, fSources, fHistograms[i], frame, fDataType[i], fProgressBar);
      // fEfficiencyDatatypeTab[i]->UpdateEfficiencyGraph();
      fEfficiencyGraph[i] = fEfficiencyDatatypeTab[i]->EfficiencyGraph();
   }
   AddFrame(fDatatypeTab, new TGLayoutHints(kLHintsTop | kLHintsCenterX | kLHintsExpandX, 0, 0, 0, 0));
 
   if(fVerboseLevel > EVerbosity::kSubroutines) { std::cout << "Second interface done" << std::endl; }
}
 
void TEfficiencyCalibrator::MakeSecondConnections()
{
   if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << __PRETTY_FUNCTION__ << std::endl; }   // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay)
   // we don't need to connect the range, threshold, and degree number entries, those are automatically read when we start the calibration
   for(auto* sourceTab : fEfficiencyDatatypeTab) {
      sourceTab->MakeConnections();
   }
}
 
void TEfficiencyCalibrator::DisconnectSecond()
{
   if(fVerboseLevel > EVerbosity::kBasicFlow) { std::cout << __PRETTY_FUNCTION__ << std::endl; }   // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay)
   for(auto* sourceTab : fEfficiencyDatatypeTab) {
      sourceTab->Disconnect();
   }
}
 
#endif