TRestDataSetGainMap.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 "TRestDataSetGainMap.h"
 
ClassImp(TRestDataSetGainMap);
TRestDataSetGainMap::TRestDataSetGainMap() { Initialize(); }
 
TRestDataSetGainMap::TRestDataSetGainMap(const char* configFilename, std::string name)
    : TRestMetadata(configFilename) {
    LoadConfigFromFile(fConfigFileName, name);
    if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) PrintMetadata();
}
 
TRestDataSetGainMap::~TRestDataSetGainMap() {}
 
void TRestDataSetGainMap::Initialize() { SetSectionName(this->ClassName()); }
void TRestDataSetGainMap::InitFromConfigFile() {
    this->Initialize();
    TRestMetadata::InitFromConfigFile();
 
    // Load Module from RML
    TiXmlElement* moduleDefinition = GetElement("module");
    while (moduleDefinition != nullptr) {
        fModulesCal.push_back(Module(*this));
        fModulesCal.back().LoadConfigFromTiXmlElement(moduleDefinition);
 
        moduleDefinition = GetNextElement(moduleDefinition);
    }
 
    fCut = (TRestCut*)InstantiateChildMetadata("TRestCut");
 
    if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) PrintMetadata();
}
 
void TRestDataSetGainMap::GenerateGainMap() {
    for (auto& mod : fModulesCal) {
        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();
        }
    }
}
 
void TRestDataSetGainMap::CalibrateDataSet(const std::string& dataSetFileName, std::string outputFileName,
                                           std::vector<std::string> excludeColumns) {
    if (fModulesCal.empty()) {
        RESTError << "TRestDataSetGainMap::CalibrateDataSet: No modules defined." << RESTendl;
        return;
    }
 
    TRestDataSet dataSet;
    dataSet.EnableMultiThreading(true);
    dataSet.Import(dataSetFileName);
    auto dataFrame = dataSet.GetDataFrame();
 
    // Define a new column with the identifier (pmID) of the module for each row (event)
    std::string pmIDname = (std::string)GetName() + "_pmID";
    std::string modCut = fModulesCal[0].GetModuleDefinitionCut();
    int pmID = fModulesCal[0].GetPlaneId() * 10 + fModulesCal[0].GetModuleId();
 
    auto columnList = dataFrame.GetColumnNames();
    if (std::find(columnList.begin(), columnList.end(), pmIDname) == columnList.end())
        dataFrame = dataFrame.Define(pmIDname, modCut + " ? " + std::to_string(pmID) + " : -1");
    else
        dataFrame = dataFrame.Redefine(pmIDname, modCut + " ? " + std::to_string(pmID) + " : -1");
 
    for (size_t n = 1; n < fModulesCal.size(); n++) {
        modCut = fModulesCal[n].GetModuleDefinitionCut();
        pmID = fModulesCal[n].GetPlaneId() * 10 + fModulesCal[n].GetModuleId();
        dataFrame = dataFrame.Redefine(pmIDname, (modCut + " ? " + std::to_string(pmID) + " : " + pmIDname));
    }
 
    // Define a new column with the calibrated observable
    auto calibrate = [this](double val, double x, double y, int pmID) {
        for (auto& m : fModulesCal) {
            if (pmID == m.GetPlaneId() * 10 + m.GetModuleId())
                return m.GetSlope(x, y) * val + m.GetIntercept(x, y);
        }
        // RESTError << "TRestDataSetGainMap::CalibrateDataSet: Module not found for pmID " << pmID <<
        // RESTendl;
        return std::numeric_limits<double>::quiet_NaN();
    };
    std::string calibObsName = (std::string)GetName() + "_";
    calibObsName += GetObservable().erase(0, GetObservable().find("_") + 1);  // remove the "rawAna_" part
    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<double>::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
        std::string gmName = GetName();
        outputFileName = outputFileName.substr(0, outputFileName.find_last_of("."));  // remove extension
        outputFileName += "_" + gmName;
        outputFileName += TRestTools::GetFileNameExtension(dataSetFileName);
    }
 
    // Export dataset. Exclude columns if requested.
    std::set<std::string> 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 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
    excludeCol.erase(calibObsName);
    excludeCol.erase(calibObsNameFullSpc);
    excludeCol.erase(pmIDname);
 
    RESTDebug << "Excluding columns: ";
    for (auto& c : excludeCol) RESTDebug << c << ", ";
    RESTDebug << RESTendl;
 
    dataSet.Export(outputFileName, std::vector<std::string>(excludeCol.begin(), excludeCol.end()));
 
    // Add this TRestDataSetGainMap metadata to the output file
    TFile* f = TFile::Open(outputFileName.c_str(), "UPDATE");
    this->Write();
    f->Close();
    delete f;
}
 
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;
}
 
