source: fact/tools/pyscripts/sandbox/kraehenb/calfits.h@ 13428

//********************************
//
// Class Calfits
// Wrapper class for fits.h or 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 pyfits.h++
// root [2] .L ../sandbox/kraehenb/calfits.h++
//
// Usage in Python:
// See CalFitsTest.py
//
//********************************

//ToDo: alle Event-Parameter zugänglich, shared library creation debuggen

#ifndef CALFITS_H
#define CALFITS_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

#ifdef __EXCEPTIONS
#include <stdexcept>
#endif

#include "../../pyfact/pyfits.h"

class CalFits
{
public:
        //No standard constructor CalFits()!
        
        //Direct handlers of the files
        fits datafile, calibfile; //Class name should be PyFits or better FactPyFits...
        
        //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;
        
        //Calibration variables
        float* calib_baselinemean;
        float* calib_gainmean;
        float* calib_triggeroffsetmean;
        //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;
        
        CalFits(const string &datafilename, const string &calibfilename) //Constructor with two filenames
                : datafile(datafilename),
                        calibfile(calibfilename),
                        npcaldata(NULL)
        {
                //Read basic parameters of the two files
//              std::cout << "...Reading basic file parameters..." << std::endl;
                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;
                
//              std::cout << "Column sizes: " << calib_blm_size << " " << calib_gm_size << " " << calib_tom_size << std::endl;
                
                //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.";
#ifdef __EXCEPTIONS
                        throw runtime_error(str.str());
#else
                        gLog << ___err___ << "ERROR - " << str.str() << endl;
                        return;
#endif
                }
                nevents = datafile.GetNumRows();
                
                //Read the calibration data
//              std::cout << "...Reading calibration data..." << std::endl;
                calib_baselinemean = new float[calibfile.GetN("BaselineMean")];
                calibfile.SetPtrAddress("BaselineMean", calib_baselinemean, calibfile.GetN("BaselineMean"));
                calib_gainmean = new float[calibfile.GetN("GainMean")];
                calibfile.SetPtrAddress("GainMean", calib_gainmean, calibfile.GetN("GainMean"));
                calib_triggeroffsetmean = new float[calibfile.GetN("TriggerOffsetMean")];
                calibfile.SetPtrAddress("TriggerOffsetMean", calib_triggeroffsetmean, calibfile.GetN("TriggerOffsetMean"));
                calibfile.GetRow(0);
                
                //Set the event pointers
//              std::cout << "...Setting event pointers..." << std::endl;
                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"));
        }
        
        ~CalFits()
        {
                delete[] calib_baselinemean;
                delete[] calib_gainmean;
                delete[] calib_triggeroffsetmean;
                delete[] event_data;
                delete[] event_offset;
                delete[] event_boardtimes;
        }
        
        bool GetCalEvent()
        {
                if(!npcaldata) {
                        ostringstream str;
                        str << "Pointer to the calibrated data not initialized!";
#ifdef __EXCEPTIONS
                        throw runtime_error(str.str());
#else
                        gLog << ___err___ << "ERROR - " << str.str() << endl;
                        return false;
#endif
                }
                if(datafile.GetNextRow() == false) {
//                      std::cout << "Last event reached..." << std::endl;
                        return false;
                }
                else {
                        UInt_t drs_calib_offset;
                        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])%calib_blm_size;
                                        npcaldata[pixel*data_nroi+slice] = double((event_data[pixel*data_nroi+slice]*2000./4096.-calib_baselinemean[pixel*calib_blm_size+drs_calib_offset]-calib_triggeroffsetmean[pixel*data_nroi+slice])/calib_gainmean[pixel*calib_blm_size+drs_calib_offset]*1907.35);
                                        //Note: data_nroi=calib_nroi, calib_blm_size=calib_gm_size
                                }
                        }
                }
                return true;
        }
        
        void SetNpcaldataPtr(double *ptr)
        {
                npcaldata = ptr;
                return;
        }
};
#endif