Calibrator.hpp

Go to the documentation of this file.

 
#ifndef H_CALIBRATOR
#define H_CALIBRATOR
 
 
#include "DetectorEvent.hpp"
#include "ContextBase.hpp"
#include "StageDetectorEvent.hpp"
 
#include <fmt/format.h>
 
#include <cstdlib>
#include <cstdio>
#include <random>
#include <array>
#include <string>
#include <cstdint>
#include <cassert>
#include <fstream>
#include <stdexcept>
#include <iostream>
#include <bitset>
 
 
 
namespace pottu {
 
    class Calibrator : public StageDetectorEvent {
    public:
 
    struct item_t {
        float    c0{0};
        float    c1{1};
        float    c2{0};
        uint32_t time_correction{0};
    };
 
    Calibrator()
        : _ctxh( ContextBase::getActive().createHandle("Calibrator") )
        {}
 
    void readFromGainsFile( const std::string &filename ) {
        std::fstream fsgains;
        fsgains.open( filename, std::fstream::in );
        if( !fsgains.is_open() ) {
        _ctxh->logError( fmt::format("read_from_gains_file(): Couldn't open file \'{}\'!", filename ) );
        throw std::runtime_error( "File could not be opened" );
        }
 
        _ctxh->logInfo( fmt::format("Opening gains file (grain compatible) {}.", filename ) );
        
        unsigned int row = 0;
            unsigned int channel;
        double       cc[3];
        int          timecorr;
        int          tmp;
        uint32_t     numImported{0};
        std::string  s_line, s_tmp;
        std::bitset<4096>  foundChannels;
        while( !fsgains.eof() ) {
        ++row;
        std::getline( fsgains, s_line );
        tmp = sscanf( s_line.c_str(), " %d = %lf %lf %lf %d", &channel, cc, cc+1, cc+2, &timecorr );
        if( tmp == 5 ) {
            if( channel >= _items.size() ) {
            _ctxh->logError( fmt::format(
                         "Calibrator::readFromGainsFile():\n"
                         "  Too big channel number!\n"
                         "  Error on gains file \'{}\' on line {}!\n"
                         "  Line: \'{}\'"
                         , filename, row, s_line ) );
            throw std::runtime_error( "Too big channel number. Max is 4095." );
            }
            if( foundChannels[channel] ) {
            _ctxh->logWarning( fmt::format("Overwriting existing gains for channel {}.", channel ) );
            } else {
            ++numImported;
            foundChannels.set(channel);
            }
            _items[channel].c0 = cc[0];
            _items[channel].c1 = cc[1];
            _items[channel].c2 = cc[2];
            _items[channel].time_correction = timecorr;
        }
        if( tmp != 5 && tmp != 0 && tmp != EOF ) {
            _ctxh->logError( fmt::format(
                    "Calibrator::readFromGainsFile():\n"
                    "  Error on gains file \'{}\' on line {}!\n"
                    "  Line \'{}\'",
                    filename, row, s_line ) );
            throw std::runtime_error( "Syntax error when reading gains file" );
        }
        }
        _ctxh->logInfo( fmt::format("Set gains for {} unique channels.", numImported ) );
    }
 
 
    inline double getCalibratedValue( unsigned int ch, int data ) noexcept {
        const double d_data = double(data) + _udist( _re );
        const item_t &item = _items[ch];
        return item.c0 + item.c1*d_data + item.c2*d_data*d_data;
    }
    
    inline void fillCalibratedEnergy( DetectorEvent &event ) noexcept {
        event.e = getCalibratedValue( event.ch, event.adc );
    }
 
 
 
    virtual void process( std::vector<DetectorEvent> &data ) noexcept {
        auto proci = _ctxh->processInstance();
 
        for( auto &event : data ) {
        const double d_data = double(event.adc) + _udist( _re );
        const item_t &item = _items[event.ch];
        event.e = item.c0 + item.c1*d_data + item.c2*d_data*d_data;
        }
    }
 
    
    void setCoefficients( unsigned int ch, float c0, float c1, float c2 ) noexcept {
        if( ch > 4095 )
        exit(EXIT_FAILURE);
        _items[ch].c0 = c0;
        _items[ch].c1 = c1;
        _items[ch].c2 = c2;
        _items[ch].time_correction = 0;
    }
 
 
    private:
    ContextHandle                   *_ctxh;
    
    std::array<item_t,4096>          _items;
    std::minstd_rand                 _re;
    std::uniform_real_distribution<> _udist {-0.5, 0.5};
    };
 
    
}
 
 
#endif