TRestDataSetGainMap::Module* TRestDataSetGainMap::GetModule(const int planeID, const int moduleID) {
    for (auto& i : fModulesCal) {
        if (i.GetPlaneId() == planeID && i.GetModuleId() == moduleID) return &i;
    }
    RESTError << "No ModuleCalibration with planeID " << planeID << " and moduleID " << moduleID << RESTendl;
    return nullptr;
}
 
double TRestDataSetGainMap::GetSlopeParameter(const int planeID, const int moduleID, const double x,
                                              const double y) {
    Module* moduleCal = GetModule(planeID, moduleID);
    if (moduleCal == nullptr) return 0;  // return numeric_limits<double>::quiet_NaN()
    return moduleCal->GetSlope(x, y);
}
 
double TRestDataSetGainMap::GetSlopeParameterFullSpc(const int planeID, const int moduleID) {
    Module* moduleCal = GetModule(planeID, moduleID);
    if (moduleCal == nullptr) return 0;  // return numeric_limits<double>::quiet_NaN()
    return moduleCal->GetSlopeFullSpc();
}
 
double TRestDataSetGainMap::GetInterceptParameter(const int planeID, const int moduleID, const double x,
                                                  const double y) {
    Module* moduleCal = GetModule(planeID, moduleID);
    if (moduleCal == nullptr) return 0;  // return numeric_limits<double>::quiet_NaN()
    return moduleCal->GetIntercept(x, y);
}
 
double TRestDataSetGainMap::GetInterceptParameterFullSpc(const int planeID, const int moduleID) {
    Module* moduleCal = GetModule(planeID, moduleID);
    if (moduleCal == nullptr) return 0;  // return numeric_limits<double>::quiet_NaN()
    return moduleCal->GetInterceptFullSpc();
}
 
std::set<int> TRestDataSetGainMap::GetPlaneIDs() const {
    std::set<int> planeIDs;
    for (const auto& mc : fModulesCal) planeIDs.insert(mc.GetPlaneId());
    return planeIDs;
}
 
std::set<int> TRestDataSetGainMap::GetModuleIDs(const int planeId) const {
    std::set<int> moduleIDs;
    for (const auto& mc : fModulesCal)
        if (mc.GetPlaneId() == planeId) moduleIDs.insert(mc.GetModuleId());
    return moduleIDs;
}
 
std::map<int, std::set<int>> TRestDataSetGainMap::GetModuleIDs() const {
    std::map<int, std::set<int>> moduleIds;
    for (const int planeId : GetPlaneIDs())
        moduleIds.insert(std::pair<int, std::set<int>>(planeId, GetModuleIDs(planeId)));
    return moduleIds;
}
 
TRestDataSetGainMap& TRestDataSetGainMap::operator=(TRestDataSetGainMap& src) {
    SetName(src.GetName());
    fOutputFileName = src.GetOutputFileName();
    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));
    return *this;
}
 
void TRestDataSetGainMap::SetModule(const Module& moduleCal) {
    for (auto& i : fModulesCal) {
        if (i.GetPlaneId() == moduleCal.GetPlaneId() && i.GetModuleId() == moduleCal.GetModuleId()) {
            i = moduleCal;
            return;
        }
    }
    fModulesCal.push_back(moduleCal);
}
 
void TRestDataSetGainMap::Import(const std::string& fileName) {
    if (fileName.empty()) {
        RESTError << "No input calibration file defined" << RESTendl;
        return;
    }
 
    if (TRestTools::isRootFile(fileName)) {
        RESTInfo << "Opening " << fileName << RESTendl;
        TFile* f = TFile::Open(fileName.c_str(), "READ");
        if (f == nullptr) {
            RESTError << "Cannot open calibration file " << fileName << RESTendl;
            return;
        }
 
        TRestDataSetGainMap* cal = nullptr;
        if (f != nullptr) {
            TIter nextkey(f->GetListOfKeys());
            TKey* key;
            while ((key = (TKey*)nextkey())) {
                std::string kName = key->GetClassName();
                if (REST_Reflection::GetClassQuick(kName.c_str()) != nullptr &&
                    REST_Reflection::GetClassQuick(kName.c_str())->InheritsFrom("TRestDataSetGainMap")) {
                    cal = f->Get<TRestDataSetGainMap>(key->GetName());
                    *this = *cal;
                }
            }
        }
    } else
        RESTError << "File extension not supported for " << fileName << RESTendl;
    if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) PrintMetadata();
}
 
void TRestDataSetGainMap::Export(const std::string& fileName) {
    if (!fileName.empty()) fOutputFileName = fileName;
    if (fOutputFileName.empty()) {
        RESTError << "No output file defined" << RESTendl;
        return;
    }
 
    if (TRestTools::GetFileNameExtension(fOutputFileName) == "root") {
        TFile* f = TFile::Open(fOutputFileName.c_str(), "UPDATE");
        this->Write(GetName());
        f->Close();
        delete f;
        RESTInfo << "Calibration saved to " << fOutputFileName << RESTendl;
    } else
        RESTError << "File extension for " << fOutputFileName << "is not supported." << RESTendl;
}
 
