TRestResponse.cxx

/*************************************************************************
 * 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.                  *
 *************************************************************************/
 
#include "TRestResponse.h"
 
ClassImp(TRestResponse);
 
TRestResponse::TRestResponse() { Initialize(); }
 
TRestResponse::TRestResponse(const char* cfgFileName, const std::string& name) : TRestMetadata(cfgFileName) {
    Initialize();
 
    LoadConfigFromFile(fConfigFileName, name);
 
    if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) PrintMetadata();
}
 
TRestResponse::~TRestResponse() {}
 
void TRestResponse::Initialize() { SetSectionName(this->ClassName()); }
 
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;
}
 
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;
    }
 
    if (!fInterpolation) {
        Int_t bin = (Int_t)((input - fOrigin.X()) / fBinSize);
 
        for (std::size_t n = 0; n < fResponseMatrix[bin].size(); n++) {
            Double_t output = fOrigin.Y() + ((double)n + 0.5) * fBinSize;
            Double_t value = fResponseMatrix[bin][n];
 
            std::pair<Double_t, Double_t> outp{output, value};
 
            response.push_back(outp);
        }
 
        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;
}
 
void TRestResponse::PrintResponseMatrix(Int_t fromRow = 0, Int_t toRow = 0) {
    TRestTools::PrintTable(fResponseMatrix, fromRow, toRow);
}
 
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;
    }
    if (fInterpolation) {
        RESTMetadata << " " << RESTendl;
        RESTMetadata << "Interpolation is enabled" << RESTendl;
    } else {
        RESTMetadata << " " << RESTendl;
        RESTMetadata << "Interpolation is disabled" << RESTendl;
    }
    RESTMetadata << "----" << RESTendl;
}