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

//********************************
//
// 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:
// ...
//
// Usage in Python:
// ...handles etc...
//
//********************************

//ToDo: Code aufräumen, ToDo's abarbeiten, alle Event-Parameter zugänglich, shared library creation debuggen

#include <cstdio>
#include <vector>
#include <string>
#include "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...
        
        //Calibration variables
        UInt_t calib_nroi, calib_npix;
        vector<float> calib_basemean;
        vector<float> calib_gainmean;
        vector<float> calib_triggeroffsetmean;
        //ToDo: use arrays of size 1440x1024 (x2 for wrap-arounds) and read all variables into those
        
        //Data variables
        UInt_t data_nroi, data_npix, data_ndata;
        
        //Common variables
        UInt_t nroi, npix, nevents;
        
        //Event variables
        UInt_t data_eventid;
        vector<int16_t> data_data;
        vector<int16_t> data_offset;
        
        //Vector for calibrated event
        vector<float> calevent;
        
        CalFits(const string &datafilename, const string &calibfilename) //Construct with two filenames
                : datafile(datafilename),
                        calibfile(calibfilename)
        {
                //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");
                nevents = datafile.GetNumRows();
                data_ndata = datafile.GetN("Data");
                
                //ToDo: check nroi*npix = ndata, calib_nroi=data_nroi, calib_npix=data_npix
                
                //Read the calibration data
                std::cout << "...Reading calibration data..." << std::endl;
                calib_basemean.resize(calibfile.GetN("BaselineMean"),0);
                calibfile.SetVecAddress("BaselineMean", calib_basemean);
                calib_gainmean.resize(calibfile.GetN("GainMean"),0);
                calibfile.SetVecAddress("GainMean", calib_gainmean);
                calib_triggeroffsetmean.resize(calibfile.GetN("TriggerOffsetMean"),0);
                calibfile.SetVecAddress("TriggerOffsetMean", calib_triggeroffsetmean);
                calibfile.GetRow(0);
                
                //Set the data vectors
                std::cout << "...Setting datafile pointers..." << std::endl;
                datafile.SetRefAddress("EventNum", data_eventid);
                data_data.resize(data_ndata,0);
                datafile.SetVecAddress("Data", data_data);
                calevent.resize(data_ndata,0);
                data_offset.resize(datafile.GetN("StartCellData"),0);
                datafile.SetVecAddress("StartCellData", data_offset);
        }
        
        //double* GetCalEvent() { }
        int GetCalEvent()
        {
                if(datafile.GetNextRow() == false) {
                        std::cout << "Last event reached..." << std::endl;
                        //ToDo: raise exception or stop or whatever GetNextRow does
                }
                else {
                        UInt_t drs_calib_offset;
                        for(UInt_t pixel=0;pixel<data_npix;pixel++) {
                                for(UInt_t slice=0;slice<data_nroi;slice++) {
                                        //ToDo: replace 1024 by column size/npix
                                        drs_calib_offset = (slice+data_offset[pixel])%1024;
                                        calevent[pixel*data_nroi+slice] = float((data_data[pixel*data_nroi+slice]*2000./4096.-calib_basemean[pixel*1024+drs_calib_offset]-calib_triggeroffsetmean[pixel*calib_nroi+slice])/calib_gainmean[pixel*1024+drs_calib_offset]*1907.35);
                                        //Note: data_nroi=calib_nroi
                                }
                        }
                }
                return 0;
        }
};