void TRestDataSetGainMap::PrintMetadata() {
    TRestMetadata::PrintMetadata();
    RESTMetadata << " Calibration dataset: " << fCalibFileName << RESTendl;
    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;
    RESTMetadata << " Spatial observable X: " << fSpatialObservableX << RESTendl;
    RESTMetadata << " Spatial observable Y: " << fSpatialObservableY << RESTendl;
    RESTMetadata << "-----------------------------------------------" << RESTendl;
    for (auto& i : fModulesCal) i.Print();
    RESTMetadata << "***********************************************" << RESTendl;
    RESTMetadata << RESTendl;
}
 
std::pair<int, int> TRestDataSetGainMap::Module::GetIndexMatrix(const double x, const double y) const {
    int index_x = -1, index_y = -1;
 
    if (fSplitX.upper_bound(x) != fSplitX.end()) {
        index_x = std::distance(fSplitX.begin(), fSplitX.upper_bound(x)) - 1;
        if (index_x < 0) {
            RESTWarning << "index_x < 0 for x = " << x << " and fSplitX[0]=" << *fSplitX.begin()
                        << p->RESTendl;
            index_x = 0;
        }
    } else {
        RESTWarning << "x is out of split for x = " << x << p->RESTendl;
        index_x = fSplitX.size() - 2;
    }
 
    if (fSplitY.upper_bound(y) != fSplitY.end()) {
        index_y = std::distance(fSplitY.begin(), fSplitY.upper_bound(y)) - 1;
        if (index_y < 0) {
            RESTWarning << "index_y < 0 for y = " << y << " and fSplitY[0]=" << *fSplitY.begin()
                        << p->RESTendl;
            index_y = 0;
        }
    } else {
        RESTWarning << "y is out of split for y = " << y << p->RESTendl;
        index_y = fSplitY.size() - 2;
    }
 
    return std::make_pair(index_x, index_y);
}
double TRestDataSetGainMap::Module::GetSlope(const double x, const double y) const {
    auto [index_x, index_y] = GetIndexMatrix(x, y);
    if (fSlope.empty()) {
        RESTError << "Calibration slope matrix is empty. Returning 0" << p->RESTendl;
        return 0;
    }
 
    if (index_x > (int)fSlope.size() || index_y > (int)fSlope.at(0).size()) {
        RESTError << "Index out of range. Returning 0" << p->RESTendl;
        return 0;
    }
 
    return fSlope[index_x][index_y];
}
 
double TRestDataSetGainMap::Module::GetIntercept(const double x, const double y) const {
    auto [index_x, index_y] = GetIndexMatrix(x, y);
    if (fIntercept.empty()) {
        RESTError << "Calibration constant matrix is empty. Returning 0" << p->RESTendl;
        return 0;
    }
 
    if (index_x > (int)fIntercept.size() || index_y > (int)fIntercept.at(0).size()) {
        RESTError << "Index out of range. Returning 0" << p->RESTendl;
        return 0;
    }
 
    return fIntercept[index_x][index_y];
}
 
void TRestDataSetGainMap::Module::SetSplits() {
    SetSplitX();
    SetSplitY();
}
void TRestDataSetGainMap::Module::SetSplitX() {
    if (fNumberOfSegmentsX < 1) {
        RESTError << "SetSplitX: fNumberOfSegmentsX must be >=1." << p->RESTendl;
        return;
    }
    std::set<double> split;
    for (int i = 0; i <= fNumberOfSegmentsX; i++) {  // <= so that the last segment is included
        double x =
            fReadoutRange.X() + ((fReadoutRange.Y() - fReadoutRange.X()) / (float)fNumberOfSegmentsX) * i;
        split.insert(x);
    }
    SetSplitX(split);
}
 
void TRestDataSetGainMap::Module::SetSplitX(const std::set<double>& splitX) {
    if (splitX.size() < 2) {
        RESTError << "SetSplitX: split size must be >=2 (start and end of range must be included)."
                  << p->RESTendl;
        return;
    }
    if (!fSlope.empty())
        RESTWarning << "SetSplitX: changing split but current gain map and calibration paremeters correspond "
                       "to previous splitting. Use GenerateGainMap() to update them."
                    << p->RESTendl;
    fSplitX = splitX;
    fNumberOfSegmentsX = fSplitX.size() - 1;
}
void TRestDataSetGainMap::Module::SetSplitY() {
    if (fNumberOfSegmentsY < 1) {
        RESTError << "SetSplitY: fNumberOfSegmentsY must be >=1." << p->RESTendl;
        return;
    }
    std::set<double> split;
    for (int i = 0; i <= fNumberOfSegmentsY; i++) {  // <= so that the last segment is included
        double y =
            fReadoutRange.X() + ((fReadoutRange.Y() - fReadoutRange.X()) / (float)fNumberOfSegmentsY) * i;
        split.insert(y);
    }
    SetSplitY(split);
}
 
