Jungfraujoch/jungfrau/JFConversionFloatingPoint.cpp

// Copyright (2019-2022) Paul Scherrer Institute
// SPDX-License-Identifier: GPL-3.0-or-later

#include "JFConversionFloatingPoint.h"
#include "../common/RawToConvertedGeometry.h"

JFConversionFloatingPoint::JFConversionFloatingPoint()
: gain_g0(RAW_MODULE_SIZE), gain_g1(RAW_MODULE_SIZE), gain_g2(RAW_MODULE_SIZE),
pedestal_g0(RAW_MODULE_SIZE), pedestal_g1(RAW_MODULE_SIZE), pedestal_g2(RAW_MODULE_SIZE) {}

inline float one_over_gain_energy(double gain_factor, double energy) {
    double tmp = gain_factor * energy;
    if (!std::isfinite(tmp) || (tmp == 0.0))
        return std::numeric_limits<float>::infinity();
    else
        return static_cast<float>(1.0 / (gain_factor * energy));
}

void JFConversionFloatingPoint::Setup(const JFModuleGainCalibration &gain_calibration,
                                 const JFModulePedestal &in_pedestal_g0,
                                 const JFModulePedestal &in_pedestal_g1,
                                 const JFModulePedestal &in_pedestal_g2,
                                 double energy) {
    auto &gain_arr = gain_calibration.GetGainCalibration();
    memcpy(pedestal_g0.data(), in_pedestal_g0.GetPedestal(), RAW_MODULE_SIZE * sizeof(uint16_t));
    memcpy(pedestal_g1.data(), in_pedestal_g1.GetPedestal(), RAW_MODULE_SIZE * sizeof(uint16_t));
    memcpy(pedestal_g2.data(), in_pedestal_g2.GetPedestal(), RAW_MODULE_SIZE * sizeof(uint16_t));

    for (int i = 0; i < RAW_MODULE_SIZE; i++) {
        gain_g0[i] = one_over_gain_energy(gain_arr[i], energy);
        gain_g1[i] = one_over_gain_energy(gain_arr[i + RAW_MODULE_SIZE], energy);
        gain_g2[i] = one_over_gain_energy(gain_arr[i + 2 * RAW_MODULE_SIZE], energy);
    }
}

void JFConversionFloatingPoint::ConvertPacket(int16_t* dest, const uint16_t* source, uint16_t packet_number)  {
    for (int i = 0; i < 4; i++)
        ConvertLine(dest + i * RAW_MODULE_COLS, source  + i * RAW_MODULE_COLS, 4 * packet_number + i);
}

void JFConversionFloatingPoint::ConvertFP(float *dest, const uint16_t *source) {
    for (size_t i = 0; i < RAW_MODULE_SIZE; i++) {
        uint16_t gainbits = source[i] & 0xc000;
        uint16_t adc = source[i] & 0x3fff;
        int16_t pedestal_subtracted_adu;

        float expected = PIXEL_OUT_LOST;

        bool special_val = false;
        switch (gainbits) {
            case 0:
                pedestal_subtracted_adu = adc - pedestal_g0[i];
                expected = static_cast<float>(pedestal_subtracted_adu) * gain_g0[i];
                break;
            case 0x4000:
                pedestal_subtracted_adu = adc - pedestal_g1[i];
                expected = static_cast<float>(pedestal_subtracted_adu) * gain_g1[i];
                if (adc == 0) expected = PIXEL_OUT_G1_SATURATION, special_val = true;
                break;
            case 0xc000:
                pedestal_subtracted_adu = adc - pedestal_g2[i];
                expected = static_cast<float>(pedestal_subtracted_adu) * gain_g2[i];
                if (adc == 0)           expected = PIXEL_OUT_SATURATION, special_val = true;
                else if (adc == 0x3fff) expected = PIXEL_OUT_0xFFFF, special_val = true;
                break;
            default:
                expected = PIXEL_OUT_GAINBIT_2, special_val = true;;
                break;
        }

        if (std::isinf(expected))
            expected = PIXEL_OUT_LOST;
        else if ((expected > PIXEL_OUT_SATURATION) && !special_val)
            expected = PIXEL_OUT_SATURATION;
        dest[i] = expected;
    }
}

void JFConversionFloatingPoint::ConvertAdjustGeom(int16_t *dest, const uint16_t *source, int64_t slow_dir_step,
                                                  int64_t fast_dir_step, int64_t offset) {
    for (int line = 0; line < RAW_MODULE_LINES; line++) {
        int16_t tmp[RAW_MODULE_COLS];
        ConvertLine(tmp, source + line * RAW_MODULE_COLS, line);

        if ((line != 255) && (line != 256)) {
            LineCopyAndAdjustMultipixel<int16_t, int16_t>(dest + slow_dir_step * (line + ((line > 255) ? 2 : 0)),
                                                          tmp, INT16_MIN, INT16_MAX, fast_dir_step, offset);
        } else {
            LineCopyAndAdjustMultipixelMidRow<int16_t, int16_t>(dest + slow_dir_step * (line + 1), tmp,
                                                                INT16_MIN, INT16_MAX,
                                                                fast_dir_step, offset);
            LineCopyAndAdjustMultipixelMidRow<int16_t, int16_t>(dest + slow_dir_step * ((line > 255) ? 2 : 0), tmp,
                                                                INT16_MIN, INT16_MAX,
                                                                fast_dir_step, offset);
        }
    }
}

void JFConversionFloatingPoint::ConvertLine(int16_t *dest, const uint16_t *source, int line) {
    for (size_t i = 0; i < RAW_MODULE_COLS; i++) {
        uint16_t gainbits = source[i] & 0xc000;
        uint16_t adc = source[i] & 0x3fff;
        int16_t pedestal_subtracted_adu;

        float expected = PIXEL_OUT_LOST;

        switch (gainbits) {
            case 0:
                [[likely]]
                        pedestal_subtracted_adu = adc - pedestal_g0[i + line * RAW_MODULE_COLS];
                expected = static_cast<float>(pedestal_subtracted_adu) * gain_g0[i + line * RAW_MODULE_COLS];
                break;
            case 0x4000:
                pedestal_subtracted_adu = adc - pedestal_g1[i + line * RAW_MODULE_COLS];
                expected = static_cast<float>(pedestal_subtracted_adu) * gain_g1[i + line * RAW_MODULE_COLS];
                if (adc == 0) [[unlikely]] expected = PIXEL_OUT_G1_SATURATION;
                break;
            case 0xc000:
                pedestal_subtracted_adu = adc - pedestal_g2[i + line * RAW_MODULE_COLS];
                expected = static_cast<float>(pedestal_subtracted_adu) * gain_g2[i + line * RAW_MODULE_COLS];
                if (adc == 0) [[unlikely]] expected = PIXEL_OUT_SATURATION;
                else if (adc == 0x3fff) [[unlikely]] expected = PIXEL_OUT_0xFFFF;
                break;
            default:
                expected = PIXEL_OUT_GAINBIT_2;
                break;
        }

        if (std::isinf(expected) || (expected == INT16_MIN))
            dest[i] = PIXEL_OUT_LOST;
        else if (expected > INT16_MAX)
            dest[i] = INT16_MAX;
        else if (expected >= 0.0)
            [[likely]]
                    dest[i] = static_cast<uint16_t>(expected + 0.5f);
        else
            dest[i] = static_cast<uint16_t>(expected - 0.5f);
    }
}