source: trunk/Mars/mcore/factofits.h@ 17226

/*
 * factofits.h
 *
 *  Created on: Oct 16, 2013
 *      Author: lyard
 */

#ifndef FACTOFITS_H_
#define FACTOFITS_H_

#include "zofits.h"

#ifndef __MARS__
namespace std
{
#else
using namespace std;
#endif

class factofits : public zofits
{

    public:

        factofits(uint32_t numTiles=1000,
                  uint32_t rowPerTile=100,
                  uint64_t maxUsableMem=0) : zofits(numTiles, rowPerTile, maxUsableMem)
        {
            fStartCellsOffset = -1;
            fDataOffset       = -1;
        }
        factofits(const char* fname,
                  uint32_t numTiles=1000,
                  uint32_t rowPerTile=100,
                  uint64_t maxUsableMem=0) : zofits(fname, numTiles, rowPerTile, maxUsableMem)
        {
            fStartCellsOffset = -1;
            fDataOffset       = -1;
        }
        virtual ~factofits()
        {
        }

        //whether or not a calibration was given to the file writer
        virtual bool IsOffsetCalibration()
        {
            return (fOffsetCalibration.size() != 0);
        }

        //assign a given drs offset calibration
        void SetDrsCalibration(const vector<float> &calib)
        {
            if (!IsOffsetCalibration())
                fOffsetCalibration.resize(1440*1024);

            for (uint32_t i=0; i<1440*1024; i++)
                fOffsetCalibration[i] = (int16_t)(calib[i]*4096/2000);
        }

        void SetDrsCalibration(const vector<float>& calib)
        {
            if (calib.size() != 1440*1024)
    #ifdef __EXCEPTIONS
            throw runtime_error("Cannot load calibration with anything else than 1024 samples per pixel");
    #else
            gLog << ___err___ << "ERROR - Cannot load calibration with anything else than 1024 samples per pixel");
    #endif
            SetDrsCalibration(calib.data());
        }

        void SetDrsCalibration(const vector<int16_t>& vec)
        {
            if (!IsOffsetCalibration())
                fOffsetCalibration.resize(1440*1024);

            for (uint32_t i=0; i<1440*1024; i++)
                fOffsetCalibration[i] = vec[i];
        }

        void SetDrsCalibration(const DrsCalibration& drs)
        {
            if (drs.fNumOffset==0)
                return;

            if (!IsOffsetCalibration())
                fOffsetCalibration.resize(1440*1024);

            for (uint32_t i=0; i<1024*1440; i++)
                fOffsetCalibration[i] = drs.fOffset[i]/drs.fNumOffset;
        }

        bool WriteTableHeader(const char* name="DATA")
        {
            if (!zofits::WriteTableHeader(name))
                return false;

            if (IsOffsetCalibration())
            {//retrieve the column storing the start cell offsets, if required.

                for (auto it=fRealColumns.begin(); it!=fRealColumns.end(); it++)//Table.cols.begin(); it!= fTable.cols.end(); it++)
                {
                    if (it->col.name == "StartCellData")
                        fStartCellsOffset = it->col.offset;
                    if (it->col.name == "Data")
                    {
                        fNumSlices = it->col.num;
                        fDataOffset = it->col.offset;
                        if (fNumSlices % 1440 != 0)
                        {
#ifdef __EXCEPTIONS
                            throw runtime_error("Number of data samples not a multiple of 1440.");
#else
                            gLog << ___err___ << "ERROR - Number of data samples not a multiple of 1440. Doing it uncalibrated." << endl;
#endif
                            fOffsetCalibration.resize(0);
                        }
                        fNumSlices /= 1440;
                    }
                }
                if (fStartCellsOffset < 0)
                {
#ifdef __EXCEPTIONS
                    throw runtime_error("FACT Calibration requested, but \"StartCellData\" column not found.");
#else
                    gLog << ___err___ << "ERROR - FACT Calibration requested, but \"StartCellData\" column not found. Doing it uncalibrated." << endl;
#endif
                    //throw away the calibration data
                    fOffsetCalibration.resize(0);
                }
                if (fDataOffset < 0)
                {
#ifdef __EXCEPTIONS
                    throw runtime_error("FACT Calibration requested, but \"Data\" column not found.");
#else
                    gLog << ___err___ << "ERROR - FACT Calibration requested, but \"Data\" column not found. Doing it uncalibrated." << endl;
#endif
                    //throw away the calibration data
                    fOffsetCalibration.resize(0);
                }
            }

            return true;
        }