void TRestDataSetGainMap::Module::SetSplitY(const std::set<double>& splitY) {
    if (splitY.size() < 2) {
        RESTError << "SetSplitY: split size must be >=2 (start and end of range must be included)."
                  << p->RESTendl;
        return;
    }
    if (!fSlope.empty())
        RESTWarning << "SetSplitY: changing split but current gain map and calibration paremeters correspond "
                       "to previous splitting. Use GenerateGainMap() to update them."
                    << p->RESTendl;
    fSplitY = splitY;
    fNumberOfSegmentsY = fSplitY.size() - 1;
}
 
void TRestDataSetGainMap::Module::GenerateGainMap() {
    //--- Initial checks and settings ---
    std::string dsFileName = fDataSetFileName;
    if (dsFileName.empty()) dsFileName = p->GetCalibrationFileName();
    if (dsFileName.empty()) {
        RESTError << "No calibration file defined" << p->RESTendl;
        return;
    }
 
    if (!TRestTools::fileExists(dsFileName)) {
        RESTError << "Calibration file " << dsFileName << " does not exist." << p->RESTendl;
        return;
    }
    if (!TRestTools::isDataSet(dsFileName)) RESTWarning << dsFileName << " is not a dataset." << p->RESTendl;
    TRestDataSet dataSet;
    dataSet.EnableMultiThreading(true);
    dataSet.Import(dsFileName);
    fDataSetFileName = dsFileName;
 
    dataSet.SetDataFrame(dataSet.MakeCut(p->GetCut()));
 
    if (fSplitX.empty()) SetSplitX();
    if (fSplitY.empty()) SetSplitY();
 
    //--- Get the calibration range if not provided (default is 0,0) ---
    if (fCalibRange.X() >= fCalibRange.Y()) {
        // Get spectrum for this file
        std::string cut = fDefinitionCut;
        if (cut.empty()) cut = "1";
        auto histo = dataSet.GetDataFrame().Filter(cut).Histo1D({"temp", "", fNBins, 0, 0}, GetObservable());
        std::unique_ptr<TH1F> hpunt = std::unique_ptr<TH1F>(static_cast<TH1F*>(histo->Clone()));
        // double xMin = hpunt->GetXaxis()->GetXmin();
        double xMax = hpunt->GetXaxis()->GetXmax();
 
        // Reduce the range to avoid the possible empty (nCounts<1%) end part of the spectrum
        double fraction = 1, nAtEndSpc = 0, step = 0.66;
        while (nAtEndSpc * 1. / hpunt->Integral() < 0.01 && fraction > 0.001) {
            fraction *= step;
            nAtEndSpc = hpunt->Integral(hpunt->FindFixBin(hpunt->GetXaxis()->GetXmax() * fraction),
                                        hpunt->FindFixBin(hpunt->GetXaxis()->GetXmax()));
        }
        xMax = hpunt->GetXaxis()->GetXmax() * fraction /
               step;  // previous step is the last one that nAtEndSpc<1%
        // Set the calibration range if needed
        // fCalibRange.SetX(xMin);
        fCalibRange.SetY(xMax);
        hpunt.reset();  // delete hpunt;
        RESTDebug << "Calibration range (auto)set to (" << fCalibRange.X() << "," << fCalibRange.Y() << ")"
                  << p->RESTendl;
    }
 
    // --- 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
    std::string cut = fDefinitionCut;
    if (cut.empty()) cut = "1";
    auto histoMod = dataSet.GetDataFrame().Filter(cut).Histo1D(
        {"tempMod", "", fNBins, fCalibRange.X(), fCalibRange.Y()}, GetObservable());
    std::unique_ptr<TH1F> hpuntMod = std::unique_ptr<TH1F>(static_cast<TH1F*>(histoMod->Clone()));
    fFullSpectrum->Add(hpuntMod.get());
 
    //--- Definition of histogram matrix ---
    std::vector<std::vector<TH1F*>> h(fNumberOfSegmentsX, std::vector<TH1F*>(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);
            h[i][j] = new TH1F(name.c_str(), "", fNBins, fCalibRange.X(),
                               fCalibRange.Y());  // h[column][row] equivalent to h[x][y]
        }
    }
 
    // build the spectrum for each segment
    auto itX = fSplitX.begin();
    for (size_t i = 0; i < h.size(); i++) {
        auto itY = fSplitY.begin();
        for (size_t j = 0; j < h.at(0).size(); j++) {
            // Get the segment limits from the splits
            auto xLower = *itX;
            auto xUpper = *std::next(itX);
            auto yLower = *itY;
            auto yUpper = *std::next(itY);
 
            std::string segment_cut = "";
            if (!GetSpatialObservableX().empty())
                segment_cut += GetSpatialObservableX() + ">=" + std::to_string(xLower) + "&&" +
                               GetSpatialObservableX() + "<" + std::to_string(xUpper);
            if (!GetSpatialObservableY().empty())
                segment_cut += "&&" + GetSpatialObservableY() + ">=" + std::to_string(yLower) + "&&" +
                               GetSpatialObservableY() + "<" + std::to_string(yUpper);
            if (!fDefinitionCut.empty()) segment_cut += "&&" + fDefinitionCut;
            if (segment_cut.empty()) segment_cut = "1";
            RESTExtreme << "Segment[" << i << "][" << j << "] cut: " << segment_cut << p->RESTendl;
            auto histo = dataSet.GetDataFrame()
                             .Filter(segment_cut)
                             .Histo1D({"temp", "", h[i][j]->GetNbinsX(), h[i][j]->GetXaxis()->GetXmin(),
                                       h[i][j]->GetXaxis()->GetXmax()},
                                      GetObservable());
            std::unique_ptr<TH1F> hpunt = std::unique_ptr<TH1F>(static_cast<TH1F*>(histo->Clone()));
            h[i][j]->Add(hpunt.get());
            hpunt.reset();  // delete hpunt;
            itY++;
        }
        itX++;
    }
 
    //--- Fit every peak energy for every segment ---
    std::vector<std::vector<double>> calParSlope(fNumberOfSegmentsX,
                                                 std::vector<double>(fNumberOfSegmentsY, 0));
    std::vector<std::vector<double>> calParIntercept(fNumberOfSegmentsX,
                                                     std::vector<double>(fNumberOfSegmentsY, 0));
    fSegLinearFit = std::vector(h.size(), std::vector<TGraph*>(h.at(0).size(), nullptr));
    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;
 
    //--- Fit every peak energy for the whole module ---
    delete fFullLinearFit;
    fFullLinearFit = new TGraph();
    auto [intercept, slope] = FitPeaks(fFullSpectrum, fFullLinearFit);
    fFullSlope = slope;
    fFullIntercept = intercept;
}
 
