TRestResponse implementation (WIP)

source/framework/sensitivity/inc/TRestComponent.h

@@ -26,6 +26,7 @@
 #include <TCanvas.h>
 #include <THn.h>
 
+#include "TRestResponse.h"
 #include "TRestDataSet.h"
 #include "TRestMetadata.h"
 
@@ -50,6 +51,9 @@ class TRestComponent : public TRestMetadata {
     /// It is used to define the node that will be accessed for rate retrieval
     Int_t fActiveNode = -1;  //<
 
+    ///
+    TRestResponse* fResponse = nullptr;  //<
+
     /// A canvas for drawing the active node component
     TCanvas* fCanvas = nullptr;  //!
 
@@ -69,15 +73,19 @@ class TRestComponent : public TRestMetadata {
 
    public:
     Double_t GetNormalizedRate(std::vector<Double_t> point);
+    Double_t GetResponseRate(std::vector<Double_t> point);
+
     virtual Double_t GetRate(std::vector<Double_t> point) = 0;
     virtual Double_t GetTotalRate() = 0;
 
     size_t GetDimensions() { return fVariables.size(); }
 
-    Int_t SetActiveNode(Double_t node);
     Int_t GetActiveNode() { return fActiveNode; }
+    Int_t SetActiveNode(Double_t node);
     Double_t GetActiveNodeValue() { return fParameterizationNodes[fActiveNode]; }
 
+    void LoadResponse(const TRestResponse& resp);
+
     void PrintMetadata() override;
 
     void PrintStatistics();
@@ -88,6 +96,6 @@ class TRestComponent : public TRestMetadata {
     TRestComponent();
     ~TRestComponent();
 
-    ClassDefOverride(TRestComponent, 1);
+    ClassDefOverride(TRestComponent, 2);
 };
 #endif

source/framework/sensitivity/inc/TRestResponse.h

@@ -25,16 +25,59 @@
 
 #include "TRestMetadata.h"
 
-/// A response matrix that might be applied to a given signal variable inside TRestResponse
+/// A response matrix that might be applied to a given component inside a TRestComponent
 class TRestResponse : public TRestMetadata {
    private:
-    // TODO Add here the response matrix. Probably a TH2D
+    /// The filename used to import the response matrix
+    std::string fFilename = "";
+
+    /// It defines the variable name for which the response should be applied to
+    std::string fVariable = "";
+
+    /// First element of the response matrix (input/incident, output/detected)
+    TVector2 fOrigin = TVector2(0, 0);
+
+    /// The resolution of the response matrix (binning)
+    Double_t fBinSize = 0.1;  //<
+
+    /// The response matrix
+    std::vector<std::vector<Float_t>> fResponseMatrix;  //<
+
+    /// Determines if the response matrix has been transposed
+    Bool_t fTransposed = false;
 
    public:
+    void SetBinSize(Double_t bSize) { fBinSize = bSize; }
+    void SetResponseFilename(std::string responseFile) { fFilename = responseFile; }
+    void SetOrigin(const TVector2& v) { fOrigin = v; }
+    void SetVariable(const std::string& var) { fVariable = var; }
+
+    Double_t GetBinSize() const { return fBinSize; }
+    std::string GetResponseFilename() const { return fFilename; }
+    TVector2 GetOrigin() const { return fOrigin; }
+    std::string GetVariable() const { return fVariable; }
+
+    TVector2 GetInputRange() const {
+        return TVector2(fOrigin.X(), fOrigin.X() + fResponseMatrix[0].size() * fBinSize);
+    }
+
+    TVector2 GetOutputRange() const {
+        return TVector2(fOrigin.Y(), fOrigin.Y() + fResponseMatrix.size() * fBinSize);
+    }
+
     void Initialize() override;
 
+    void LoadResponse(Bool_t transpose = true);
+
+    std::vector<std::pair<Double_t, Double_t>> GetResponse(Double_t input);
+
+    void PrintResponseMatrix(Int_t fromRow, Int_t toRow);
+
     void PrintMetadata() override;
 
+    std::vector<std::vector<Float_t>> GetMatrix() const { return fResponseMatrix; }
+
+    TRestResponse(const char* cfgFileName, const std::string& name = "");
     TRestResponse();
     ~TRestResponse();
 

source/framework/sensitivity/src/TRestComponent.cxx

@@ -113,6 +113,37 @@ Double_t TRestComponent::GetNormalizedRate(std::vector<Double_t> point) {
     return normFactor * GetRate(point);
 }
 
+///////////////////////////////////////////////
+/// \brief It returns the intensity/rate (in seconds) corresponding to the
+/// generated distribution or formula evaluated at the position of the parameter
+/// space given by point.
+///
+/// The rate returned by the TRestComponent::GetRate method will be normalized
+/// to the corresponding parameter space. Thus, if the parameter consists of
+/// 2-spatial dimensions and 1-energy dimension, the returned rate will be
+/// expressed in standard REST units as, s-1 mm-2 keV-1.
+///
+Double_t TRestComponent::GetResponseRate(std::vector<Double_t> point) {
+    Double_t normFactor = 1;
+    for (size_t n = 0; n < GetDimensions(); n++) normFactor *= fNbins[n] / (fRanges[n].Y() - fRanges[n].X());
+
+    return normFactor * GetRate(point);
+}
+
+///////////////////////////////////////////////
+/// \brief
+///
+void TRestComponent::LoadResponse(const TRestResponse& resp) {
+    if (fResponse) {
+        delete fResponse;
+        fResponse = nullptr;
+    }
+
+    fResponse = (TRestResponse*)resp.Clone("response");
+
+    fResponse->PrintMetadata();
+}
+
 /////////////////////////////////////////////
 /// \brief Prints on screen the information about the metadata members of TRestAxionSolarFlux
 ///

source/framework/sensitivity/src/TRestResponse.cxx

@@ -30,7 +30,7 @@
 ///
 /// History of developments:
 ///
-/// 2022-December: First implementation of TRestResponse
+/// 2023-December: First implementation of TRestResponse
 /// Javier Galan
 ///
 /// \class TRestResponse
@@ -47,6 +47,28 @@ ClassImp(TRestResponse);
 ///
 TRestResponse::TRestResponse() { Initialize(); }
 
+/////////////////////////////////////////////
+/// \brief Constructor loading data from a config file
+///
+/// If no configuration path is defined using TRestMetadata::SetConfigFilePath
+/// the path to the config file must be specified using full path, absolute or
+/// relative.
+///
+/// The default behaviour is that the config file must be specified with
+/// full path, absolute or relative.
+///
+/// \param cfgFileName A const char* giving the path to an RML file.
+/// \param name The name of the specific metadata. It will be used to find the
+/// corresponding TRestAxionMagneticField section inside the RML.
+///
+TRestResponse::TRestResponse(const char* cfgFileName, const std::string& name) : TRestMetadata(cfgFileName) {
+    Initialize();
+
+    LoadConfigFromFile(fConfigFileName, name);
+
+    if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) PrintMetadata();
+}
+
 ///////////////////////////////////////////////
 /// \brief Default destructor
 ///
@@ -58,11 +80,146 @@ TRestResponse::~TRestResponse() {}
 ///
 void TRestResponse::Initialize() { SetSectionName(this->ClassName()); }
 
+///////////////////////////////////////////////
+/// \brief It loads into the fResponseMatrix data member the response from a file.
+///
+/// For the moment only binary data files with format .N???f have been implemented.
+///
+/// The response file should be arranged in a way that the more internal `std::vector`
+/// (row) inside the `std::vector <std::vector>` table corresponds to a specific
+/// transformed energy (output). For example, if we have the incident particle energy
+/// (input) and the expected detected energy response (output) for that energy, then
+/// each row in the matrix corresponds to a detected energy and each element of that
+/// row defines the fraction of incident energies that contributed to that detected
+/// energy.
+///
+/// Thats why we may need to transpose the matrix, so that when we can extract a row
+/// (detected energy) from the matrix directly, such as fMatrix[n], and we just get
+/// the vector that should convolute with the Signal(Ei) that is a vector of signals
+/// as a function of incident energy. The resulting scalar product will give the
+/// expected signal at the detection energy.
+///
+void TRestResponse::LoadResponse(Bool_t transpose) {
+    if (fFilename == "") {
+        RESTError << "TRestResponse::LoadResponse. The response filename was not defined" << RESTendl;
+        return;
+    }
+
+    std::string fullFilename = SearchFile(fFilename);
+    if (fullFilename.empty()) {
+        RESTError << "TRestResponse::LoadResponse. The response filename was not found!" << RESTendl;
+        RESTError << "Filename : " << fFilename << RESTendl;
+        RESTError << "You may want to define <seachPath inside <globals> definition" << RESTendl;
+        return;
+    }
+
+    std::string extension = TRestTools::GetFileNameExtension(fFilename);
+    if (!extension.empty() && extension[0] == 'N' && extension.back() == 'f') {
+        TRestTools::ReadBinaryTable(fullFilename, fResponseMatrix);
+
+        fTransposed = false;
+        if (transpose) {
+            fTransposed = transpose;
+            TRestTools::TransposeTable(fResponseMatrix);
+        }
+
+        return;
+    }
+
+    RESTError << "Extension format - " << extension << " - not recognized!" << RESTendl;
+}
+
+/////////////////////////////////////////////
+/// \brief This method will return a vector of std::pair, each pair will contain the
+/// output/frequency/energy value for the corresponding response.
+///
+/// The output value will be mapped following the binning and the origin given on the
+/// metadata members.
+///
+std::vector<std::pair<Double_t, Double_t>> TRestResponse::GetResponse(Double_t input) {
+
+    std::vector<std::pair<Double_t, Double_t>> response;
+
+    if (fResponseMatrix.empty()) {
+        RESTError << "TRestResponse::GetResponse. Response matrix has not been loaded yet!" << RESTendl;
+        return response;
+    }
+
+    if (input < GetInputRange().X() || input > GetInputRange().Y()) {
+        RESTError << "TRestResponse::GetResponse. The input value " << input << " is outside range!"
+                  << RESTendl;
+        return response;
+    }
+
+    Int_t binLeft = (Int_t)((input - fBinSize / 2. - fOrigin.X()) / fBinSize);
+    Int_t binRight = binLeft + 1;
+
+    Double_t distLeft = (input - fBinSize / 2. + fOrigin.X()) - binLeft * fBinSize;
+
+    if (input <= GetInputRange().X() + fBinSize / 2. || input >= GetInputRange().Y() - fBinSize / 2.)
+        binRight = binLeft;
+
+    /*
+    std::cout << "Top : " << GetInputRange().Y() - fBinSize/2. << std::endl;
+    std::cout << "binLeft : " << binLeft << std::endl;
+    std::cout << "binRight : " << binRight << std::endl;
+    std::cout << "dLeft : " << distLeft << std::endl;
+    std::cout << "dLeft/fBinSize : " << distLeft/fBinSize << std::endl;
+    std::cout << "1 - distLeft/fBinSize : " << 1 - distLeft/fBinSize << std::endl;
+    */
+
+    for (std::size_t n = 0; n < fResponseMatrix[binLeft].size(); n++) {
+        Double_t output = fOrigin.Y() + ((double)n + 0.5) * fBinSize;
+        Double_t value = fResponseMatrix[binLeft][n] * (1 - distLeft / fBinSize) +
+                         fResponseMatrix[binRight][n] * distLeft / fBinSize;
+
+        std::pair<Double_t, Double_t> outp{output, value};
+
+        response.push_back(outp);
+
+        /*
+        std::cout << "n: " << n << " output : " << output << std::endl;
+        std::cout << "response: " << response << std::endl;
+        */
+    }
+
+    return response;
+}
+
+/////////////////////////////////////////////
+/// \brief Prints on screen the information about the metadata members of TRestAxionSolarFlux
+///
+void TRestResponse::PrintResponseMatrix(Int_t fromRow = 0, Int_t toRow = 0) {
+    TRestTools::PrintTable(fResponseMatrix, fromRow, toRow);
+}
+
 /////////////////////////////////////////////
 /// \brief Prints on screen the information about the metadata members of TRestAxionSolarFlux
 ///
 void TRestResponse::PrintMetadata() {
     TRestMetadata::PrintMetadata();
 
+    RESTMetadata << "Response file : " << fFilename << RESTendl;
+    RESTMetadata << "Variable : " << fVariable << RESTendl;
+    RESTMetadata << "Bin size : " << fBinSize << RESTendl;
+    RESTMetadata << " " << RESTendl;
+
+    if (!fResponseMatrix.empty()) {
+        RESTMetadata << "Response matrix has been loaded" << RESTendl;
+        RESTMetadata << " - Number of columns: " << fResponseMatrix[0].size() << RESTendl;
+        RESTMetadata << " - Number of rows : " << fResponseMatrix.size() << RESTendl;
+        RESTMetadata << " - Input range : " << GetInputRange().X() << " - " << GetInputRange().Y()
+                     << RESTendl;
+        RESTMetadata << " - Output range : " << GetOutputRange().X() << " - " << GetOutputRange().Y()
+                     << RESTendl;
+
+        if (fTransposed) {
+            RESTMetadata << " " << RESTendl;
+            RESTMetadata << "Original matrix was transposed" << RESTendl;
+        }
+    } else {
+        RESTMetadata << "Response matrix has NOT been loaded" << RESTendl;
+        RESTMetadata << "Try calling TRestResponse::LoadResponse()" << RESTendl;
+    }
     RESTMetadata << "----" << RESTendl;
 }