        virtual bool WriteDrsOffsetsTable()
        {
            if (!IsOffsetCalibration())
                return false;

            ofits c;
            c.SetStr("XTENSION", "BINTABLE"            , "binary table extension");
            c.SetInt("BITPIX"  , 8                     , "8-bit bytes");
            c.SetInt("NAXIS"   , 2                     , "2-dimensional binary table");
            c.SetInt("NAXIS1"  , 1024*1440*2           , "width of table in bytes");
            c.SetInt("NAXIS2"  , 1                     , "number of rows in table");
            c.SetInt("PCOUNT"  , 0                     , "size of special data area");
            c.SetInt("GCOUNT"  , 1                     , "one data group (required keyword)");
            c.SetInt("TFIELDS" , 1                     , "number of fields in each row");
            c.SetStr("CHECKSUM", "0000000000000000"    , "Checksum for the whole HDU");
            c.SetStr("DATASUM" ,  "         0"         , "Checksum for the data block");
            c.SetStr("EXTNAME" , "ZDrsCellOffsets"     , "name of this binary table extension");
            c.SetStr("TTYPE1"  , "OffsetCalibration"   , "label for field   1");
            c.SetStr("TFORM1"  , "1474560I"            , "data format of field: 2-byte INTEGER");
            c.End();

            vector<char> swappedOffsets;
            swappedOffsets.resize(1024*1440*sizeof(int16_t));
            revcpy<sizeof(int16_t)>(swappedOffsets.data(), (char*)(fOffsetCalibration.data()), 1024*1440);

            Checksum datasum;
            datasum.add(swappedOffsets.data(), sizeof(int16_t)*1024*1440);

            ostringstream dataSumStr;
            dataSumStr << datasum.val();
            c.SetStr("DATASUM", dataSumStr.str());

            datasum += c.WriteHeader(*this);

            const off_t here_I_am = tellp();

            c.SetStr("CHECKSUM", datasum.str());
            c.WriteHeader(*this);

            seekp(here_I_am);

            write(swappedOffsets.data(), swappedOffsets.size());

            AlignTo2880Bytes();

            return good();
        }

        virtual void DrsOffsetCalibrate(char* target_location)
        {
            if (IsOffsetCalibration())
            {
                int16_t* startCell = reinterpret_cast<int16_t*>(target_location + fStartCellsOffset);
                int16_t* data      = reinterpret_cast<int16_t*>(target_location + fDataOffset);

                for (uint32_t ch=0; ch<1440; ch++)
                {
                    if (startCell[ch] < 0)
                    {
                        data += fNumSlices;
                        continue;
                    }

                    const int16_t modStart = startCell[ch]%1024;
                    const int16_t *off     = fOffsetCalibration.data() + ch*1024;

                    const int16_t* cal        = off+modStart;
                    const int16_t* end_stride = data+fNumSlices;

                    if (modStart+fNumSlices > 1024)
                    {
                        while (cal < off+1024)
                            *data++ -= *cal++;
                        cal = off;
                    }

                    while (data<end_stride)
                        *data++ -= *cal++;
                }
            }
        }

private:
        //Offsets calibration stuff.
        vector<int16_t> fOffsetCalibration; ///< The calibration itself
        int32_t         fStartCellsOffset;  ///< Offset in bytes for the startcell data
        int32_t         fDataOffset;        ///< Offset in bytes for the data
        int32_t         fNumSlices;         ///< Number of samples per pixel per event

}; //class factofits

#ifndef __MARS
}; //namespace std
#endif

#endif /* FACTOFITS_H_ */