Jungfraujoch/jungfrau/JFModuleGainCalibration.cpp

// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only

#include <fstream>
#include <cmath>

#include "JFModuleGainCalibration.h"
#include "../common/JFJochException.h"

JFModuleGainCalibration::JFModuleGainCalibration() : gain(RAW_MODULE_SIZE * 4) {
    for (int i = 0; i < RAW_MODULE_SIZE; i++) {
        gain[i + offset_g0] = DEFAULT_G0_FACTOR;
        gain[i + offset_g1] = DEFAULT_G1_FACTOR;
        gain[i + offset_g2] = DEFAULT_G2_FACTOR;
        gain[i + offset_hg0] = DEFAULT_HG0_FACTOR;
    }
}

JFModuleGainCalibration::JFModuleGainCalibration(const std::string &filename) : gain(RAW_MODULE_SIZE * 4) {
    std::ifstream file(filename.c_str(), std::fstream::binary);
    if (!file.is_open())
        throw JFJochException(JFJochExceptionCategory::GainFileOpenError, "Gain file cannot be opened");
    try {
        file.read((char *) gain.data(), gain.size() * sizeof(double));
    } catch (...) {
        throw JFJochException(JFJochExceptionCategory::GainFileOpenError, "Gain file cannot be read");
    }
}

JFModuleGainCalibration::JFModuleGainCalibration(const std::vector<double> &vec) {
    if (vec.size() != 3 * RAW_MODULE_SIZE) {
        throw JFJochException(JFJochExceptionCategory::GainFileOpenError,
                              "Wrong size of input vector for gain calibration");
    }
    gain = vec;
}

const std::vector<double> &JFModuleGainCalibration::GetGainCalibration() const {
    return gain;
}

void JFModuleGainCalibration::ExportG0(DeviceOutput *output) const {
    for (int i = 0; i < RAW_MODULE_SIZE; i++)
        ((float *) output->pixels)[i] = static_cast<float>(1.0 / gain[offset_g0 + i]);
}

void JFModuleGainCalibration::ExportFixedG1(DeviceOutput *output) const {
    for (int i = 0; i < RAW_MODULE_SIZE; i++)
        ((float *) output->pixels)[i] = static_cast<float>(fixed_g1_gain_coeff / gain[offset_g1 + i]);
}

void JFModuleGainCalibration::ExportG1(DeviceOutput *output) const {
    for (int i = 0; i < RAW_MODULE_SIZE; i++)
        ((float *) output->pixels)[i] = static_cast<float>(1.0 / gain[offset_g1 + i]);
}

void JFModuleGainCalibration::ExportG2(DeviceOutput *output) const {
    for (int i = 0; i < RAW_MODULE_SIZE; i++)
        ((float *) output->pixels)[i] = static_cast<float>(1.0 / gain[offset_g2 + i]);
}

void JFModuleGainCalibration::ExportHG0(DeviceOutput *output) const {
    for (int i = 0; i < RAW_MODULE_SIZE; i++)
        ((float *) output->pixels)[i] = static_cast<float>(1.0 / gain[offset_hg0 + i]);
}

double JFModuleGainCalibration::GetMean(size_t offset) const{
    double ret = 0;
    for (int i = 0; i < RAW_MODULE_SIZE; i++)
        ret += gain.at(offset + i);
    return ret / static_cast<double>(RAW_MODULE_SIZE);
}

double JFModuleGainCalibration::GetStdDev(size_t offset) const {
    double mean = GetMean(offset);

    double ret = 0;
    for (int i = 0; i < RAW_MODULE_SIZE; i++) {
        double val = gain.at(offset + i);
        ret += (val - mean) * (val - mean) / (RAW_MODULE_SIZE - 1) ;
    }
    return sqrt(ret);
}

double JFModuleGainCalibration::GetG0Mean() const {
    return GetMean(offset_g0);
}

double JFModuleGainCalibration::GetG1Mean() const {
    return GetMean(offset_g1);
}

double JFModuleGainCalibration::GetG2Mean() const {
    return GetMean(offset_g2);
}

double JFModuleGainCalibration::GetHG0Mean() const {
    return GetMean(offset_hg0);
}

double JFModuleGainCalibration::GetG0StdDev() const {
    return GetStdDev(offset_g0);
}

double JFModuleGainCalibration::GetG1StdDev() const {
    return GetStdDev(offset_g1);
}

double JFModuleGainCalibration::GetG2StdDev() const {
    return GetStdDev(offset_g2);
}

double JFModuleGainCalibration::GetHG0StdDev() const {
    return GetStdDev(offset_hg0);
}

void JFGainCalibration::LoadGain(const std::vector<std::string>& filenames) {
    c.clear();
    for (const auto& filename : filenames)
        c.emplace_back(filename);
}

const std::vector<JFModuleGainCalibration>& JFGainCalibration::GetCalibration() const {
    return c;
}