source: fact/tools/pyscripts/pyfact/next_steps/calfactfits_new.h@ 17689

//********************************
//
// Class CalFactFits
// Wrapper class for factfits.h (or fits.h, pyfits.h)
// Provides fast access to calibrated events of FACT raw files
//
// written by Thomas Kraehenbuehl, ETH Zurich
// tpk@phys.ethz.ch, +41 44 633 3973
// April 2012
//
//********************************
//
// Compilation (root in pyfact directory)
// root [0] gSystem->Load("/usr/lib64/libz.so");
// root [1] .L izstream.h++
// root [2] .L factfits.h++
// root [3] .L calfactfits.h++
// 
// Try /usr/lib/x86_64-linux-gnu/libz.so at ETH for the libz.so file.
//
// Note: the produced izstream_h.so must be in LD_LIBRARY_PATH, e.g.
// export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/home/guest/kraehenb/software/PythonIncludes/
//
// Usage in Python:
// See pyscripts/examples/CalFitsTest.py
//
//********************************

//ToDo: shared library creation debuggen

#ifndef CALFACTFITS_H
#define CALFACTFITS_H

#include <cstdio>
#include <string>

#ifndef __MARS__
#include <vector>
#include <iomanip>
#include <iostream>
#define gLog cerr
#define ___err___ ""
#define ___all___ ""
#else
#include "MLog.h"
#include "MLogManip.h"
#define ___err___ err
#define ___all___ all
#endif

#include <stdexcept>

#include "extern_Mars_mcore/factfits.h"

class CalFactFits
{
public:
    //No standard constructor CalFactFits()!

    //Direct handlers of the files
    factfits datafile, calibfile;

    //Basic file parameters
    UInt_t calib_nroi, calib_npix;
    UInt_t calib_blm_size, calib_gm_size, calib_tom_size;
    UInt_t data_nroi, data_npix, data_ndata;

    //Common variables
    UInt_t nroi, npix, nevents;

    // The same variables as above, but this time
    // they are doubled in memory, so the modulo operation is not necessary.
    // idea by TPK again
    // since we like to work with doubles for some reason, I don't know yet
    // I cast the float to double already in this stage.
    
    double* baseline;
    double* gain;
    double* trigger_offset;
    //Using <vector> instead of arrays makes no visible difference
    //ToDo: use arrays of size 1440x1024 (x2 for wrap-arounds) and read all variables into those

    //Event variables
    UInt_t event_id;
    UShort_t event_triggertype;
    int16_t* event_data;
    int16_t* event_offset;
    int32_t* event_boardtimes;
    double* npcaldata;

    // Switches, to control if or if not to omit the actual calibration and return 
    // the uncalibrated data instead of the calibrated
    bool do_calibration;          // default: true
    bool do_despiking;            // default: true
    double single_spike_limit;       // default: 25
    double double_spike_outer_limit; // default: 22
    double double_spike_inner_limit; // default: 4
    
    void SetDoCalibration( bool do_it){
        do_calibration = do_it;
    }
    void SetDoDespiking( bool do_it){
        do_despiking = do_it;
    }
    void SetSingleSpikeLimit( double limit ){
        single_spike_limit = limit;
    }
    void SetDoubleSpikeOuterLimit( double limit ){
        double_spike_outer_limit = limit;
    }
    void SetDoubleSpikeInnerLimit( double limit ){
        double_spike_inner_limit = limit;
    }



