Merge branch 'master' into jgalan_fix_typo

rest-for-physics · Nov 16, 2023 · 5d6e1ff · 5d6e1ff
2 parents 09ccda1 + 7434500
commit 5d6e1ff
Show file tree

Hide file tree

Showing 7 changed files with 1,465 additions and 6 deletions.
diff --git a/cmake/thisREST.cmake b/cmake/thisREST.cmake
@@ -16,22 +16,25 @@ execute_process(
 string(REGEX REPLACE "\n$" "" GEANT4_PATH "${GEANT4_PATH}")
 get_filename_component(GEANT4_BIN_DIR "${GEANT4_PATH}/bin/" REALPATH)
 
+set(g4LibPath "")
+set(loadG4 "")
 if (${REST_G4} MATCHES "ON")
     # https://github.com/rest-for-physics/framework/issues/331
     set(g4LibPath ":${GEANT4_PATH}/lib/")
     set(loadG4
         "\# if geant4.sh script is found we load the same Geant4 version as used in compilation\nif [[ -f \\\"${GEANT4_BIN_DIR}/geant4.sh\\\" ]]; then
     [[ -n \\\"\\\${ZSH_VERSION}\\\" ]] && pushd ${GEANT4_BIN_DIR} > /dev/null\n    source ${GEANT4_BIN_DIR}/geant4.sh\n    [[ -n \\\"\\\${ZSH_VERSION}\\\" ]] && popd > /dev/null\nfi\n"
     )
-else ()
-    set(g4LibPath "")
-    set(loadG4 "")
-endif (${REST_G4} MATCHES "ON")
+endif ()
 
+set(loadMPFR "")
 if (DEFINED MPFR_PATH)
     set(loadMPFR "export LD_LIBRARY_PATH=${MPFR_PATH}/lib:\$LD_LIBRARY_PATH")
-else ()
-    set(loadMPFR "")
+endif ()
+
+set(loadCRY "")
+if (DEFINED REST_CRY_PATH)
+    set(loadCRY "export LD_LIBRARY_PATH=${REST_CRY_PATH}/lib:\$LD_LIBRARY_PATH")
 endif ()
 
 set(loadGarfield "")
@@ -111,6 +114,7 @@ fi
 
 ${loadG4}
 ${loadMPFR}
+${loadCRY}
 ${loadGarfield}
 
 if [ \\\$REST_PATH ] ; then

diff --git a/doc/doxygen/images/drawGainMap.png b/doc/doxygen/images/drawGainMap.png
diff --git a/doc/doxygen/images/drawSpectrum.png b/doc/doxygen/images/drawSpectrum.png
diff --git a/doc/doxygen/images/gainCorrectionComparison.png b/doc/doxygen/images/gainCorrectionComparison.png
diff --git a/examples/calibrationCorrection.rml b/examples/calibrationCorrection.rml
@@ -0,0 +1,24 @@
+<?xml version="1.0"?>
+<TRestDataSetGainMap name="calib" verboseLevel="info">
+        <parameter name="calibFileName" value="myDataSet.root"/>
+        <parameter name="outputFileName" value = "myCalibration.root"/>
+        <parameter name = "observable" value="rawAna_ThresholdIntegral"/>
+        <parameter name = "spatialObservableX" value="hitsAna_xMean"/>
+        <parameter name = "spatialObservableY" value="hitsAna_yMean"/>
+        <module planeId="0" moduleId="0"
+            moduleDefinitionCut="TREXsides_tagId==1"
+            numberOfSegmentsX="5"
+            numberOfSegmentsY="5"
+            readoutRange="(-1,246.24)">
+            <peak energy="22.5" range=""/>
+            <peak energy="8.0"  range=""/>
+        </module>
+        <module planeId="1" moduleId="0"
+            moduleDefinitionCut="TREXsides_tagId==2"
+            numberOfSegmentsX="5"
+            numberOfSegmentsY="5"
+            readoutRange="(-1,246.24)">
+            <peak energy="22.5" range=""/>
+            <peak energy="8.0"  range=""/>
+        </module>
+</TRestDataSetGainMap>
diff --git a/source/framework/analysis/inc/TRestDataSetGainMap.h b/source/framework/analysis/inc/TRestDataSetGainMap.h
@@ -0,0 +1,215 @@
+/*************************************************************************
+ * This file is part of the REST software framework.                     *
+ *                                                                       *
+ * Copyright (C) 2016 GIFNA/TREX (University of Zaragoza)                *
+ * For more information see https://gifna.unizar.es/trex                 *
+ *                                                                       *
+ * REST is free software: you can redistribute it and/or modify          *
+ * it under the terms of the GNU General Public License as published by  *
+ * the Free Software Foundation, either version 3 of the License, or     *
+ * (at your option) any later version.                                   *
+ *                                                                       *
+ * REST is distributed in the hope that it will be useful,               *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of        *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the          *
+ * GNU General Public License for more details.                          *
+ *                                                                       *
+ * You should have a copy of the GNU General Public License along with   *
+ * REST in $REST_PATH/LICENSE.                                           *
+ * If not, see https://www.gnu.org/licenses/.                            *
+ * For the list of contributors see $REST_PATH/CREDITS.                  *
+ *************************************************************************/
+
+#ifndef REST_TRestDataSetGainMap
+#define REST_TRestDataSetGainMap
+
+#include <TCanvas.h>
+#include <TFile.h>
+#include <TGraph.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <TRestStringOutput.h>
+#include <TSpectrum.h>
+#include <TTree.h>
+#include <TVector2.h>
+
+#include "TRestDataSet.h"
+#include "TRestMetadata.h"
+
+/// Metadata class to calculate,store and apply the gain corrected calibration of a group of detectors.
+class TRestDataSetGainMap : public TRestMetadata {
+   public:
+    class Module;
+
+   private:
+    std::string fObservable = "";          //"rawAna_ThresholdIntegral"; //<
+    std::string fSpatialObservableX = "";  //"hitsAna_xMean"; //<
+    std::string fSpatialObservableY = "";  //"hitsAna_yMean"; //<
+
+    std::vector<Module> fModulesCal = {};
+    std::string fCalibFileName = "";
+    std::string fOutputFileName = "";
+
+    void Initialize() override;
+    void InitFromConfigFile() override;
+
+   public:
+    std::set<int> GetPlaneIDs() const;
+    std::set<int> GetModuleIDs(const int planeId) const;
+    std::map<int, std::set<int>> GetModuleIDs() const;
+    Int_t GetNumberOfPlanes() const { return GetPlaneIDs().size(); }
+    Int_t GetNumberOfModules() const {
+        int sum = 0;
+        for (auto pID : GetPlaneIDs()) sum += GetModuleIDs(pID).size();
+        return sum;
+    }
+
+    std::string GetCalibrationFileName() const { return fCalibFileName; }
+    std::string GetOutputFileName() const { return fOutputFileName; }
+    std::string GetObservable() const { return fObservable; }
+    std::string GetSpatialObservableX() const { return fSpatialObservableX; }
+    std::string GetSpatialObservableY() const { return fSpatialObservableY; }
+
+    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);
+
+    void SetCalibrationFileName(const std::string& fileName) { fCalibFileName = fileName; }
+    void SetOutputFileName(const std::string& fileName) { fOutputFileName = fileName; }
+    void SetModuleCalibration(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 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 = "");
+
+    TRestDataSetGainMap();
+    TRestDataSetGainMap(const char* configFilename, std::string name = "");
+    ~TRestDataSetGainMap();
+
+    ClassDefOverride(TRestDataSetGainMap, 1);
+
+    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 = {};      //<
+        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 = {};
+
+       public:
+        void AddPeak(const double& energyPeak, const TVector2& rangePeak = TVector2(0, 0)) {
+            fEnergyPeaks.push_back(energyPeak);
+            fRangePeaks.push_back(rangePeak);
+        }
+
+        void LoadConfigFromTiXmlElement(const TiXmlElement* module);
+
+        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;
+
+        Int_t GetPlaneId() const { return fPlaneId; }
+        Int_t GetModuleId() const { return fModuleId; }
+        std::string GetObservable() const { return p->fObservable; }
+        std::string GetSpatialObservableX() const { return p->fSpatialObservableX; }
+        std::string GetSpatialObservableY() const { return p->fSpatialObservableY; }
+        inline std::string GetModuleDefinitionCut() const { return fDefinitionCut; }
+        Int_t GetNumberOfSegmentsX() const { return fNumberOfSegmentsX; }
+        Int_t GetNumberOfSegmentsY() const { return fNumberOfSegmentsY; }
+
+        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; }
+
+        void DrawSpectrum(bool drawFits = true, int color = -1, TCanvas* c = nullptr);
+        void DrawSpectrum(const double x, const double y, bool drawFits = true, int color = -1,
+                          TCanvas* c = nullptr);
+        void DrawSpectrum(const size_t index_x, const size_t index_y, bool drawFits = true, int color = -1,
+                          TCanvas* c = nullptr);
+        void DrawFullSpectrum();
+
+        void DrawLinearFit();
+        void DrawLinearFit(const double x, const double y, TCanvas* c = nullptr);
+        void DrawLinearFit(const size_t index_x, const size_t index_y, TCanvas* c = nullptr);
+
+        void DrawGainMap(const int peakNumber = 0);
+
+        void Refit(const double x, const double y, const double energy, const TVector2& range);
+        void Refit(const int x, const int y, const int peakNumber, const TVector2& range);
+
+        void SetPlaneId(const Int_t& planeId) { fPlaneId = planeId; }
+        void SetModuleId(const Int_t& moduleId) { fModuleId = moduleId; }
+        void SetModuleDefinitionCut(const std::string& moduleDefinitionCut) {
+            fDefinitionCut = moduleDefinitionCut;
+        }
+        void SetCalibrationRange(const TVector2& calibrationRange) { fCalibRange = calibrationRange; }
+        void SetNBins(const Int_t& nBins) { fNBins = nBins; }
+        void SetSplitX();
+        void SetSplitY();
+        void SetSplitX(const std::set<double>& splitX);
+        void SetSplitY(const std::set<double>& splitY);
+        void SetSplits();
+        void SetSplits(const std::set<double>& splitXandY) {
+            SetSplitX(splitXandY);
+            SetSplitY(splitXandY);
+        }
+
+        void SetNumberOfSegmentsX(const Int_t& numberOfSegmentsX) { fNumberOfSegmentsX = numberOfSegmentsX; }
+        void SetNumberOfSegmentsY(const Int_t& numberOfSegmentsY) { fNumberOfSegmentsY = numberOfSegmentsY; }
+
+        void SetDataSetFileName(const std::string& dataSetFileName) { fDataSetFileName = dataSetFileName; }
+        void SetReadoutRange(const TVector2& readoutRange) { fReadoutRange = readoutRange; }
+        void SetZeroPoint(const bool& ZeroPoint) { fZeroPoint = ZeroPoint; }
+        void SetAutoRangePeaks(const bool& autoRangePeaks) { fAutoRangePeaks = autoRangePeaks; }
+
+        void Print() const;
+
+        void GenerateGainMap();
+        void Initialize();
+
+        Module() {}
+        Module(const TRestDataSetGainMap& parent) : p(&parent){};
+        Module(const TRestDataSetGainMap& parent, const Int_t planeId, const Int_t moduleId) : p(&parent) {
+            SetPlaneId(planeId);
+            SetModuleId(moduleId);
+        };
+        ~Module(){};
+    };
+};
+#endif