Skip to content

Commit

Permalink
TRestDataSetGainMap improvements and bug fixes (#512)
Browse files Browse the repository at this point in the history 
* Documentation fix at header

* EnableMultiThreading when handling TRestDataSet

* Add a TRestCut

* Fix bugs in Draw methods

* Adding full module spectrum (new attributes and methods) and use it as reference for DrawGainMap

* Add GetModule(index), GetParent() and minor changes (methods names, messages and nullptr checkings)

* Fixing full module spectrum. Move fitting to new method FitPeaks.

* CalibrateDataSet: add excludeColumns and calibrated observable with no segmentation.

* Print cuts with metadata and better debug messages

* Fixing duplicate SetSplits in GenerateGainMap

* Option for changing reference in DrawGainMap

* Avoid potential memory leaks and duplicated code for UpdateCalibrationX
  • Loading branch information
@AlvaroEzq
AlvaroEzq authored Apr 8, 2024
1 parent 75b29bb commit 7a4a132
Show file tree
Hide file tree
Showing 2 changed files with 526 additions and 187 deletions.
154 changes: 116 additions & 38 deletions source/framework/analysis/inc/TRestDataSetGainMap.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,7 @@
#include <TTree.h>
#include <TVector2.h>

#include "TRestCut.h"
#include "TRestDataSet.h"
#include "TRestMetadata.h"

Expand All @@ -42,13 +43,26 @@ 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 = ""; //<

std::vector<Module> fModulesCal = {};
std::string fCalibFileName = "";
std::string fOutputFileName = "";
/// 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<Module> 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 = ""; //<

/// Cut to be applied to the calibration data
TRestCut* fCut = nullptr; //<

void Initialize() override;
void InitFromConfigFile() override;
Expand All @@ -69,67 +83,125 @@ 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 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);
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; }
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;
}
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 = "");

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<std::string> excludeColumns = {});

TRestDataSetGainMap();
TRestDataSetGainMap(const char* configFilename, std::string name = "");
~TRestDataSetGainMap();

ClassDefOverride(TRestDataSetGainMap, 1);
ClassDefOverride(TRestDataSetGainMap, 2);

class Module {
private:
const TRestDataSetGainMap* p = nullptr; //<! Pointer to the parent class
public: /// Members that will be written to the ROOT file.
Int_t fPlaneId = -1; //< // Plane ID
Int_t fModuleId = -1; //< // Module ID

std::vector<double> fEnergyPeaks = {};
std::vector<TVector2> 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<double> fSplitX = {}; //<
std::set<double> fSplitY = {}; //<

std::string fDataSetFileName = ""; //< // File name for the calibration dataset

std::vector<std::vector<double>> fSlope = {}; //<
/// Pointer to the parent class
const TRestDataSetGainMap* p = nullptr; //<!

std::pair<double, double> FitPeaks(TH1F* hSeg, TGraph* gr);
std::pair<double, double> UpdateCalibrationFits(TH1F* hSeg, TGraph* gr);

public:
/// Plane ID (unique identifier). Although it is not linked to any TRestDetectorReadout it is
/// recommended to use the same.
Int_t fPlaneId = -1; //<

// Module ID (unique identifier). Although it is not linked to any TRestDetectorReadout it is
// recommended to use the same.
Int_t fModuleId = -1; //<

/// Energy of the peaks to be used for the calibration.
std::vector<double> fEnergyPeaks = {}; //<

/// Range of the peaks to be used for the calibration. If empty it will be automatically calculated.
std::vector<TVector2> 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<double> fSplitX = {}; //<

/// Split points in the y direction.
std::set<double> 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<std::vector<double>> 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<std::vector<double>> 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<std::vector<TH1F*>> fSegSpectra = {};
std::vector<std::vector<TGraph*>> 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<std::vector<TH1F*>> fSegSpectra = {}; //<

/// Array containing the calibration linear fit for each segment.
std::vector<std::vector<TGraph*>> fSegLinearFit = {}; //<

/// Spectrum of the observable for the whole module.
TH1F* fFullSpectrum = nullptr; //<

/// Calibration linear fit for the whole module.
TGraph* fFullLinearFit = nullptr; //<

public:
void AddPeak(const double& energyPeak, const TVector2& rangePeak = TVector2(0, 0)) {
Expand All @@ -139,9 +211,12 @@ class TRestDataSetGainMap : public TRestMetadata {

void LoadConfigFromTiXmlElement(const TiXmlElement* module);

TRestDataSetGainMap* GetParent() const { return const_cast<TRestDataSetGainMap*>(p); }
std::pair<int, int> 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; }
Expand All @@ -155,24 +230,27 @@ class TRestDataSetGainMap : public TRestMetadata {
std::set<double> GetSplitX() const { return fSplitX; }
std::set<double> 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,
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();
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);

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);
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; }
Expand Down
Loading

0 comments on commit 7a4a132

Please sign in to comment.