std::pair<double, double> 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);
    }
    std::shared_ptr<TGraph> graph = std::shared_ptr<TGraph>(new TGraph());
    RESTExtreme << "Fitting peaks for " << hSeg->GetName() << p->RESTendl;
 
    // Search for peaks --> peakPos
    std::unique_ptr<TSpectrum> s(new TSpectrum(2 * fEnergyPeaks.size() + 1));
    std::vector<double> 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<double>());
    const double ratio =
        peakPos.size() == 0 ? 1 : peakPos.front() / fEnergyPeaks.front();  // to estimate peak position
 
    // Initialize graph for linear fit
    graph->SetName("grFit");
    graph->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();
        }
 
        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);
                            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);
                    }
                }
            }
 
            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
        graph->SetPoint(c++, mu, energy);
    }
    s.reset();  // delete s;
 
    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;
    }
 
    // Linear fit
    std::unique_ptr<TF1> linearFit;
    linearFit = std::unique_ptr<TF1>(new TF1("linearFit", "pol1"));
    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));
}
 
void TRestDataSetGainMap::Module::Refit(const TVector2& position, const double energyPeak,
                                        const TVector2& range) {
    auto [index_x, index_y] = GetIndexMatrix(position.X(), position.Y());
    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;
    }
    Refit((size_t)index_x, (size_t)index_y, (size_t)peakNumber, range);
}
 
void TRestDataSetGainMap::Module::Refit(const size_t x, const size_t y, const size_t peakNumber,
                                        const TVector2& range) {
    if (fSegSpectra.empty()) {
        RESTError << "No gain map found. Use GenerateGainMap() first." << p->RESTendl;
        return;
    }
    if (x >= fSegSpectra.size() || y >= fSegSpectra.at(0).size()) {
        RESTError << "Segment with index (" << x << ", " << y << ") not found" << 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());
    TH1F* h = fSegSpectra.at(x).at(y);
    while (h->GetFunction(name.c_str()))  // clear previous fits for this peakNumber
        h->GetListOfFunctions()->Remove(h->GetFunction(name.c_str()));
    h->Fit(g, "R+Q0");  // use 0 to not draw the fit but save it
 
    // Change the point of the graph
    UpdateCalibrationFits(x, y);
}
 
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);
}
 
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();
}
 
