From ecec98760061a50b061fe120caa4da89e6de84fc Mon Sep 17 00:00:00 2001 From: Alvaro Ezquerro Date: Fri, 15 Mar 2024 11:33:15 +0100 Subject: [PATCH 01/12] Documentation fix at header --- .../analysis/inc/TRestDataSetGainMap.h | 114 +++++++++++++----- 1 file changed, 81 insertions(+), 33 deletions(-) diff --git a/source/framework/analysis/inc/TRestDataSetGainMap.h b/source/framework/analysis/inc/TRestDataSetGainMap.h index 50000a01b..adbe5f655 100644 --- a/source/framework/analysis/inc/TRestDataSetGainMap.h +++ b/source/framework/analysis/inc/TRestDataSetGainMap.h @@ -42,13 +42,24 @@ class TRestDataSetGainMap : public TRestMetadata { class Module; private: - std::string fObservable = ""; //"rawAna_ThresholdIntegral"; //< - std::string fSpatialObservableX = ""; //"hitsAna_xMean"; //< - std::string fSpatialObservableY = ""; //"hitsAna_yMean"; //< + /// Observable that will be used to calculate the gain map + std::string fObservable = ""; //< + + /// Observable that will be used to segmentize the gain map in the x direction + std::string fSpatialObservableX = ""; //< + + /// Observable that will be used to segmentize the gain map in the y direction + std::string fSpatialObservableY = ""; //< + + /// List of modules + std::vector fModulesCal = {}; //< + + /// Name of the file that contains the calibration data + std::string fCalibFileName = ""; //< + + /// Name of the file where the gain map was (or will be) exported + std::string fOutputFileName = ""; //< - std::vector fModulesCal = {}; - std::string fCalibFileName = ""; - std::string fOutputFileName = ""; void Initialize() override; void InitFromConfigFile() override; @@ -90,46 +101,83 @@ class TRestDataSetGainMap : public TRestMetadata { TRestDataSetGainMap& operator=(TRestDataSetGainMap& src); - public: void PrintMetadata() override; void GenerateGainMap(); - void CalibrateDataSet(const std::string& dataSetFileName, std::string outputFileName = ""); + void CalibrateDataSet(const std::string& dataSetFileName, std::string outputFileName = "", + std::vector excludeColumns = {}); TRestDataSetGainMap(); TRestDataSetGainMap(const char* configFilename, std::string name = ""); ~TRestDataSetGainMap(); - ClassDefOverride(TRestDataSetGainMap, 1); + ClassDefOverride(TRestDataSetGainMap, 2); class Module { private: - const TRestDataSetGainMap* p = nullptr; // fEnergyPeaks = {}; - std::vector fRangePeaks = {}; //{TVector2(230000, 650000), TVector2(40000, 230000)}; - TVector2 fCalibRange = TVector2(0, 0); //< // Calibration range - Int_t fNBins = 100; //< // Number of bins for the spectrum histograms - std::string fDefinitionCut = ""; //"TREXsides_tagId == 2"; //< - - Int_t fNumberOfSegmentsX = 1; //< - Int_t fNumberOfSegmentsY = 1; //< - TVector2 fReadoutRange = TVector2(-1, 246.24); //< // Readout dimensions - std::set fSplitX = {}; //< - std::set fSplitY = {}; //< - - std::string fDataSetFileName = ""; //< // File name for the calibration dataset - - std::vector> fSlope = {}; //< + /// Pointer to the parent class + const TRestDataSetGainMap* p = nullptr; // fEnergyPeaks = {}; //< + + /// Range of the peaks to be used for the calibration. If empty it will be automatically calculated. + std::vector fRangePeaks = {}; //< + + /// Calibration range. If fCalibRange.X()>=fCalibRange.Y() the range will be automatically calculated. + TVector2 fCalibRange = TVector2(0, 0); //< + + /// Number of bins for the spectrum histograms. + Int_t fNBins = 100; //< + + /// Cut that defines which events are from this module. + std::string fDefinitionCut = ""; //< + + /// Number of segments in the x direction. + Int_t fNumberOfSegmentsX = 1; //< + + /// Number of segments in the y direction. + Int_t fNumberOfSegmentsY = 1; //< + + /// Readout dimensions + TVector2 fReadoutRange = TVector2(-1, 246.24); //< + + /// Split points in the x direction. + std::set fSplitX = {}; //< + + /// Split points in the y direction. + std::set fSplitY = {}; //< + + /// Name of the file that contains the calibration data. If empty, it will use its parent + /// TRestDataSetGainMap::fCalibFileName. + std::string fDataSetFileName = ""; //< + + /// Array containing the slope of the linear fit for each segment. + std::vector> fSlope = {}; //< + + /// Array containing the intercept of the linear fit for each segment. std::vector> fIntercept = {}; //< - bool fZeroPoint = false; //< Zero point will be automatically added if there are less than 2 peaks - bool fAutoRangePeaks = true; //< Automatic range peaks - std::vector> fSegSpectra = {}; - std::vector> fSegLinearFit = {}; + /// Add zero point to the calibration linear fit. Zero point will be automatically added if there are + /// less than 2 points. + bool fZeroPoint = false; //< + + /// Automatic range for the peaks fitting. See GenerateGainMap() for more information of the logic. + bool fAutoRangePeaks = true; //< + + /// Array containing the observable spectrum for each segment. + std::vector> fSegSpectra = {}; //< + + /// Array containing the calibration linear fit for each segment. + std::vector> fSegLinearFit = {}; //< public: void AddPeak(const double& energyPeak, const TVector2& rangePeak = TVector2(0, 0)) { From 6ea5197f2e186c009ff08840bf398c7061fd5bb6 Mon Sep 17 00:00:00 2001 From: Alvaro Ezquerro Date: Fri, 15 Mar 2024 11:35:33 +0100 Subject: [PATCH 02/12] EnableMultiThreading when handling TRestDataSet --- source/framework/analysis/src/TRestDataSetGainMap.cxx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/source/framework/analysis/src/TRestDataSetGainMap.cxx b/source/framework/analysis/src/TRestDataSetGainMap.cxx index 53a5fa0a0..c9ed305d7 100644 --- a/source/framework/analysis/src/TRestDataSetGainMap.cxx +++ b/source/framework/analysis/src/TRestDataSetGainMap.cxx @@ -218,6 +218,7 @@ void TRestDataSetGainMap::CalibrateDataSet(const std::string& dataSetFileName, s } TRestDataSet dataSet; + dataSet.EnableMultiThreading(true); dataSet.Import(dataSetFileName); auto dataFrame = dataSet.GetDataFrame(); @@ -628,6 +629,7 @@ void TRestDataSetGainMap::Module::GenerateGainMap() { } if (!TRestTools::isDataSet(dsFileName)) RESTWarning << dsFileName << " is not a dataset." << p->RESTendl; TRestDataSet dataSet; + dataSet.EnableMultiThreading(true); dataSet.Import(dsFileName); fDataSetFileName = dsFileName; From 76dd1adfcb34c8ea5c01e98d91c93037bfb5c5f6 Mon Sep 17 00:00:00 2001 From: Alvaro Ezquerro Date: Fri, 15 Mar 2024 11:55:21 +0100 Subject: [PATCH 03/12] Add a TRestCut --- source/framework/analysis/inc/TRestDataSetGainMap.h | 5 +++++ source/framework/analysis/src/TRestDataSetGainMap.cxx | 6 ++++++ 2 files changed, 11 insertions(+) diff --git a/source/framework/analysis/inc/TRestDataSetGainMap.h b/source/framework/analysis/inc/TRestDataSetGainMap.h index adbe5f655..035cdf294 100644 --- a/source/framework/analysis/inc/TRestDataSetGainMap.h +++ b/source/framework/analysis/inc/TRestDataSetGainMap.h @@ -33,6 +33,7 @@ #include #include +#include "TRestCut.h" #include "TRestDataSet.h" #include "TRestMetadata.h" @@ -60,6 +61,8 @@ class TRestDataSetGainMap : public TRestMetadata { /// Name of the file where the gain map was (or will be) exported std::string fOutputFileName = ""; //< + /// Cut to be applied to the calibration data + TRestCut* fCut = nullptr; //< void Initialize() override; void InitFromConfigFile() override; @@ -80,6 +83,7 @@ class TRestDataSetGainMap : public TRestMetadata { std::string GetObservable() const { return fObservable; } std::string GetSpatialObservableX() const { return fSpatialObservableX; } std::string GetSpatialObservableY() const { return fSpatialObservableY; } + TRestCut* GetCut() const { return fCut; } Module* GetModule(const int planeID, const int moduleID); double GetSlopeParameter(const int planeID, const int moduleID, const double x, const double y); @@ -95,6 +99,7 @@ class TRestDataSetGainMap : public TRestMetadata { void SetSpatialObservableY(const std::string& spatialObservableY) { fSpatialObservableY = spatialObservableY; } + void SetCut(TRestCut* cut) { fCut = cut; } void Import(const std::string& fileName); void Export(const std::string& fileName = ""); diff --git a/source/framework/analysis/src/TRestDataSetGainMap.cxx b/source/framework/analysis/src/TRestDataSetGainMap.cxx index c9ed305d7..29e8ea69c 100644 --- a/source/framework/analysis/src/TRestDataSetGainMap.cxx +++ b/source/framework/analysis/src/TRestDataSetGainMap.cxx @@ -186,6 +186,8 @@ void TRestDataSetGainMap::InitFromConfigFile() { moduleDefinition = GetNextElement(moduleDefinition); } + fCut = (TRestCut*)InstantiateChildMetadata("TRestCut"); + if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) PrintMetadata(); } @@ -347,6 +349,7 @@ TRestDataSetGainMap& TRestDataSetGainMap::operator=(TRestDataSetGainMap& src) { fObservable = src.GetObservable(); fSpatialObservableX = src.GetSpatialObservableX(); fSpatialObservableY = src.GetSpatialObservableY(); + fCut = src.GetCut(); fModulesCal.clear(); for (auto pID : src.GetPlaneIDs()) for (auto mID : src.GetModuleIDs(pID)) fModulesCal.push_back(*src.GetModule(pID, mID)); @@ -432,6 +435,7 @@ void TRestDataSetGainMap::Export(const std::string& fileName) { void TRestDataSetGainMap::PrintMetadata() { TRestMetadata::PrintMetadata(); RESTMetadata << " Calibration dataset: " << fCalibFileName << RESTendl; + if (fCut) fCut->Print(); RESTMetadata << " Output file: " << fOutputFileName << RESTendl; RESTMetadata << " Number of planes: " << GetNumberOfPlanes() << RESTendl; RESTMetadata << " Number of modules: " << GetNumberOfModules() << RESTendl; @@ -633,6 +637,8 @@ void TRestDataSetGainMap::Module::GenerateGainMap() { dataSet.Import(dsFileName); fDataSetFileName = dsFileName; + dataSet.SetDataFrame(dataSet.MakeCut(p->GetCut())); + SetSplits(); if (fSplitX.empty()) SetSplitX(); From 377d72777db47d9ed252e4fd7aaa1bdeb8c2f55f Mon Sep 17 00:00:00 2001 From: Alvaro Ezquerro Date: Fri, 15 Mar 2024 12:02:16 +0100 Subject: [PATCH 04/12] Fix bugs in Draw methods --- .../analysis/inc/TRestDataSetGainMap.h | 2 +- .../analysis/src/TRestDataSetGainMap.cxx | 50 ++++++++++++++----- 2 files changed, 38 insertions(+), 14 deletions(-) diff --git a/source/framework/analysis/inc/TRestDataSetGainMap.h b/source/framework/analysis/inc/TRestDataSetGainMap.h index 035cdf294..08df49944 100644 --- a/source/framework/analysis/inc/TRestDataSetGainMap.h +++ b/source/framework/analysis/inc/TRestDataSetGainMap.h @@ -217,7 +217,7 @@ class TRestDataSetGainMap : public TRestMetadata { TCanvas* c = nullptr); void DrawFullSpectrum(); - void DrawLinearFit(); + void DrawLinearFit(TCanvas* c = nullptr); void DrawLinearFit(const TVector2& position, TCanvas* c = nullptr); void DrawLinearFit(const int index_x, const int index_y, TCanvas* c = nullptr); diff --git a/source/framework/analysis/src/TRestDataSetGainMap.cxx b/source/framework/analysis/src/TRestDataSetGainMap.cxx index 29e8ea69c..9469ce1ca 100644 --- a/source/framework/analysis/src/TRestDataSetGainMap.cxx +++ b/source/framework/analysis/src/TRestDataSetGainMap.cxx @@ -1103,8 +1103,9 @@ void TRestDataSetGainMap::Module::DrawSpectrum(const bool drawFits, const int co for (size_t i = 0; i < fSegSpectra.size(); i++) { for (size_t j = 0; j < fSegSpectra[i].size(); j++) { - c->cd(i + 1 + fSegSpectra[i].size() * j); - DrawSpectrum(i, fSegSpectra[i].size() - 1 - j, drawFits, color, c); + int pad = fSegSpectra.size() * (fSegSpectra[i].size() - 1) + 1 + i - fSegSpectra.size() * j; + c->cd(pad); + DrawSpectrum(i, j, drawFits, color, c); } } } @@ -1160,15 +1161,31 @@ void TRestDataSetGainMap::Module::DrawLinearFit(const int index_x, const int ind fSegLinearFit[index_x][index_y]->Draw("AL*"); } -void TRestDataSetGainMap::Module::DrawLinearFit() { - std::string t = "linearFits_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId); - TCanvas* myCanvas = new TCanvas(t.c_str(), t.c_str()); - myCanvas->Divide(fSegLinearFit.size(), fSegLinearFit.at(0).size()); +void TRestDataSetGainMap::Module::DrawLinearFit(TCanvas* c) { + if (fSegLinearFit.size() == 0) { + RESTError << "Spectra matrix is empty." << p->RESTendl; + return; + } + if (!c) { + std::string t = "linearFits_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId); + c = new TCanvas(t.c_str(), t.c_str()); + } + + size_t nPads = 0; + for (const auto& object : *c->GetListOfPrimitives()) + if (object->InheritsFrom(TVirtualPad::Class())) ++nPads; + if (nPads != 0 && nPads != fSegLinearFit.size() * fSegLinearFit.at(0).size()) { + RESTError << "Canvas " << c->GetName() << " has " << nPads << " pads, but " + << fSegLinearFit.size() * fSegLinearFit.at(0).size() << " are needed." << p->RESTendl; + return; + } else if (nPads == 0) + c->Divide(fSegLinearFit.size(), fSegLinearFit.at(0).size()); + for (size_t i = 0; i < fSegLinearFit.size(); i++) { for (size_t j = 0; j < fSegLinearFit[i].size(); j++) { - myCanvas->cd(i + 1 + fSegLinearFit[i].size() * j); - // fSegLinearFit[i][j]->Draw("AL*"); - DrawLinearFit(i, fSegSpectra[i].size() - 1 - j, myCanvas); + int pad = fSegLinearFit.size() * (fSegLinearFit[i].size() - 1) + 1 + i - fSegLinearFit.size() * j; + c->cd(pad); + DrawLinearFit(i, j, c); } } } @@ -1179,6 +1196,11 @@ void TRestDataSetGainMap::Module::DrawGainMap(const int peakNumber) { << fEnergyPeaks.size() - 1 << " )" << p->RESTendl; return; } + if (fSegLinearFit.size() == 0) { + RESTError << "Linear fit matrix is empty." << p->RESTendl; + return; + } + double peakEnergy = fEnergyPeaks[peakNumber]; std::string title = "Gain map for energy " + DoubleToString(peakEnergy, "%g") + ";" + GetSpatialObservableX() + ";" + GetSpatialObservableY(); // + " keV"; @@ -1186,8 +1208,8 @@ void TRestDataSetGainMap::Module::DrawGainMap(const int peakNumber) { std::to_string(fModuleId); TCanvas* gainMap = new TCanvas(t.c_str(), t.c_str()); gainMap->cd(); - TH2F* hGainMap = new TH2F(("h" + t).c_str(), title.c_str(), fNumberOfSegmentsY, fReadoutRange.X(), - fReadoutRange.Y(), fNumberOfSegmentsX, fReadoutRange.X(), fReadoutRange.Y()); + TH2F* hGainMap = new TH2F(("h" + t).c_str(), title.c_str(), fNumberOfSegmentsX, fReadoutRange.X(), + fReadoutRange.Y(), fNumberOfSegmentsY, fReadoutRange.X(), fReadoutRange.Y()); const double peakPosRef = fSegLinearFit[(fNumberOfSegmentsX - 1) / 2][(fNumberOfSegmentsY - 1) / 2]->GetPointX(peakNumber); @@ -1200,8 +1222,10 @@ void TRestDataSetGainMap::Module::DrawGainMap(const int peakNumber) { auto xUpper = *std::next(itX); auto yLower = *itY; auto yUpper = *std::next(itY); - hGainMap->Fill((xUpper + xLower) / 2., (yUpper + yLower) / 2., - fSegLinearFit[i][j]->GetPointX(peakNumber) / peakPosRef); + float xMean = (xUpper + xLower) / 2.; + float yMean = (yUpper + yLower) / 2.; + auto [index_x, index_y] = GetIndexMatrix(xMean, yMean); + hGainMap->Fill(xMean, yMean, fSegLinearFit[index_x][index_y]->GetPointX(peakNumber) / peakPosRef); itY++; } itX++; From 85d256370e5b7b0f74add4f9fdcda3a82626c564 Mon Sep 17 00:00:00 2001 From: Alvaro Ezquerro Date: Sat, 16 Mar 2024 21:18:09 +0100 Subject: [PATCH 05/12] Adding full module spectrum (new attributes and methods) and use it as reference for DrawGainMap --- .../analysis/inc/TRestDataSetGainMap.h | 21 +- .../analysis/src/TRestDataSetGainMap.cxx | 216 +++++++++++++++--- 2 files changed, 207 insertions(+), 30 deletions(-) diff --git a/source/framework/analysis/inc/TRestDataSetGainMap.h b/source/framework/analysis/inc/TRestDataSetGainMap.h index 08df49944..27a426137 100644 --- a/source/framework/analysis/inc/TRestDataSetGainMap.h +++ b/source/framework/analysis/inc/TRestDataSetGainMap.h @@ -88,6 +88,8 @@ class TRestDataSetGainMap : public TRestMetadata { Module* GetModule(const int planeID, const int moduleID); double GetSlopeParameter(const int planeID, const int moduleID, const double x, const double y); double GetInterceptParameter(const int planeID, const int moduleID, const double x, const double y); + double GetSlopeParameterFullSpc(const int planeID, const int moduleID); + double GetInterceptParameterFullSpc(const int planeID, const int moduleID); void SetCalibrationFileName(const std::string& fileName) { fCalibFileName = fileName; } void SetOutputFileName(const std::string& fileName) { fOutputFileName = fileName; } @@ -168,6 +170,12 @@ class TRestDataSetGainMap : public TRestMetadata { /// Array containing the slope of the linear fit for each segment. std::vector> fSlope = {}; //< + /// Slope of the calibration linear fit of whole module + double fFullSlope = 0; //< + + /// Intercept of the calibration linear fit of whole module + double fFullIntercept = 0; //< + /// Array containing the intercept of the linear fit for each segment. std::vector> fIntercept = {}; //< @@ -181,9 +189,15 @@ class TRestDataSetGainMap : public TRestMetadata { /// Array containing the observable spectrum for each segment. std::vector> fSegSpectra = {}; //< + /// Spectrum of the observable for the whole module. + TH1F* fFullSpectrum = nullptr; //< + /// Array containing the calibration linear fit for each segment. std::vector> fSegLinearFit = {}; //< + /// Calibration linear fit for the whole module. + TGraph* fFullLinearFit = nullptr; //< + public: void AddPeak(const double& energyPeak, const TVector2& rangePeak = TVector2(0, 0)) { fEnergyPeaks.push_back(energyPeak); @@ -195,6 +209,8 @@ class TRestDataSetGainMap : public TRestMetadata { std::pair GetIndexMatrix(const double x, const double y) const; double GetSlope(const double x, const double y) const; double GetIntercept(const double x, const double y) const; + double GetSlopeFullSpc() const { return fFullSlope; }; + double GetInterceptFullSpc() const { return fFullIntercept; }; Int_t GetPlaneId() const { return fPlaneId; } Int_t GetModuleId() const { return fModuleId; } @@ -215,7 +231,7 @@ class TRestDataSetGainMap : public TRestMetadata { TCanvas* c = nullptr); void DrawSpectrum(const int index_x, const int index_y, bool drawFits = true, int color = -1, TCanvas* c = nullptr); - void DrawFullSpectrum(); + void DrawFullSpectrum(const bool drawFits = true, const int color = -1, TCanvas* c = nullptr); void DrawLinearFit(TCanvas* c = nullptr); void DrawLinearFit(const TVector2& position, TCanvas* c = nullptr); @@ -225,7 +241,10 @@ class TRestDataSetGainMap : public TRestMetadata { void Refit(const TVector2& position, const double energy, const TVector2& range); void Refit(const size_t x, const size_t y, const size_t peakNumber, const TVector2& range); + void RefitFullSpc(const double energy, const TVector2& range); + void RefitFullSpc(const size_t peakNumber, const TVector2& range); void UpdateCalibrationFits(const size_t x, const size_t y); + void UpdateCalibrationFitsFullSpc(); void SetPlaneId(const Int_t& planeId) { fPlaneId = planeId; } void SetModuleId(const Int_t& moduleId) { fModuleId = moduleId; } diff --git a/source/framework/analysis/src/TRestDataSetGainMap.cxx b/source/framework/analysis/src/TRestDataSetGainMap.cxx index 9469ce1ca..2d8883646 100644 --- a/source/framework/analysis/src/TRestDataSetGainMap.cxx +++ b/source/framework/analysis/src/TRestDataSetGainMap.cxx @@ -300,6 +300,16 @@ double TRestDataSetGainMap::GetSlopeParameter(const int planeID, const int modul return moduleCal->GetSlope(x, y); } +///////////////////////////////////////////// +/// \brief Function to get the slope parameter of the whole module with +/// planeID and moduleID. +/// +double TRestDataSetGainMap::GetSlopeParameterFullSpc(const int planeID, const int moduleID) { + Module* moduleCal = GetModule(planeID, moduleID); + if (moduleCal == nullptr) return 0; // return numeric_limits::quiet_NaN() + return moduleCal->GetSlopeFullSpc(); +} + ///////////////////////////////////////////// /// \brief Function to get the intercept parameter of the module with /// planeID and moduleID at physical position (x,y) @@ -312,6 +322,16 @@ double TRestDataSetGainMap::GetInterceptParameter(const int planeID, const int m return moduleCal->GetIntercept(x, y); } +///////////////////////////////////////////// +/// \brief Function to get the intercept parameter of the whole module with +/// planeID and moduleID. +/// +double TRestDataSetGainMap::GetInterceptParameterFullSpc(const int planeID, const int moduleID) { + Module* moduleCal = GetModule(planeID, moduleID); + if (moduleCal == nullptr) return 0; // return numeric_limits::quiet_NaN() + return moduleCal->GetInterceptFullSpc(); +} + ///////////////////////////////////////////// /// \brief Function to get a list (set) of the plane IDs /// @@ -671,12 +691,22 @@ void TRestDataSetGainMap::Module::GenerateGainMap() { << p->RESTendl; } + // --- Definition of histogram whole module --- + std::string hModuleName = "hSpc_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId); + TH1F* hModule = new TH1F(hModuleName.c_str(), "", fNBins, fCalibRange.X(), fCalibRange.Y()); + + // build the spectrum for the whole module + std::string cut = fDefinitionCut; + if (cut.empty()) cut = "1"; + auto histoMod = dataSet.GetDataFrame().Filter(cut).Histo1D( + {"tempMod", "", fNBins, fCalibRange.X(), fCalibRange.Y()}, GetObservable()); + hModule = (TH1F*)histoMod->Clone(hModuleName.c_str()); + //--- Definition of histogram matrix --- std::vector> h(fNumberOfSegmentsX, std::vector(fNumberOfSegmentsY, nullptr)); for (size_t i = 0; i < h.size(); i++) { for (size_t j = 0; j < h.at(0).size(); j++) { - std::string name = "hSpc_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId) + "_" + - std::to_string(i) + "_" + std::to_string(j); + std::string name = hModuleName + std::to_string(i) + "_" + std::to_string(j); h[i][j] = new TH1F(name.c_str(), "", fNBins, fCalibRange.X(), fCalibRange.Y()); // h[column][row] equivalent to h[x][y] } @@ -723,12 +753,23 @@ void TRestDataSetGainMap::Module::GenerateGainMap() { //--- Fit every peak energy for every segment --- fSegLinearFit = std::vector(h.size(), std::vector(h.at(0).size(), nullptr)); for (size_t i = 0; i < h.size(); i++) { - for (size_t j = 0; j < h.at(0).size(); j++) { - RESTExtreme << "Segment[" << i << "][" << j << "]" << p->RESTendl; + for (int j = -1; j < (int)h.at(0).size(); j++) // -1 for the whole module + { + TH1F* hSeg = nullptr; + if (i > 0 && j == -1) + continue; + else if (i == 0 && j == -1) { + hSeg = hModule; + RESTExtreme << "Whole module" << p->RESTendl; + } else { + hSeg = h[i][j]; + RESTExtreme << "Segment[" << i << "][" << j << "]" << p->RESTendl; + } + // Search for peaks --> peakPos std::unique_ptr s(new TSpectrum(2 * fEnergyPeaks.size() + 1)); std::vector peakPos; - s->Search(h[i][j], 2, "goff", 0.1); + s->Search(hSeg, 2, "goff", 0.1); for (int k = 0; k < s->GetNPeaks(); k++) peakPos.push_back(s->GetPositionX()[k]); std::sort(peakPos.begin(), peakPos.end(), std::greater()); const double ratio = peakPos.size() == 0 @@ -790,10 +831,10 @@ void TRestDataSetGainMap::Module::GenerateGainMap() { << DoubleToString(start, "%.3g") << ", " << DoubleToString(end, "%.3g") << ")" << p->RESTendl; - if (h[i][j]->GetFunction(name.c_str())) // remove previous fit - h[i][j]->GetListOfFunctions()->Remove(h[i][j]->GetFunction(name.c_str())); + if (hSeg->GetFunction(name.c_str())) // remove previous fit + hSeg->GetListOfFunctions()->Remove(hSeg->GetFunction(name.c_str())); - h[i][j]->Fit(g, "R+Q0"); // use 0 to not draw the fit but save it + hSeg->Fit(g, "R+Q0"); // use 0 to not draw the fit but save it mu = g->GetParameter(1); RESTExtreme << "\t\tgaus mean " << DoubleToString(mu, "%g") << p->RESTendl; } while (fAutoRangePeaks && peakPos.size() > 0 && @@ -814,16 +855,22 @@ void TRestDataSetGainMap::Module::GenerateGainMap() { std::unique_ptr linearFit; linearFit = std::unique_ptr(new TF1("linearFit", "pol1")); gr->Fit("linearFit", "SQ"); // Q for quiet mode - - fSegLinearFit.at(i).at(j) = (TGraph*)gr->Clone(); - const double slope = linearFit->GetParameter(1); const double intercept = linearFit->GetParameter(0); - calParSlope.at(i).at(j) = slope; - calParIntercept.at(i).at(j) = intercept; + + if (j == -1) { + fFullLinearFit = (TGraph*)gr->Clone(); + fFullSlope = slope; + fFullIntercept = intercept; + } else { + fSegLinearFit.at(i).at(j) = (TGraph*)gr->Clone(); + calParSlope.at(i).at(j) = slope; + calParIntercept.at(i).at(j) = intercept; + } } } fSegSpectra = h; + fFullSpectrum = hModule; fSlope = calParSlope; fIntercept = calParIntercept; } @@ -887,6 +934,55 @@ void TRestDataSetGainMap::Module::Refit(const size_t x, const size_t y, const si UpdateCalibrationFits(x, y); } +///////////////////////////////////////////// +/// \brief Function to fit again manually a peak for the whole module spectrum. The +/// calibration curve is updated after the fit. +/// +/// \param energyPeak The energy of the peak to be fitted (in physical units). +/// \param range The range for the fitting of the peak (in the observables corresponding units). +/// +void TRestDataSetGainMap::Module::RefitFullSpc(const double energyPeak, const TVector2& range) { + int peakNumber = -1; + for (size_t i = 0; i < fEnergyPeaks.size(); i++) + if (fEnergyPeaks.at(i) == energyPeak) { + peakNumber = i; + break; + } + if (peakNumber == -1) { + RESTError << "Energy " << energyPeak << " not found in the list of energy peaks" << p->RESTendl; + return; + } + RefitFullSpc((size_t)peakNumber, range); +} + +///////////////////////////////////////////// +/// \brief Function to fit again manually a peak for the whole module spectrum. The +/// calibration curve is updated after the fit. +/// +/// \param peakNumber The index of the peak to be fitted. +/// \param range The range for the fitting of the peak (in the observables corresponding units). +/// +void TRestDataSetGainMap::Module::RefitFullSpc(const size_t peakNumber, const TVector2& range) { + if (!fFullSpectrum) { + RESTError << "No gain map found. Use GenerateGainMap() first." << p->RESTendl; + return; + } + if (peakNumber >= fEnergyPeaks.size()) { + RESTError << "Peak with index " << peakNumber << " not found" << p->RESTendl; + return; + } + + // Refit the desired peak + std::string name = "g" + std::to_string(peakNumber); + TF1* g = new TF1(name.c_str(), "gaus", range.X(), range.Y()); + while (fFullSpectrum->GetFunction(name.c_str())) // clear previous fits for this peakNumber + fFullSpectrum->GetListOfFunctions()->Remove(fFullSpectrum->GetFunction(name.c_str())); + fFullSpectrum->Fit(g, "R+Q0"); // use 0 to not draw the fit but save it + + // Change the point of the graph + UpdateCalibrationFitsFullSpc(); +} + ///////////////////////////////////////////// /// \brief Function to update the calibration curve for a given segment of the module. The calibration /// curve is cleared and then the means of the gaussian fits for each energy peak are added. If there are @@ -941,6 +1037,51 @@ void TRestDataSetGainMap::Module::UpdateCalibrationFits(const size_t x, const si fIntercept.at(x).at(y) = lf->GetParameter(0); } +///////////////////////////////////////////// +/// \brief Function to update the calibration curve for the whole module. The calibration +/// curve is cleared and then the means of the gaussian fits for each energy peak are added. If there are +/// less than 2 fits, zero points are added. Then, the calibration curve is refitted (linearFit). +/// +void TRestDataSetGainMap::Module::UpdateCalibrationFitsFullSpc() { + if (!fFullSpectrum) { + RESTError << "No gain map found. Use GenerateGainMap() first." << p->RESTendl; + return; + } + + TGraph* gr = fFullLinearFit; + + // Clear the points of the graph + for (size_t i = 0; i < fEnergyPeaks.size(); i++) gr->RemovePoint(i); + // Add the new points to the graph + int c = 0; + for (size_t i = 0; i < fEnergyPeaks.size(); i++) { + std::string fitName = (std::string) "g" + std::to_string(i); + TF1* g = fFullSpectrum->GetFunction(fitName.c_str()); + if (!g) { + RESTWarning << "No fit ( " << fitName << " ) found for energy peak " << fEnergyPeaks[i] + << " in whole module" << p->RESTendl; + continue; + } + gr->SetPoint(c++, g->GetParameter(1), fEnergyPeaks[i]); + } + + // Add zero points if needed (if there are less than 2 points) + while (gr->GetN() < 2) { + gr->SetPoint(c++, 0, 0); + RESTWarning << "Not enough points for linear fit at whole module. Adding zero point." << p->RESTendl; + } + + // Refit the calibration curve + TF1* lf = nullptr; + if (gr->GetFunction("linearFit")) + lf = gr->GetFunction("linearFit"); + else + lf = new TF1("linearFit", "pol1"); + gr->Fit(lf, "SQ"); // Q for quiet mode + fFullSlope = lf->GetParameter(1); + fFullIntercept = lf->GetParameter(0); +} + ///////////////////////////////////////////// /// \brief Function to read the parameters from the RML element (TiXmlElement) /// and set those class members. @@ -1109,25 +1250,35 @@ void TRestDataSetGainMap::Module::DrawSpectrum(const bool drawFits, const int co } } } -void TRestDataSetGainMap::Module::DrawFullSpectrum() { - if (fSegSpectra.size() == 0) { - RESTError << "Spectra matrix is empty." << p->RESTendl; +void TRestDataSetGainMap::Module::DrawFullSpectrum(const bool drawFits, const int color, TCanvas* c) { + if (!fFullSpectrum) { + RESTError << "Spectrum is empty." << p->RESTendl; return; } - TH1F* sumHist = - new TH1F("fullSpc", "", fSegSpectra[0][0]->GetNbinsX(), fSegSpectra[0][0]->GetXaxis()->GetXmin(), - fSegSpectra[0][0]->GetXaxis()->GetXmax()); - sumHist->SetTitle(("Full spectrum;" + GetObservable() + ";counts").c_str()); - for (size_t i = 0; i < fSegSpectra.size(); i++) { - for (size_t j = 0; j < fSegSpectra.at(0).size(); j++) { - sumHist->Add(fSegSpectra[i][j]); - } + if (!c) { + std::string t = "fullSpc_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId); + c = new TCanvas(t.c_str(), t.c_str()); } - std::string t = "fullSpc_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId); - TCanvas* c = new TCanvas(t.c_str(), t.c_str()); c->cd(); - sumHist->Draw(); + + fFullSpectrum->SetTitle(("Full spectrum;" + GetObservable() + ";counts").c_str()); + + if (color > 0) fFullSpectrum->SetLineColor(color); + size_t colorT = fFullSpectrum->GetLineColor(); + fFullSpectrum->Draw("same"); + + if (drawFits) + for (size_t c = 0; c < fEnergyPeaks.size(); c++) { + auto fit = fFullSpectrum->GetFunction(("g" + std::to_string(c)).c_str()); + if (!fit) RESTWarning << "Fit for energy peak" << fEnergyPeaks[c] << " not found." << p->RESTendl; + if (!fit) continue; + fit->SetLineColor(c + 2 != colorT ? c + 2 : ++colorT); /* does not work with kRed, kBlue, etc. + as they are not defined with the same number as the first 10 basic colors. See + https://root.cern.ch/doc/master/classTColor.html#C01 and + https://root.cern.ch/doc/master/classTColor.html#C02 */ + fit->Draw("same"); + } } void TRestDataSetGainMap::Module::DrawLinearFit(const TVector2& position, TCanvas* c) { @@ -1200,6 +1351,10 @@ void TRestDataSetGainMap::Module::DrawGainMap(const int peakNumber) { RESTError << "Linear fit matrix is empty." << p->RESTendl; return; } + if (!fFullLinearFit) { + RESTError << "Full linear fit is empty." << p->RESTendl; + return; + } double peakEnergy = fEnergyPeaks[peakNumber]; std::string title = "Gain map for energy " + DoubleToString(peakEnergy, "%g") + ";" + @@ -1211,8 +1366,7 @@ void TRestDataSetGainMap::Module::DrawGainMap(const int peakNumber) { TH2F* hGainMap = new TH2F(("h" + t).c_str(), title.c_str(), fNumberOfSegmentsX, fReadoutRange.X(), fReadoutRange.Y(), fNumberOfSegmentsY, fReadoutRange.X(), fReadoutRange.Y()); - const double peakPosRef = - fSegLinearFit[(fNumberOfSegmentsX - 1) / 2][(fNumberOfSegmentsY - 1) / 2]->GetPointX(peakNumber); + const double peakPosRef = fFullLinearFit->GetPointX(peakNumber); auto itX = fSplitX.begin(); for (size_t i = 0; i < fSegLinearFit.size(); i++) { @@ -1313,5 +1467,9 @@ void TRestDataSetGainMap::Module::Print() const { } RESTMetadata << p->RESTendl; } + + RESTMetadata << " Full slope: " << DoubleToString(fFullSlope, "%.3e") << p->RESTendl; + RESTMetadata << " Full intercept: " << DoubleToString(fFullIntercept, "%+.3e") << p->RESTendl; + RESTMetadata << "-----------------------------------------------" << p->RESTendl; } From 114e4f47149280b7b1b59781eb78ca0fc9436939 Mon Sep 17 00:00:00 2001 From: Alvaro Ezquerro Date: Sat, 30 Mar 2024 22:36:38 +0100 Subject: [PATCH 06/12] Add GetModule(index), GetParent() and minor changes (methods names, messages and nullptr checkings) --- .../analysis/inc/TRestDataSetGainMap.h | 12 ++++--- .../analysis/src/TRestDataSetGainMap.cxx | 36 ++++++++++++++++--- 2 files changed, 39 insertions(+), 9 deletions(-) diff --git a/source/framework/analysis/inc/TRestDataSetGainMap.h b/source/framework/analysis/inc/TRestDataSetGainMap.h index 27a426137..4015b944d 100644 --- a/source/framework/analysis/inc/TRestDataSetGainMap.h +++ b/source/framework/analysis/inc/TRestDataSetGainMap.h @@ -85,6 +85,7 @@ class TRestDataSetGainMap : public TRestMetadata { std::string GetSpatialObservableY() const { return fSpatialObservableY; } TRestCut* GetCut() const { return fCut; } + Module* GetModule(const size_t index = 0); Module* GetModule(const int planeID, const int moduleID); double GetSlopeParameter(const int planeID, const int moduleID, const double x, const double y); double GetInterceptParameter(const int planeID, const int moduleID, const double x, const double y); @@ -93,7 +94,7 @@ class TRestDataSetGainMap : public TRestMetadata { void SetCalibrationFileName(const std::string& fileName) { fCalibFileName = fileName; } void SetOutputFileName(const std::string& fileName) { fOutputFileName = fileName; } - void SetModuleCalibration(const Module& moduleCal); + void SetModule(const Module& moduleCal); void SetObservable(const std::string& observable) { fObservable = observable; } void SetSpatialObservableX(const std::string& spatialObservableX) { fSpatialObservableX = spatialObservableX; @@ -189,12 +190,12 @@ class TRestDataSetGainMap : public TRestMetadata { /// Array containing the observable spectrum for each segment. std::vector> fSegSpectra = {}; //< - /// Spectrum of the observable for the whole module. - TH1F* fFullSpectrum = nullptr; //< - /// Array containing the calibration linear fit for each segment. std::vector> fSegLinearFit = {}; //< + /// Spectrum of the observable for the whole module. + TH1F* fFullSpectrum = nullptr; //< + /// Calibration linear fit for the whole module. TGraph* fFullLinearFit = nullptr; //< @@ -206,6 +207,7 @@ class TRestDataSetGainMap : public TRestMetadata { void LoadConfigFromTiXmlElement(const TiXmlElement* module); + TRestDataSetGainMap* GetParent() const { return const_cast(p); } std::pair GetIndexMatrix(const double x, const double y) const; double GetSlope(const double x, const double y) const; double GetIntercept(const double x, const double y) const; @@ -224,7 +226,7 @@ class TRestDataSetGainMap : public TRestMetadata { std::set GetSplitX() const { return fSplitX; } std::set GetSplitY() const { return fSplitY; } std::string GetDataSetFileName() const { return fDataSetFileName; } - TVector2 GetReadoutRangeVar() const { return fReadoutRange; } + TVector2 GetReadoutRange() const { return fReadoutRange; } void DrawSpectrum(const bool drawFits = true, const int color = -1, TCanvas* c = nullptr); void DrawSpectrum(const TVector2& position, bool drawFits = true, int color = -1, diff --git a/source/framework/analysis/src/TRestDataSetGainMap.cxx b/source/framework/analysis/src/TRestDataSetGainMap.cxx index 2d8883646..e96c209b7 100644 --- a/source/framework/analysis/src/TRestDataSetGainMap.cxx +++ b/source/framework/analysis/src/TRestDataSetGainMap.cxx @@ -196,10 +196,12 @@ void TRestDataSetGainMap::InitFromConfigFile() { /// void TRestDataSetGainMap::GenerateGainMap() { for (auto& mod : fModulesCal) { - RESTInfo << "Calibrating plane " << mod.GetPlaneId() << " module " << mod.GetModuleId() << RESTendl; + RESTInfo << "Generating gain map of plane " << mod.GetPlaneId() << " module " << mod.GetModuleId() + << RESTendl; mod.GenerateGainMap(); if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) { mod.DrawSpectrum(); + mod.DrawFullSpectrum(); mod.DrawGainMap(); } } @@ -276,6 +278,21 @@ void TRestDataSetGainMap::CalibrateDataSet(const std::string& dataSetFileName, s delete f; } +///////////////////////////////////////////// +/// \brief Function to retrieve the module calibration by index. Default is 0. +/// +/// +TRestDataSetGainMap::Module* TRestDataSetGainMap::GetModule(const size_t index) { + if (index < fModulesCal.size()) return &fModulesCal[index]; + + RESTError << "No ModuleCalibration with index " << index; + if (fModulesCal.empty()) + RESTError << ". There are no modules defined." << RESTendl; + else + RESTError << ". Max index is " << fModulesCal.size() - 1 << RESTendl; + return nullptr; +} + ///////////////////////////////////////////// /// \brief Function to retrieve the module calibration with planeID and moduleID /// @@ -380,7 +397,7 @@ TRestDataSetGainMap& TRestDataSetGainMap::operator=(TRestDataSetGainMap& src) { /// \brief Function to set a module calibration. If the module calibration /// already exists (same planeId and moduleId), it will be replaced. /// -void TRestDataSetGainMap::SetModuleCalibration(const Module& moduleCal) { +void TRestDataSetGainMap::SetModule(const Module& moduleCal) { for (auto& i : fModulesCal) { if (i.GetPlaneId() == moduleCal.GetPlaneId() && i.GetModuleId() == moduleCal.GetModuleId()) { i = moduleCal; @@ -1170,6 +1187,10 @@ void TRestDataSetGainMap::Module::DrawSpectrum(const int index_x, const int inde RESTError << "Index out of range." << p->RESTendl; return; } + if (!fSegSpectra[index_x][index_y]) { + RESTError << "No Spectrum for segment (" << index_x << ", " << index_y << ")." << p->RESTendl; + return; + } if (!c) { std::string t = "spectrum_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId) + "_" + @@ -1193,7 +1214,8 @@ void TRestDataSetGainMap::Module::DrawSpectrum(const int index_x, const int inde if (drawFits) for (size_t c = 0; c < fEnergyPeaks.size(); c++) { auto fit = fSegSpectra[index_x][index_y]->GetFunction(("g" + std::to_string(c)).c_str()); - if (!fit) RESTWarning << "Fit for energy peak" << fEnergyPeaks[c] << " not found." << p->RESTendl; + if (!fit) + RESTWarning << "Fit for energy peak " << fEnergyPeaks[c] << " not found." << p->RESTendl; if (!fit) continue; fit->SetLineColor(c + 2 != colorT ? c + 2 : ++colorT); /* does not work with kRed, kBlue, etc. as they are not defined with the same number as the first 10 basic colors. See @@ -1296,6 +1318,11 @@ void TRestDataSetGainMap::Module::DrawLinearFit(const int index_x, const int ind RESTError << "Index out of range." << p->RESTendl; return; } + if (!fSegLinearFit[index_x][index_y]) { + RESTError << "No linear fit for segment (" << index_x << ", " << index_y << ")." << p->RESTendl; + return; + } + if (!c) { std::string t = "linearFit_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId) + "_" + std::to_string(index_x) + "_" + std::to_string(index_y); @@ -1379,6 +1406,7 @@ void TRestDataSetGainMap::Module::DrawGainMap(const int peakNumber) { float xMean = (xUpper + xLower) / 2.; float yMean = (yUpper + yLower) / 2.; auto [index_x, index_y] = GetIndexMatrix(xMean, yMean); + if (!fSegLinearFit[index_x][index_y]) continue; hGainMap->Fill(xMean, yMean, fSegLinearFit[index_x][index_y]->GetPointX(peakNumber) / peakPosRef); itY++; } @@ -1467,7 +1495,7 @@ void TRestDataSetGainMap::Module::Print() const { } RESTMetadata << p->RESTendl; } - + RESTMetadata << p->RESTendl; RESTMetadata << " Full slope: " << DoubleToString(fFullSlope, "%.3e") << p->RESTendl; RESTMetadata << " Full intercept: " << DoubleToString(fFullIntercept, "%+.3e") << p->RESTendl; From 74df7df1bdbf45cf683f164f7c76a7e275088d2f Mon Sep 17 00:00:00 2001 From: Alvaro Ezquerro Date: Sat, 30 Mar 2024 22:40:36 +0100 Subject: [PATCH 07/12] Fixing full module spectrum. Move fitting to new method FitPeaks. --- .../analysis/inc/TRestDataSetGainMap.h | 3 + .../analysis/src/TRestDataSetGainMap.cxx | 237 +++++++++--------- 2 files changed, 120 insertions(+), 120 deletions(-) diff --git a/source/framework/analysis/inc/TRestDataSetGainMap.h b/source/framework/analysis/inc/TRestDataSetGainMap.h index 4015b944d..e5ddf6e9b 100644 --- a/source/framework/analysis/inc/TRestDataSetGainMap.h +++ b/source/framework/analysis/inc/TRestDataSetGainMap.h @@ -125,6 +125,9 @@ class TRestDataSetGainMap : public TRestMetadata { private: /// Pointer to the parent class const TRestDataSetGainMap* p = nullptr; // FitPeaks(TH1F* hSeg, TGraph* gr); + public: /// Plane ID (unique identifier). Although it is not linked to any TRestDetectorReadout it is /// recommended to use the same. diff --git a/source/framework/analysis/src/TRestDataSetGainMap.cxx b/source/framework/analysis/src/TRestDataSetGainMap.cxx index e96c209b7..bcf313139 100644 --- a/source/framework/analysis/src/TRestDataSetGainMap.cxx +++ b/source/framework/analysis/src/TRestDataSetGainMap.cxx @@ -710,28 +710,25 @@ void TRestDataSetGainMap::Module::GenerateGainMap() { // --- Definition of histogram whole module --- std::string hModuleName = "hSpc_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId); - TH1F* hModule = new TH1F(hModuleName.c_str(), "", fNBins, fCalibRange.X(), fCalibRange.Y()); + fFullSpectrum = new TH1F(hModuleName.c_str(), "", fNBins, fCalibRange.X(), fCalibRange.Y()); // build the spectrum for the whole module std::string cut = fDefinitionCut; if (cut.empty()) cut = "1"; auto histoMod = dataSet.GetDataFrame().Filter(cut).Histo1D( {"tempMod", "", fNBins, fCalibRange.X(), fCalibRange.Y()}, GetObservable()); - hModule = (TH1F*)histoMod->Clone(hModuleName.c_str()); + std::unique_ptr hpuntMod = std::unique_ptr(static_cast(histoMod->Clone())); + fFullSpectrum->Add(hpuntMod.get()); //--- Definition of histogram matrix --- std::vector> h(fNumberOfSegmentsX, std::vector(fNumberOfSegmentsY, nullptr)); for (size_t i = 0; i < h.size(); i++) { for (size_t j = 0; j < h.at(0).size(); j++) { - std::string name = hModuleName + std::to_string(i) + "_" + std::to_string(j); + std::string name = hModuleName + "_" + std::to_string(i) + "_" + std::to_string(j); h[i][j] = new TH1F(name.c_str(), "", fNBins, fCalibRange.X(), fCalibRange.Y()); // h[column][row] equivalent to h[x][y] } } - std::vector> calParSlope(fNumberOfSegmentsX, - std::vector(fNumberOfSegmentsY, 0)); - std::vector> calParIntercept(fNumberOfSegmentsX, - std::vector(fNumberOfSegmentsY, 0)); // build the spectrum for each segment auto itX = fSplitX.begin(); @@ -768,128 +765,128 @@ void TRestDataSetGainMap::Module::GenerateGainMap() { } //--- Fit every peak energy for every segment --- + std::vector> calParSlope(fNumberOfSegmentsX, + std::vector(fNumberOfSegmentsY, 0)); + std::vector> calParIntercept(fNumberOfSegmentsX, + std::vector(fNumberOfSegmentsY, 0)); fSegLinearFit = std::vector(h.size(), std::vector(h.at(0).size(), nullptr)); - for (size_t i = 0; i < h.size(); i++) { - for (int j = -1; j < (int)h.at(0).size(); j++) // -1 for the whole module - { - TH1F* hSeg = nullptr; - if (i > 0 && j == -1) - continue; - else if (i == 0 && j == -1) { - hSeg = hModule; - RESTExtreme << "Whole module" << p->RESTendl; - } else { - hSeg = h[i][j]; - RESTExtreme << "Segment[" << i << "][" << j << "]" << p->RESTendl; - } + for (size_t i = 0; i < h.size(); i++) + for (int j = 0; j < (int)h.at(0).size(); j++) { + fSegLinearFit[i][j] = new TGraph(); + auto [intercept, slope] = FitPeaks(h[i][j], fSegLinearFit[i][j]); + calParSlope[i][j] = slope; + calParIntercept[i][j] = intercept; + } + fSlope = calParSlope; + fIntercept = calParIntercept; + fSegSpectra = h; - // Search for peaks --> peakPos - std::unique_ptr s(new TSpectrum(2 * fEnergyPeaks.size() + 1)); - std::vector peakPos; - s->Search(hSeg, 2, "goff", 0.1); - for (int k = 0; k < s->GetNPeaks(); k++) peakPos.push_back(s->GetPositionX()[k]); - std::sort(peakPos.begin(), peakPos.end(), std::greater()); - const double ratio = peakPos.size() == 0 - ? 1 - : peakPos.front() / fEnergyPeaks.front(); // to estimate peak position - - // Initialize graph for linear fit - std::shared_ptr gr; - gr = std::shared_ptr(new TGraph()); - gr->SetName("grFit"); - gr->SetTitle((";" + GetObservable() + ";energy").c_str()); - - // Fit every energy peak - int c = 0; - double mu = 0; - for (const auto& energy : fEnergyPeaks) { - RESTExtreme << "\t fitting energy " << DoubleToString(energy, "%g") << p->RESTendl; - // estimation of the peak position is between start and end - double pos = energy * ratio; - double start = pos * 0.8; - double end = pos * 1.2; - if (fRangePeaks.at(c).X() < fRangePeaks.at(c).Y()) { // if range is defined use it - start = fRangePeaks.at(c).X(); - end = fRangePeaks.at(c).Y(); - } + //--- Fit every peak energy for the whole module --- + fFullLinearFit = new TGraph(); + auto [intercept, slope] = FitPeaks(fFullSpectrum, fFullLinearFit); + fFullSlope = slope; + fFullIntercept = intercept; +} + +std::pair TRestDataSetGainMap::Module::FitPeaks(TH1F* hSeg, TGraph* gr) { + if (!hSeg) { + RESTError << "No histogram for fitting" << p->RESTendl; + return std::make_pair(0, 0); + } + if (hSeg->Integral() == 0) { + RESTError << "Empty spectrum " << hSeg->GetName() << p->RESTendl; + return std::make_pair(0, 0); + } + if (!gr) gr = new TGraph(); + RESTExtreme << "Fitting peaks for " << hSeg->GetName() << p->RESTendl; + + // Search for peaks --> peakPos + std::unique_ptr s(new TSpectrum(2 * fEnergyPeaks.size() + 1)); + std::vector peakPos; + s->Search(hSeg, 2, "goff", 0.1); + for (int k = 0; k < s->GetNPeaks(); k++) peakPos.push_back(s->GetPositionX()[k]); + std::sort(peakPos.begin(), peakPos.end(), std::greater()); + const double ratio = + peakPos.size() == 0 ? 1 : peakPos.front() / fEnergyPeaks.front(); // to estimate peak position + + // Initialize graph for linear fit + gr->SetName("grFit"); + gr->SetTitle((";" + GetObservable() + ";energy").c_str()); - do { - if (fAutoRangePeaks) { - if (peakPos.size() > 0) { - // Find the peak position that is between start and end + // Fit every energy peak + int c = 0; + double mu = 0; + for (const auto& energy : fEnergyPeaks) { + RESTExtreme << "\t fitting energy " << DoubleToString(energy, "%g") << p->RESTendl; + // estimation of the peak position is between start and end + double pos = energy * ratio; + double start = pos * 0.8; + double end = pos * 1.2; + if (fRangePeaks.at(c).X() < fRangePeaks.at(c).Y()) { // if range is defined use it + start = fRangePeaks.at(c).X(); + end = fRangePeaks.at(c).Y(); + } + + do { + if (fAutoRangePeaks) { + if (peakPos.size() > 0) { + // Find the peak position that is between start and end + pos = peakPos.at(0); + while (!(start < pos && pos < end)) { + // if none of the peak position is + // between start and end, use the greater one. + if (pos == peakPos.back()) { pos = peakPos.at(0); - while (!(start < pos && pos < end)) { - // if none of the peak position is - // between start and end, use the greater one. - if (pos == peakPos.back()) { - pos = peakPos.at(0); - break; - } - pos = *std::next(std::find(peakPos.begin(), peakPos.end(), - pos)); // get next peak position - } - peakPos.erase(std::find(peakPos.begin(), peakPos.end(), - pos)); // remove this peak position from the list - // final estimation of the peak range (idem fitting window) with this peak - // position pos - start = pos * 0.8; - end = pos * 1.2; - const double relDist = peakPos.size() > 0 ? (pos - peakPos.front()) / pos : 999; - if (relDist < 0.2) { // if the next peak is too close reduce the window width - start = pos * (1 - relDist / 2); - end = pos * (1 + relDist / 2); - } + break; } + pos = *std::next(std::find(peakPos.begin(), peakPos.end(), + pos)); // get next peak position } - - std::string name = "g" + std::to_string(c); - TF1* g = new TF1(name.c_str(), "gaus", start, end); - RESTExtreme << "\t\tat " << DoubleToString(pos, "%.3g") << ". Range(" - << DoubleToString(start, "%.3g") << ", " << DoubleToString(end, "%.3g") << ")" - << p->RESTendl; - - if (hSeg->GetFunction(name.c_str())) // remove previous fit - hSeg->GetListOfFunctions()->Remove(hSeg->GetFunction(name.c_str())); - - hSeg->Fit(g, "R+Q0"); // use 0 to not draw the fit but save it - mu = g->GetParameter(1); - RESTExtreme << "\t\tgaus mean " << DoubleToString(mu, "%g") << p->RESTendl; - } while (fAutoRangePeaks && peakPos.size() > 0 && - !(start < mu && mu < end)); // avoid small peaks on main peak tail - gr->SetPoint(c++, mu, energy); - } - s.reset(); // delete s; - - if (fZeroPoint) gr->SetPoint(c++, 0, 0); - while (gr->GetN() < 2) { // minimun 2 points needed for linear fit - gr->SetPoint(c++, 0, 0); - SetZeroPoint(true); - RESTDebug << "Not enough points for linear fit. Adding and setting zero point to true" - << p->RESTendl; + peakPos.erase(std::find(peakPos.begin(), peakPos.end(), + pos)); // remove this peak position from the list + // final estimation of the peak range (idem fitting window) with this peak + // position pos + start = pos * 0.8; + end = pos * 1.2; + const double relDist = peakPos.size() > 0 ? (pos - peakPos.front()) / pos : 999; + if (relDist < 0.2) { // if the next peak is too close reduce the window width + start = pos * (1 - relDist / 2); + end = pos * (1 + relDist / 2); + } + } } - // Linear fit - std::unique_ptr linearFit; - linearFit = std::unique_ptr(new TF1("linearFit", "pol1")); - gr->Fit("linearFit", "SQ"); // Q for quiet mode - const double slope = linearFit->GetParameter(1); - const double intercept = linearFit->GetParameter(0); - - if (j == -1) { - fFullLinearFit = (TGraph*)gr->Clone(); - fFullSlope = slope; - fFullIntercept = intercept; - } else { - fSegLinearFit.at(i).at(j) = (TGraph*)gr->Clone(); - calParSlope.at(i).at(j) = slope; - calParIntercept.at(i).at(j) = intercept; - } - } + std::string name = "g" + std::to_string(c); + TF1* g = new TF1(name.c_str(), "gaus", start, end); + RESTExtreme << "\t\tat " << DoubleToString(pos, "%.3g") << ". Range(" + << DoubleToString(start, "%.3g") << ", " << DoubleToString(end, "%.3g") << ")" + << p->RESTendl; + + if (hSeg->GetFunction(name.c_str())) // remove previous fit + hSeg->GetListOfFunctions()->Remove(hSeg->GetFunction(name.c_str())); + + hSeg->Fit(g, "R+Q0"); // use 0 to not draw the fit but save it + mu = g->GetParameter(1); + RESTExtreme << "\t\tgaus mean " << DoubleToString(mu, "%g") << p->RESTendl; + } while (fAutoRangePeaks && peakPos.size() > 0 && + !(start < mu && mu < end)); // avoid small peaks on main peak tail + gr->SetPoint(c++, mu, energy); } - fSegSpectra = h; - fFullSpectrum = hModule; - fSlope = calParSlope; - fIntercept = calParIntercept; + s.reset(); // delete s; + + if (fZeroPoint) gr->SetPoint(c++, 0, 0); + while (gr->GetN() < 2) { // minimun 2 points needed for linear fit + gr->SetPoint(c++, 0, 0); + SetZeroPoint(true); + RESTDebug << "Not enough points for linear fit. Adding and setting zero point to true" << p->RESTendl; + } + + // Linear fit + std::unique_ptr linearFit; + linearFit = std::unique_ptr(new TF1("linearFit", "pol1")); + gr->Fit("linearFit", "SQ"); // Q for quiet mode + + return std::make_pair(linearFit->GetParameter(0), linearFit->GetParameter(1)); } ///////////////////////////////////////////// From 28193694537be88e8413ef6c3e38b17c552ec016 Mon Sep 17 00:00:00 2001 From: Alvaro Ezquerro Date: Sat, 30 Mar 2024 22:42:12 +0100 Subject: [PATCH 08/12] CalibrateDataSet: add excludeColumns and calibrated observable with no segmentation. --- .../analysis/src/TRestDataSetGainMap.cxx | 59 ++++++++++++++++--- 1 file changed, 50 insertions(+), 9 deletions(-) diff --git a/source/framework/analysis/src/TRestDataSetGainMap.cxx b/source/framework/analysis/src/TRestDataSetGainMap.cxx index bcf313139..fff16570f 100644 --- a/source/framework/analysis/src/TRestDataSetGainMap.cxx +++ b/source/framework/analysis/src/TRestDataSetGainMap.cxx @@ -208,14 +208,19 @@ void TRestDataSetGainMap::GenerateGainMap() { } ///////////////////////////////////////////// -/// \brief Function to calibrate a dataset. +/// \brief Function to calibrate a dataset with this gain map. /// /// \param dataSetFileName the name of the root file where the TRestDataSet to be calibrated is stored. /// \param outputFileName the name of the output (root) file where the calibrated TRestDataSet will be -/// exported. If empty, the output file will be named as the input file with the suffix "_cc". E.g. -/// "data/myDataSet.root" -> "data/myDataSet_cc.root". -/// -void TRestDataSetGainMap::CalibrateDataSet(const std::string& dataSetFileName, std::string outputFileName) { +/// exported. If empty, the output file will be named as the input file plus the name of the +/// TRestDataSetGainMap. E.g. "data/myDataSet.root" -> "data/myDataSet_.root". +/// \param excludeColumns a vector of strings with the names of the columns to be excluded from the +/// output file. If empty, all columns will be included. If "all" is in the list, all columns will be +/// excluded except the calibrated observable, the calibrated observable with no segmentation and +/// the plane-module identifier (pmID). +/// +void TRestDataSetGainMap::CalibrateDataSet(const std::string& dataSetFileName, std::string outputFileName, + std::vector excludeColumns) { if (fModulesCal.empty()) { RESTError << "TRestDataSetGainMap::CalibrateDataSet: No modules defined." << RESTendl; return; @@ -258,18 +263,54 @@ void TRestDataSetGainMap::CalibrateDataSet(const std::string& dataSetFileName, s dataFrame = dataFrame.Define(calibObsName, calibrate, {fObservable, fSpatialObservableX, fSpatialObservableY, pmIDname}); + // Define a new column with the calibrated observable for the whole module calibration + auto calibrateFullSpc = [this](double val, int pmID) { + for (auto& m : fModulesCal) { + if (pmID == m.GetPlaneId() * 10 + m.GetModuleId()) + return m.GetSlopeFullSpc() * val + m.GetInterceptFullSpc(); + } + // RESTError << "TRestDataSetGainMap::CalibrateDataSet: Module not found for pmID " << pmID << + // RESTendl; + return std::numeric_limits::quiet_NaN(); + }; + std::string calibObsNameFullSpc = (std::string)GetName() + "_"; + calibObsNameFullSpc += + GetObservable().erase(0, GetObservable().find("_") + 1); // remove the "rawAna_" part + calibObsNameFullSpc += "_NoSegmentation"; + dataFrame = dataFrame.Define(calibObsNameFullSpc, calibrateFullSpc, {fObservable, pmIDname}); + dataSet.SetDataFrame(dataFrame); // Format the output file name and export the dataSet if (outputFileName.empty()) outputFileName = dataSetFileName; - if (outputFileName == dataSetFileName) { // TRestDataSet cannot be overwritten - outputFileName = outputFileName.substr(0, outputFileName.find_last_of(".")) + "_cc."; + std::string gmName = GetName(); + outputFileName = outputFileName.substr(0, outputFileName.find_last_of(".")); // remove extension + outputFileName += "_" + gmName; outputFileName += TRestTools::GetFileNameExtension(dataSetFileName); } - RESTInfo << "Exporting calibrated dataSet to " << outputFileName << RESTendl; - dataSet.Export(outputFileName); + // Export dataset. Exclude columns if requested. + std::set excludeCol; // vector with the explicit column names to be excluded + auto columns = dataSet.GetDataFrame().GetColumnNames(); + // Get the columns to be excluded from the list of columns. It accepts wildcards "*" and "?" + for (auto& eC : excludeColumns) { + if (eC.find("*") != std::string::npos || eC.find("?") != std::string::npos) { + for (auto& c : columns) + if (MatchString(c, eC)) excludeCol.insert(c); + } else + excludeCol.insert(eC); + } + // Remove the calibObsName, calibObsNameFullSpc and pmIDname from the list of columns to be excluded + excludeCol.erase(calibObsName); + excludeCol.erase(calibObsNameFullSpc); + excludeCol.erase(pmIDname); + + RESTDebug << "Excluding columns: " << RESTendl; + for (auto& c : excludeCol) RESTDebug << c << ", "; + RESTDebug << RESTendl; + + dataSet.Export(outputFileName, std::vector(excludeCol.begin(), excludeCol.end())); // Add this TRestDataSetGainMap metadata to the output file TFile* f = TFile::Open(outputFileName.c_str(), "UPDATE"); From f13d7046ab78d8a3229aa2b4656bd7fad3d27bf3 Mon Sep 17 00:00:00 2001 From: Alvaro Ezquerro Date: Wed, 3 Apr 2024 14:08:09 +0200 Subject: [PATCH 09/12] Print cuts with metadata and better debug messages --- .../framework/analysis/src/TRestDataSetGainMap.cxx | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/source/framework/analysis/src/TRestDataSetGainMap.cxx b/source/framework/analysis/src/TRestDataSetGainMap.cxx index fff16570f..9d5c99ecf 100644 --- a/source/framework/analysis/src/TRestDataSetGainMap.cxx +++ b/source/framework/analysis/src/TRestDataSetGainMap.cxx @@ -298,7 +298,7 @@ void TRestDataSetGainMap::CalibrateDataSet(const std::string& dataSetFileName, s if (eC.find("*") != std::string::npos || eC.find("?") != std::string::npos) { for (auto& c : columns) if (MatchString(c, eC)) excludeCol.insert(c); - } else + } else if (std::find(columns.begin(), columns.end(), eC) != columns.end()) excludeCol.insert(eC); } // Remove the calibObsName, calibObsNameFullSpc and pmIDname from the list of columns to be excluded @@ -306,7 +306,7 @@ void TRestDataSetGainMap::CalibrateDataSet(const std::string& dataSetFileName, s excludeCol.erase(calibObsNameFullSpc); excludeCol.erase(pmIDname); - RESTDebug << "Excluding columns: " << RESTendl; + RESTDebug << "Excluding columns: "; for (auto& c : excludeCol) RESTDebug << c << ", "; RESTDebug << RESTendl; @@ -513,8 +513,16 @@ void TRestDataSetGainMap::Export(const std::string& fileName) { void TRestDataSetGainMap::PrintMetadata() { TRestMetadata::PrintMetadata(); RESTMetadata << " Calibration dataset: " << fCalibFileName << RESTendl; - if (fCut) fCut->Print(); + if (fCut) { + RESTMetadata << " Cuts applied: "; + /* print only cutStrings and paramCut because + TRestDataSet::MakeCut() uses fCut->GetCutStrings() and fCut->GetParamCut() */ + for (const auto& cut : fCut->GetCutStrings()) RESTMetadata << cut << ", " << RESTendl; + for (const auto& cut : fCut->GetParamCut()) + RESTMetadata << cut.first << " " << cut.second << ", " << RESTendl; + } RESTMetadata << " Output file: " << fOutputFileName << RESTendl; + RESTMetadata << RESTendl; RESTMetadata << " Number of planes: " << GetNumberOfPlanes() << RESTendl; RESTMetadata << " Number of modules: " << GetNumberOfModules() << RESTendl; RESTMetadata << " Calibration observable: " << fObservable << RESTendl; From 45fab637f9b38eb388edf086e0993c81f9388dda Mon Sep 17 00:00:00 2001 From: Alvaro Ezquerro Date: Wed, 3 Apr 2024 14:09:26 +0200 Subject: [PATCH 10/12] Fixing duplicate SetSplits in GenerateGainMap --- source/framework/analysis/src/TRestDataSetGainMap.cxx | 2 -- 1 file changed, 2 deletions(-) diff --git a/source/framework/analysis/src/TRestDataSetGainMap.cxx b/source/framework/analysis/src/TRestDataSetGainMap.cxx index 9d5c99ecf..6aaaeaa3a 100644 --- a/source/framework/analysis/src/TRestDataSetGainMap.cxx +++ b/source/framework/analysis/src/TRestDataSetGainMap.cxx @@ -725,8 +725,6 @@ void TRestDataSetGainMap::Module::GenerateGainMap() { dataSet.SetDataFrame(dataSet.MakeCut(p->GetCut())); - SetSplits(); - if (fSplitX.empty()) SetSplitX(); if (fSplitY.empty()) SetSplitY(); From 76af6ab05c3784100480c0cd94b73e86ce5efc3e Mon Sep 17 00:00:00 2001 From: Alvaro Ezquerro Date: Wed, 3 Apr 2024 14:27:10 +0200 Subject: [PATCH 11/12] Option for changing reference in DrawGainMap --- .../analysis/inc/TRestDataSetGainMap.h | 2 +- .../analysis/src/TRestDataSetGainMap.cxx | 17 +++++++++++++++-- 2 files changed, 16 insertions(+), 3 deletions(-) diff --git a/source/framework/analysis/inc/TRestDataSetGainMap.h b/source/framework/analysis/inc/TRestDataSetGainMap.h index e5ddf6e9b..53b20639e 100644 --- a/source/framework/analysis/inc/TRestDataSetGainMap.h +++ b/source/framework/analysis/inc/TRestDataSetGainMap.h @@ -242,7 +242,7 @@ class TRestDataSetGainMap : public TRestMetadata { void DrawLinearFit(const TVector2& position, TCanvas* c = nullptr); void DrawLinearFit(const int index_x, const int index_y, TCanvas* c = nullptr); - void DrawGainMap(const int peakNumber = 0); + void DrawGainMap(const int peakNumber = 0, const bool fullModuleAsRef = true); void Refit(const TVector2& position, const double energy, const TVector2& range); void Refit(const size_t x, const size_t y, const size_t peakNumber, const TVector2& range); diff --git a/source/framework/analysis/src/TRestDataSetGainMap.cxx b/source/framework/analysis/src/TRestDataSetGainMap.cxx index 6aaaeaa3a..9477c2e50 100644 --- a/source/framework/analysis/src/TRestDataSetGainMap.cxx +++ b/source/framework/analysis/src/TRestDataSetGainMap.cxx @@ -1412,7 +1412,15 @@ void TRestDataSetGainMap::Module::DrawLinearFit(TCanvas* c) { } } -void TRestDataSetGainMap::Module::DrawGainMap(const int peakNumber) { +///////////////////////////////////////////// +/// \brief Function to draw the relative gain map for a given energy peak of the module. +/// +/// \param peakNumber The index of the peak to be drawn (remember they are in descending order). +/// \param fullModuleAsRef If true, it will use the peak position at the full module spectrum +/// as reference for the gain map. If false, it will use the centered segment of the module +/// as reference. +/// +void TRestDataSetGainMap::Module::DrawGainMap(const int peakNumber, bool fullModuleAsRef) { if (peakNumber < 0 || peakNumber >= (int)fEnergyPeaks.size()) { RESTError << "Peak number out of range (peakNumber should be between 0 and " << fEnergyPeaks.size() - 1 << " )" << p->RESTendl; @@ -1437,7 +1445,12 @@ void TRestDataSetGainMap::Module::DrawGainMap(const int peakNumber) { TH2F* hGainMap = new TH2F(("h" + t).c_str(), title.c_str(), fNumberOfSegmentsX, fReadoutRange.X(), fReadoutRange.Y(), fNumberOfSegmentsY, fReadoutRange.X(), fReadoutRange.Y()); - const double peakPosRef = fFullLinearFit->GetPointX(peakNumber); + double peakPosRef = fFullLinearFit->GetPointX(peakNumber); + if (!fullModuleAsRef) { + int index_x = fNumberOfSegmentsX > 0 ? (fNumberOfSegmentsX - 1) / 2 : 0; + int index_y = fNumberOfSegmentsY > 0 ? (fNumberOfSegmentsY - 1) / 2 : 0; + peakPosRef = fSegLinearFit[index_x][index_y]->GetPointX(peakNumber); + } auto itX = fSplitX.begin(); for (size_t i = 0; i < fSegLinearFit.size(); i++) { From d949dde3ba04b9d3e62416176970f1015bd0c62f Mon Sep 17 00:00:00 2001 From: Alvaro Ezquerro Date: Wed, 3 Apr 2024 14:56:08 +0200 Subject: [PATCH 12/12] Avoid potential memory leaks and duplicated code for UpdateCalibrationX --- .../analysis/inc/TRestDataSetGainMap.h | 1 + .../analysis/src/TRestDataSetGainMap.cxx | 86 ++++++++----------- 2 files changed, 37 insertions(+), 50 deletions(-) diff --git a/source/framework/analysis/inc/TRestDataSetGainMap.h b/source/framework/analysis/inc/TRestDataSetGainMap.h index 53b20639e..617c5f0aa 100644 --- a/source/framework/analysis/inc/TRestDataSetGainMap.h +++ b/source/framework/analysis/inc/TRestDataSetGainMap.h @@ -127,6 +127,7 @@ class TRestDataSetGainMap : public TRestMetadata { const TRestDataSetGainMap* p = nullptr; // FitPeaks(TH1F* hSeg, TGraph* gr); + std::pair UpdateCalibrationFits(TH1F* hSeg, TGraph* gr); public: /// Plane ID (unique identifier). Although it is not linked to any TRestDetectorReadout it is diff --git a/source/framework/analysis/src/TRestDataSetGainMap.cxx b/source/framework/analysis/src/TRestDataSetGainMap.cxx index 9477c2e50..eab96f424 100644 --- a/source/framework/analysis/src/TRestDataSetGainMap.cxx +++ b/source/framework/analysis/src/TRestDataSetGainMap.cxx @@ -757,6 +757,7 @@ void TRestDataSetGainMap::Module::GenerateGainMap() { // --- Definition of histogram whole module --- std::string hModuleName = "hSpc_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId); + delete fFullSpectrum; fFullSpectrum = new TH1F(hModuleName.c_str(), "", fNBins, fCalibRange.X(), fCalibRange.Y()); // build the spectrum for the whole module @@ -829,6 +830,7 @@ void TRestDataSetGainMap::Module::GenerateGainMap() { fSegSpectra = h; //--- Fit every peak energy for the whole module --- + delete fFullLinearFit; fFullLinearFit = new TGraph(); auto [intercept, slope] = FitPeaks(fFullSpectrum, fFullLinearFit); fFullSlope = slope; @@ -844,7 +846,7 @@ std::pair TRestDataSetGainMap::Module::FitPeaks(TH1F* hSeg, TGra RESTError << "Empty spectrum " << hSeg->GetName() << p->RESTendl; return std::make_pair(0, 0); } - if (!gr) gr = new TGraph(); + std::shared_ptr graph = std::shared_ptr(new TGraph()); RESTExtreme << "Fitting peaks for " << hSeg->GetName() << p->RESTendl; // Search for peaks --> peakPos @@ -857,8 +859,8 @@ std::pair TRestDataSetGainMap::Module::FitPeaks(TH1F* hSeg, TGra peakPos.size() == 0 ? 1 : peakPos.front() / fEnergyPeaks.front(); // to estimate peak position // Initialize graph for linear fit - gr->SetName("grFit"); - gr->SetTitle((";" + GetObservable() + ";energy").c_str()); + graph->SetName("grFit"); + graph->SetTitle((";" + GetObservable() + ";energy").c_str()); // Fit every energy peak int c = 0; @@ -917,13 +919,13 @@ std::pair TRestDataSetGainMap::Module::FitPeaks(TH1F* hSeg, TGra RESTExtreme << "\t\tgaus mean " << DoubleToString(mu, "%g") << p->RESTendl; } while (fAutoRangePeaks && peakPos.size() > 0 && !(start < mu && mu < end)); // avoid small peaks on main peak tail - gr->SetPoint(c++, mu, energy); + graph->SetPoint(c++, mu, energy); } s.reset(); // delete s; - if (fZeroPoint) gr->SetPoint(c++, 0, 0); - while (gr->GetN() < 2) { // minimun 2 points needed for linear fit - gr->SetPoint(c++, 0, 0); + if (fZeroPoint) graph->SetPoint(c++, 0, 0); + while (graph->GetN() < 2) { // minimun 2 points needed for linear fit + graph->SetPoint(c++, 0, 0); SetZeroPoint(true); RESTDebug << "Not enough points for linear fit. Adding and setting zero point to true" << p->RESTendl; } @@ -931,8 +933,9 @@ std::pair TRestDataSetGainMap::Module::FitPeaks(TH1F* hSeg, TGra // Linear fit std::unique_ptr linearFit; linearFit = std::unique_ptr(new TF1("linearFit", "pol1")); - gr->Fit("linearFit", "SQ"); // Q for quiet mode + graph->Fit("linearFit", "SQ"); // Q for quiet mode + if (gr) *gr = *(TGraph*)graph->Clone(); // if nullptr is passed, do not copy the graph return std::make_pair(linearFit->GetParameter(0), linearFit->GetParameter(1)); } @@ -1065,6 +1068,25 @@ void TRestDataSetGainMap::Module::UpdateCalibrationFits(const size_t x, const si TH1F* h = fSegSpectra.at(x).at(y); TGraph* gr = fSegLinearFit.at(x).at(y); + auto [intercept, slope] = UpdateCalibrationFits(h, gr); + fSlope[x][y] = slope; + fIntercept[x][y] = intercept; +} + +std::pair TRestDataSetGainMap::Module::UpdateCalibrationFits(TH1F* h, TGraph* gr) { + if (!h) { + RESTError << "No histogram for updating fits" << p->RESTendl; + return std::make_pair(0, 0); + } + if (!gr) { + RESTError << "No graph for updating fits" << p->RESTendl; + return std::make_pair(0, 0); + } + if (h->Integral() == 0) { + RESTError << "Empty spectrum " << h->GetName() << p->RESTendl; + return std::make_pair(0, 0); + } + // Clear the points of the graph for (size_t i = 0; i < fEnergyPeaks.size(); i++) gr->RemovePoint(i); // Add the new points to the graph @@ -1074,7 +1096,7 @@ void TRestDataSetGainMap::Module::UpdateCalibrationFits(const size_t x, const si TF1* g = h->GetFunction(fitName.c_str()); if (!g) { RESTWarning << "No fit ( " << fitName << " ) found for energy peak " << fEnergyPeaks[i] - << " in segment " << x << "," << y << p->RESTendl; + << " in histogram " << h->GetName() << p->RESTendl; continue; } gr->SetPoint(c++, g->GetParameter(1), fEnergyPeaks[i]); @@ -1083,8 +1105,6 @@ void TRestDataSetGainMap::Module::UpdateCalibrationFits(const size_t x, const si // Add zero points if needed (if there are less than 2 points) while (gr->GetN() < 2) { gr->SetPoint(c++, 0, 0); - RESTWarning << "Not enough points for linear fit at segment (" << x << ", " << y - << "). Adding zero point." << p->RESTendl; } // Refit the calibration curve @@ -1094,8 +1114,8 @@ void TRestDataSetGainMap::Module::UpdateCalibrationFits(const size_t x, const si else lf = new TF1("linearFit", "pol1"); gr->Fit(lf, "SQ"); // Q for quiet mode - fSlope.at(x).at(y) = lf->GetParameter(1); - fIntercept.at(x).at(y) = lf->GetParameter(0); + + return std::make_pair(lf->GetParameter(0), lf->GetParameter(1)); } ///////////////////////////////////////////// @@ -1104,43 +1124,9 @@ void TRestDataSetGainMap::Module::UpdateCalibrationFits(const size_t x, const si /// less than 2 fits, zero points are added. Then, the calibration curve is refitted (linearFit). /// void TRestDataSetGainMap::Module::UpdateCalibrationFitsFullSpc() { - if (!fFullSpectrum) { - RESTError << "No gain map found. Use GenerateGainMap() first." << p->RESTendl; - return; - } - - TGraph* gr = fFullLinearFit; - - // Clear the points of the graph - for (size_t i = 0; i < fEnergyPeaks.size(); i++) gr->RemovePoint(i); - // Add the new points to the graph - int c = 0; - for (size_t i = 0; i < fEnergyPeaks.size(); i++) { - std::string fitName = (std::string) "g" + std::to_string(i); - TF1* g = fFullSpectrum->GetFunction(fitName.c_str()); - if (!g) { - RESTWarning << "No fit ( " << fitName << " ) found for energy peak " << fEnergyPeaks[i] - << " in whole module" << p->RESTendl; - continue; - } - gr->SetPoint(c++, g->GetParameter(1), fEnergyPeaks[i]); - } - - // Add zero points if needed (if there are less than 2 points) - while (gr->GetN() < 2) { - gr->SetPoint(c++, 0, 0); - RESTWarning << "Not enough points for linear fit at whole module. Adding zero point." << p->RESTendl; - } - - // Refit the calibration curve - TF1* lf = nullptr; - if (gr->GetFunction("linearFit")) - lf = gr->GetFunction("linearFit"); - else - lf = new TF1("linearFit", "pol1"); - gr->Fit(lf, "SQ"); // Q for quiet mode - fFullSlope = lf->GetParameter(1); - fFullIntercept = lf->GetParameter(0); + auto [intercept, slope] = UpdateCalibrationFits(fFullSpectrum, fFullLinearFit); + fFullSlope = slope; + fFullIntercept = intercept; } /////////////////////////////////////////////