    CalFactFits(const string &datafilename, const string &calibfilename) //Constructor with two filenames
        : datafile(datafilename),
            calibfile(calibfilename),
            npcaldata(NULL)
    {
        calib_nroi = calibfile.GetUInt("NROI");
        calib_npix = calibfile.GetUInt("NPIX");
        data_nroi = datafile.GetUInt("NROI");
        data_npix = datafile.GetUInt("NPIX");
        data_ndata = datafile.GetN("Data");
        
        calib_blm_size = calibfile.GetN("BaselineMean")/calib_npix;
        calib_gm_size = calibfile.GetN("GainMean")/calib_npix;
        calib_tom_size = calibfile.GetN("TriggerOffsetMean")/calib_npix;
        
        do_calibration=true;
        do_despiking=true;
        single_spike_limit=25;
        double_spike_outer_limit=22;
        double_spike_inner_limit=4;
        
        //Define the common variables
        if((calib_nroi==data_nroi)&&(calib_npix==data_npix)&&(data_nroi*data_npix==data_ndata)&&(calib_blm_size==calib_gm_size)&&(calib_tom_size==calib_nroi)) {
            nroi = data_nroi;
            npix = data_npix;
        }
        else {
            ostringstream str;
            str << "Data/calib file error: NROI mismatch, NPIX mismatch, data column size wrong or calib columns mismatch.";
            throw runtime_error(str.str());
        }
        nevents = datafile.GetNumRows();
        
        //Read the calibration data
        //Calibration variables
        float* calib_baselinemean;
        float* calib_gainmean;
        float* calib_triggeroffsetmean;
        calib_baselinemean = new float[calibfile.GetN("BaselineMean")];
        calibfile.SetPtrAddress("BaselineMean", calib_baselinemean, calibfile.GetN("BaselineMean"));
        baseline = new double[calibfile.GetN("BaselineMean")*2];
            
        calib_gainmean = new float[calibfile.GetN("GainMean")];
        calibfile.SetPtrAddress("GainMean", calib_gainmean, calibfile.GetN("GainMean"));
        gain = new double[calibfile.GetN("GainMean")*2];
        
        calib_triggeroffsetmean = new float[calibfile.GetN("TriggerOffsetMean")];
        calibfile.SetPtrAddress("TriggerOffsetMean", calib_triggeroffsetmean, calibfile.GetN("TriggerOffsetMean"));
        trigger_offset = new double[calibfile.GetN("TriggerOffsetMean")*2];
        
        calibfile.GetRow(0);
        copy_calib_constants_twice( calib_baselinemean, baseline, 4096./2000., calib_npix, calib_blm_size);
        copy_calib_constants_twice( calib_gainmean, gain, 1./( 1907.35 /4096. *2000.), calib_npix, calib_gm_size);
        copy_calib_constants_twice( calib_triggeroffsetmean, trigger_offset, 4096. /2000., calib_npix, calib_tom_size);
        delete[] calib_baselinemean;
        delete[] calib_gainmean;
        delete[] calib_triggeroffsetmean;
        
        datafile.SetRefAddress("EventNum", event_id);
        datafile.SetRefAddress("TriggerType", event_triggertype);
        
        event_data = new int16_t[data_ndata];
        datafile.SetPtrAddress("Data", event_data, data_ndata);
        
        event_offset = new int16_t[datafile.GetN("StartCellData")];
        datafile.SetPtrAddress("StartCellData", event_offset, datafile.GetN("StartCellData"));
        
        event_boardtimes = new int32_t[datafile.GetN("BoardTime")];
        datafile.SetPtrAddress("BoardTime", event_boardtimes, datafile.GetN("BoardTime"));
    }

    void copy_calib_constants_twice( float *src, double *dest, double scale=1., int npix, int roi){
        for (int pix=0 ; pix < npix; ++pix) {
            for (int sl = 0; sl < roi; ++sl) {
                int orig_index = pix * roi + sl;
                int new_index1 = 2*pix * roi + sl;
                int new_index2 = (2*pix+1) * roi + sl;
                dest[new_index2] = dest[new_index1] = scale * src[orig_index];
            }
        }
    }

    ~CalFactFits() //Standard destructor
    {
        delete[] event_data;
        delete[] event_offset;
        delete[] event_boardtimes;
        delete[] baseline;
        delete[] gain;
        delete[] trigger_offset;
    }
   
    //Read calibrated event into the event variables
    bool GetCalEvent( omit_calibration=false, omit_despiking=false)
    {
        if(!npcaldata) {
            ostringstream str;
            str << "Pointer to the calibrated data not initialized!";
            throw runtime_error(str.str());
        }
        if(datafile.GetNextRow() == false) {
            return false;
        }
        else 
        {
            UInt_t drs_calib_offset;
            double raw, raw_bsl, shifted;
            if ( do_calibration) {
            for(UInt_t pixel = 0; pixel < data_npix; pixel++) {
                for(UInt_t slice = 0; slice < data_nroi; slice++) {
                    drs_calib_offset = slice + event_offset[pixel];
                    raw = (double)event_data[pixel*data_nroi+slice];
                    raw_bsl = raw - baseline[pixel*calib_blm_size*2+drs_calib_offset];
                    shifted = raw_bsl - trigger_offset[pixel*data_nroi*2+slice];
                    npcaldata[pixel*data_nroi+slice] = shifted / gain[pixel*calib_blm_size*2+drs_calib_offset];
                    //Note: data_nroi=calib_nroi, calib_blm_size=calib_gm_size
                }

                if ( do_despiking ){
                // second loop for despiking
                // stolen from TBs MARS ... DrsCalib.h or something like that
                double *p = npcaldata+(pixel*data_nroi);
                for (size_t i=1; i<data_nroi-2; i++) {
                    if (p[i]-p[i-1]>single_spike_limit && p[i]-p[i+1]>single_spike_limit) {
                        p[i] = (p[i-1]+p[i+1])/2;
                    }
                    if (p[i]-p[i-1]>double_spike_outer_limit 
                            && fabs(p[i]-p[i+1])<double_spike_inner_limit 
                            && p[i+1]-p[i+2]>double_spike_outer_limit) {
                        p[i] = (p[i-1]+p[i+2])/2;
                        p[i+1] = p[i];
                    }
                }
                } // end of if do_despiking
                
            }
            else { // not doing calibration
                for(UInt_t pixel = 0; pixel < data_npix; pixel++) {
                    for(UInt_t slice = 0; slice < data_nroi; slice++) {
                       npcaldata[pixel*data_nroi+slice] = (double)event_data[pixel*data_nroi+slice];
                    }
                }
            }
        }
        return true;
    }

    void SetNpcaldataPtr(double *numpyptr) //Set the pointer for the calibrated data to the numpy array
    {
        npcaldata = numpyptr;
        return;
    }
};
#endif /* CALFACTFITS_H */