void TRestDataSetGainMap::Module::UpdateCalibrationFits(const size_t x, const size_t y) {
    if (fSegSpectra.empty()) {
        RESTError << "No gain map found. Use GenerateGainMap() first." << p->RESTendl;
        return;
    }
    if (x >= fSegSpectra.size() || y >= fSegSpectra.at(0).size()) {
        RESTError << "Segment with index (" << x << ", " << y << ") not found" << p->RESTendl;
        return;
    }
 
    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<double, double> 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
    int c = 0;
    for (size_t i = 0; i < fEnergyPeaks.size(); i++) {
        std::string fitName = (std::string) "g" + std::to_string(i);
        TF1* g = h->GetFunction(fitName.c_str());
        if (!g) {
            RESTWarning << "No fit ( " << fitName << " ) found for energy peak " << fEnergyPeaks[i]
                        << " in histogram " << h->GetName() << 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);
    }
 
    // 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
 
    return std::make_pair(lf->GetParameter(0), lf->GetParameter(1));
}
 
void TRestDataSetGainMap::Module::UpdateCalibrationFitsFullSpc() {
    auto [intercept, slope] = UpdateCalibrationFits(fFullSpectrum, fFullLinearFit);
    fFullSlope = slope;
    fFullIntercept = intercept;
}
 
void TRestDataSetGainMap::Module::LoadConfigFromTiXmlElement(const TiXmlElement* module) {
    if (module == nullptr) {
        RESTError << "TRestDataSetGainMap::Module::LoadConfigFromTiXmlElement: module is nullptr"
                  << p->RESTendl;
        return;
    }
 
    std::string el = !module->Attribute("planeId") ? "Not defined" : module->Attribute("planeId");
    if (!(el.empty() || el == "Not defined")) this->SetPlaneId(StringToInteger(el));
    el = !module->Attribute("moduleId") ? "Not defined" : module->Attribute("moduleId");
    if (!(el.empty() || el == "Not defined")) this->SetModuleId(StringToInteger(el));
 
    el = !module->Attribute("moduleDefinitionCut") ? "Not defined" : module->Attribute("moduleDefinitionCut");
    if (!(el.empty() || el == "Not defined")) this->SetModuleDefinitionCut(el);
 
    el = !module->Attribute("numberOfSegmentsX") ? "Not defined" : module->Attribute("numberOfSegmentsX");
    if (!(el.empty() || el == "Not defined")) this->SetNumberOfSegmentsX(StringToInteger(el));
    el = !module->Attribute("numberOfSegmentsY") ? "Not defined" : module->Attribute("numberOfSegmentsY");
    if (!(el.empty() || el == "Not defined")) this->SetNumberOfSegmentsY(StringToInteger(el));
    el = !module->Attribute("readoutRange") ? "Not defined" : module->Attribute("readoutRange");
    if (!(el.empty() || el == "Not defined")) this->SetReadoutRange(StringTo2DVector(el));
 
    el = !module->Attribute("calibRange") ? "Not defined" : module->Attribute("calibRange");
    if (!(el.empty() || el == "Not defined")) this->SetCalibrationRange(StringTo2DVector(el));
    el = !module->Attribute("nBins") ? "Not defined" : module->Attribute("nBins");
    if (!(el.empty() || el == "Not defined")) this->SetNBins(StringToInteger(el));
 
    el = !module->Attribute("dataSetFileName") ? "Not defined" : module->Attribute("dataSetFileName");
    if (!(el.empty() || el == "Not defined")) this->SetDataSetFileName(el);
 
    el = !module->Attribute("zeroPoint") ? "Not defined" : module->Attribute("zeroPoint");
    if (!(el.empty() || el == "Not defined")) this->SetZeroPoint(ToLower(el) == "true");
    el = !module->Attribute("autoRangePeaks") ? "Not defined" : module->Attribute("autoRangePeaks");
    if (!(el.empty() || el == "Not defined")) this->SetAutoRangePeaks(ToLower(el) == "true");
 
    // Get peaks energy and range
    TiXmlElement* peakDefinition = (TiXmlElement*)module->FirstChildElement("peak");
    while (peakDefinition != nullptr) {
        double energy = 0;
        TVector2 range = TVector2(0, 0);
 
        std::string ell =
            !peakDefinition->Attribute("energy") ? "Not defined" : peakDefinition->Attribute("energy");
        if (ell.empty() || ell == "Not defined") {
            RESTError << "< peak variable key does not contain energy!" << p->RESTendl;
            exit(1);
        }
        energy = StringToDouble(ell);
 
        ell = !peakDefinition->Attribute("range") ? "Not defined" : peakDefinition->Attribute("range");
        if (!(ell.empty() || ell == "Not defined")) range = StringTo2DVector(ell);
 
        this->AddPeak(energy, range);
        peakDefinition = (TiXmlElement*)peakDefinition->NextSiblingElement();
    }
}
 
void TRestDataSetGainMap::Module::DrawSpectrum(const TVector2& position, bool drawFits, int color,
                                               TCanvas* c) {
    std::pair<size_t, size_t> index = GetIndexMatrix(position.X(), position.Y());
    DrawSpectrum(index.first, index.second, drawFits, color, c);
}
 
void TRestDataSetGainMap::Module::DrawSpectrum(const int index_x, const int index_y, bool drawFits, int color,
                                               TCanvas* c) {
    if (fSegSpectra.size() == 0) {
        RESTError << "Spectra matrix is empty." << p->RESTendl;
        return;
    }
    if (index_x < 0 || index_y < 0 || index_x >= (int)fSegSpectra.size() ||
        index_y >= (int)fSegSpectra.at(index_x).size()) {
        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) + "_" +
                        std::to_string(index_x) + "_" + std::to_string(index_y);
        c = new TCanvas(t.c_str(), t.c_str());
    }
 
    auto xLower = *std::next(fSplitX.begin(), index_x);
    auto xUpper = *std::next(fSplitX.begin(), index_x + 1);
    auto yLower = *std::next(fSplitY.begin(), index_y);
    auto yUpper = *std::next(fSplitY.begin(), index_y + 1);
    std::string tH = "Spectrum x=[" + DoubleToString(xLower, "%g") + ", " + DoubleToString(xUpper, "%g") +
                     ") y=[" + DoubleToString(yLower, "%g") + ", " + DoubleToString(yUpper, "%g") + ");" +
                     GetObservable() + ";counts";
    fSegSpectra[index_x][index_y]->SetTitle(tH.c_str());
 
    if (color > 0) fSegSpectra[index_x][index_y]->SetLineColor(color);
    size_t colorT = fSegSpectra[index_x][index_y]->GetLineColor();
    fSegSpectra[index_x][index_y]->Draw("same");
 
    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) 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::DrawSpectrum(const bool drawFits, const int color, TCanvas* c) {
    if (fSegSpectra.size() == 0) {
        RESTError << "Spectra matrix is empty." << p->RESTendl;
        return;
    }
    if (!c) {
        std::string t = "spectrum_" + 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 != fSegSpectra.size() * fSegSpectra.at(0).size()) {
        RESTError << "Canvas " << c->GetName() << " has " << nPads << " pads, but "
                  << fSegSpectra.size() * fSegSpectra.at(0).size() << " are needed." << p->RESTendl;
        return;
    } else if (nPads == 0)
        c->Divide(fSegSpectra.size(), fSegSpectra.at(0).size());
 
    for (size_t i = 0; i < fSegSpectra.size(); i++) {
        for (size_t j = 0; j < fSegSpectra[i].size(); j++) {
            int pad = fSegSpectra.size() * (fSegSpectra[i].size() - 1) + 1 + i - fSegSpectra.size() * j;
            c->cd(pad);
            DrawSpectrum(i, j, drawFits, color, c);
        }
    }
}
void TRestDataSetGainMap::Module::DrawFullSpectrum(const bool drawFits, const int color, TCanvas* c) {
    if (!fFullSpectrum) {
        RESTError << "Spectrum is empty." << p->RESTendl;
        return;
    }
 
    if (!c) {
        std::string t = "fullSpc_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId);
        c = new TCanvas(t.c_str(), t.c_str());
    }
    c->cd();
 
    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) {
    std::pair<size_t, size_t> index = GetIndexMatrix(position.X(), position.Y());
    DrawLinearFit(index.first, index.second, c);
}
 
void TRestDataSetGainMap::Module::DrawLinearFit(const int index_x, const int index_y, TCanvas* c) {
    if (fSegLinearFit.size() == 0) {
        RESTError << "Spectra matrix is empty." << p->RESTendl;
        return;
    }
    if (index_x < 0 || index_y < 0 || index_x >= (int)fSegLinearFit.size() ||
        index_y >= (int)fSegLinearFit.at(index_x).size()) {
        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);
        c = new TCanvas(t.c_str(), t.c_str());
    }
    auto xLower = *std::next(fSplitX.begin(), index_x);
    auto xUpper = *std::next(fSplitX.begin(), index_x + 1);
    auto yLower = *std::next(fSplitY.begin(), index_y);
    auto yUpper = *std::next(fSplitY.begin(), index_y + 1);
    std::string tH = "Linear Fit x=[" + DoubleToString(xLower, "%g") + ", " + DoubleToString(xUpper, "%g") +
                     ") y=[" + DoubleToString(yLower, "%g") + ", " + DoubleToString(yUpper, "%g") + ");" +
                     GetObservable() + ";energy";
    fSegLinearFit[index_x][index_y]->SetTitle(tH.c_str());
    fSegLinearFit[index_x][index_y]->Draw("AL*");
}
 
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++) {
            int pad = fSegLinearFit.size() * (fSegLinearFit[i].size() - 1) + 1 + i - fSegLinearFit.size() * j;
            c->cd(pad);
            DrawLinearFit(i, j, c);
        }
    }
}
 
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;
        return;
    }
    if (fSegLinearFit.size() == 0) {
        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") + ";" +
                        GetSpatialObservableX() + ";" + GetSpatialObservableY();  // + " keV";
    std::string t = "gainMap" + std::to_string(peakNumber) + "_" + std::to_string(fPlaneId) + "_" +
                    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(), fNumberOfSegmentsX, fReadoutRange.X(),
                              fReadoutRange.Y(), fNumberOfSegmentsY, fReadoutRange.X(), fReadoutRange.Y());
 
    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++) {
        auto itY = fSplitY.begin();
        for (size_t j = 0; j < fSegLinearFit.at(0).size(); j++) {
            auto xLower = *itX;
            auto xUpper = *std::next(itX);
            auto yLower = *itY;
            auto yUpper = *std::next(itY);
            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++;
        }
        itX++;
    }
    hGainMap->SetStats(0);
    hGainMap->Draw("colz");
    hGainMap->SetBarOffset(0.2);
    hGainMap->Draw("TEXT SAME");
}
 
void TRestDataSetGainMap::Module::Print() const {
    RESTMetadata << "-----------------------------------------------" << p->RESTendl;
    RESTMetadata << " Plane ID: " << fPlaneId << p->RESTendl;
    RESTMetadata << " Module ID: " << fModuleId << p->RESTendl;
    RESTMetadata << " Definition cut: " << fDefinitionCut << p->RESTendl;
    RESTMetadata << p->RESTendl;
 
    RESTMetadata << " Calibration dataset: " << fDataSetFileName << p->RESTendl;
    RESTMetadata << p->RESTendl;
 
    RESTMetadata << " Energy peaks: ";
    for (const auto& peak : fEnergyPeaks) RESTMetadata << peak << ", ";
    RESTMetadata << p->RESTendl;
    bool anyRange = false;
    for (const auto& r : fRangePeaks)
        if (r.X() < r.Y()) {
            RESTMetadata << " Range peaks: ";
            anyRange = true;
            break;
        }
    if (anyRange)
        for (const auto& r : fRangePeaks) RESTMetadata << "(" << r.X() << ", " << r.Y() << ") ";
    if (anyRange) RESTMetadata << p->RESTendl;
    RESTMetadata << " Auto range peaks: " << (fAutoRangePeaks ? "true" : "false") << p->RESTendl;
    RESTMetadata << " Zero point: " << (fZeroPoint ? "true" : "false") << p->RESTendl;
    RESTMetadata << " Calibration range: (" << fCalibRange.X() << ", " << fCalibRange.Y() << " )"
                 << p->RESTendl;
    RESTMetadata << " Number of bins: " << fNBins << p->RESTendl;
    RESTMetadata << p->RESTendl;
 
    RESTMetadata << " Number of segments X: " << fNumberOfSegmentsX << p->RESTendl;
    RESTMetadata << " Number of segments Y: " << fNumberOfSegmentsY << p->RESTendl;
    // RESTMetadata << " Draw: " << (fDrawVar ? "true" : "false") << p->RESTendl;
    RESTMetadata << " Readout range (" << fReadoutRange.X() << ", " << fReadoutRange.Y() << " )"
                 << p->RESTendl;
    RESTMetadata << "SplitX: ";
    for (auto& i : fSplitX) {
        RESTMetadata << "    " << i;
    }
    RESTMetadata << p->RESTendl;
    RESTMetadata << "SplitY: ";
    for (auto& i : fSplitY) {
        RESTMetadata << "    " << i;
    }
    RESTMetadata << p->RESTendl;
    RESTMetadata << p->RESTendl;
 
    RESTMetadata << " Slope: " << p->RESTendl;
    size_t maxSize = 0;
    for (auto& x : fSlope)
        if (maxSize < x.size()) maxSize = x.size();
    for (size_t j = 0; j < maxSize; j++) {
        for (size_t k = 0; k < fSlope.size(); k++) {
            if (j < fSlope[k].size())
                RESTMetadata << DoubleToString(fSlope[k][fSlope[k].size() - 1 - j], "%.3e") << "   ";
            else
                RESTMetadata << "         ";
        }
        RESTMetadata << p->RESTendl;
    }
    RESTMetadata << " Intercept: " << p->RESTendl;
    maxSize = 0;
    for (auto& x : fIntercept)
        if (maxSize < x.size()) maxSize = x.size();
    for (size_t j = 0; j < maxSize; j++) {
        for (size_t k = 0; k < fIntercept.size(); k++) {
            if (j < fIntercept[k].size())
                RESTMetadata << DoubleToString(fIntercept[k][fIntercept[k].size() - 1 - j], "%+.3e") << "   ";
            else
                RESTMetadata << "         ";
        }
        RESTMetadata << p->RESTendl;
    }
    RESTMetadata << p->RESTendl;
    RESTMetadata << " Full slope: " << DoubleToString(fFullSlope, "%.3e") << p->RESTendl;
    RESTMetadata << " Full intercept: " << DoubleToString(fFullIntercept, "%+.3e") << p->RESTendl;
 
    RESTMetadata << "-----------------------------------------------" << p->RESTendl;
}