source: trunk/FACT++/src/DataCalib.cc@ 11812

#include "DataCalib.h"

#include "FAD.h"
#include "FitsFile.h"
#include "DimDescriptionService.h"

#include "externals/fits.h"

using namespace std;

int  DataCalib::fStep = 0;
bool DataCalib::fProcessing = false;

vector<int32_t> DataCalib::fOffset(1440*1024, 0);
vector<int64_t> DataCalib::fGain  (1440*1024, 4096);
vector<int64_t> DataCalib::fTrgOff(1440*1024, 0);
vector<float>   DataCalib::fStats (1440*1024*6+3);

uint64_t DataCalib::fNumOffset = 1;
uint64_t DataCalib::fNumGain   = 2000;
uint64_t DataCalib::fNumTrgOff = 1;

void DataCalib::Restart()
{
    // Default gain:
    // 0.575*[45590]*2.5V / 2^16 = 0.99999 V
    fOffset.assign(1440*1024, 0);
    fGain.assign  (1440*1024, 4096);
    fTrgOff.assign(1440*1024, 0);

    fNumOffset = 1;
    fNumGain   = 2000;
    fNumTrgOff = 1;

    reinterpret_cast<uint32_t*>(fStats.data())[0] = 0;
    reinterpret_cast<uint32_t*>(fStats.data())[1] = 0;
    reinterpret_cast<uint32_t*>(fStats.data())[2] = 0;

    int i=0;
    while (i<1024*1440*2+3)  // Set mean and RMS to 0
        fStats[i++] = 0;
    while (i<1024*1440*3+3)
        fStats[i++] = 2000./4096; // Set mean to 0.5
    while (i<1024*1440*6+3)
        fStats[i++] = 0;   // Set everything else to 0

    fStep = 0;
    fProcessing = false;
}

void DataCalib::Update(DimDescribedService &dim)
{
    dim.Update(fStats);
}

bool DataCalib::Open(RUN_HEAD* h)
{
    if (h->NPix != 1440)
    {
        fMsg.Error("Number of pixels in header not 1440.");
        return false;
    }

    if (h->Nroi != 1024)
    {
        fMsg.Error("Region of interest not 1024.");
        return false;
    }

    if (fProcessing)
    {
        fMsg.Warn("Previous run not yet finished!");
        return false;
    }

    if (fStep==4)
    {
        fMsg.Warn("DRS Calibration already finished... please restart!");
        return false;
    }

    fProcessing = true;

    ostringstream name;
    name << "drs-calib-" << fStep;
    fFileName = name.str();

    Reset();
    InitSize(1440, 1024);

    return true;
}

bool DataCalib::WriteEvt(EVENT *e)
{
    // FIXME: SET StartPix to 0 if StartPix is -1

    if (fStep==0)
    {
        AddRel(e->Adc_Data, e->StartPix);
    }
    if (fStep==1)
    {
        AddRel(e->Adc_Data, e->StartPix, fOffset.data(), fNumOffset);
    }
    if (fStep==2)
    {
        AddAbs(e->Adc_Data, e->StartPix, fOffset.data(), fNumOffset);
    }

    return true;
}

bool DataCalib::ReadFits(const string &str)
{
    try
    {
        fits file(str);
        if (!file)
            return false;

        if (file.GetStr("TELESCOP")!="FACT")
        {
            Error("Reading "+str+" failed: Not a valid FACT file (TELESCOP not FACT in header)");
            return false;
        }

        if (!file.HasKey("STEP"))
        {
            Error("Reading "+str+" failed: Is not a DRS calib file (STEP not found in header)");
            return false;
        }

        if (file.GetNumRows()!=1)
        {
            Error("Reading "+str+" failed: Number of rows in table is not 1.");
            return false;
        }

        vector<float> data(1440*1024*6+3);

        float *base = data.data();

        file.SetPtrAddress("RunNumberBaseline",      base,   sizeof(float));
        file.SetPtrAddress("RunNumberGain",          base+1, sizeof(float));
        file.SetPtrAddress("RunNumberTriggerOffset", base+2, sizeof(float));
        file.SetPtrAddress("BaselineMean",           base+0*1024*1440+3, sizeof(float)*1024*1440);
        file.SetPtrAddress("BaselineRms",            base+1*1024*1440+3, sizeof(float)*1024*1440);
        file.SetPtrAddress("GainMean",               base+2*1024*1440+3, sizeof(float)*1024*1440);
        file.SetPtrAddress("GainRms",                base+3*1024*1440+3, sizeof(float)*1024*1440);
        file.SetPtrAddress("TriggerOffsetMean",      base+4*1024*1440+3, sizeof(float)*1024*1440);
        file.SetPtrAddress("TriggerOffsetRms",       base+5*1024*1440+3, sizeof(float)*1024*1440);

        if (!file.GetNextRow())
        {
            Error("Reading data from "+str+" failed.");
            return false;
        }

        fStep = file.GetUInt("STEP");

        memcpy(fStats.data(), base, fStats.size()*sizeof(float));

        return true;
    }
    catch (const runtime_error &e)
    {
        Error("Exception reading "+str+": "+e.what());
        return false;
    }
}

void DataCalib::WriteFits()
{
#ifdef HAVE_FITS
    FitsFile file(fMsg);

    file.AddColumn('I', "RunNumberBaseline");
    file.AddColumn('I', "RunNumberGain");
    file.AddColumn('I', "RunNumberTriggerOffset");
    file.AddColumn('F', "BaselineMean", 1024*1440);
    file.AddColumn('F', "BaselineRms", 1024*1440);
    file.AddColumn('F', "GainMean", 1024*1440);
    file.AddColumn('F', "GainRms", 1024*1440);
    file.AddColumn('F', "TriggerOffsetMean", 1024*1440);
    file.AddColumn('F', "TriggerOffsetRms", 1024*1440);

    fFileName = FormFileName(GetRunId(), "drs.fits");

    if (!file.OpenFile(fFileName))
        return;

    if (!file.OpenTable("DrsCalibration"))
        return;

    if (!file.WriteDefaultKeys("fadctrl"))
        return;

    if (!file.WriteKeyNT("STEP", fStep, ""))
        return;

    vector<char> buf;
    buf.reserve(fStats.size()*sizeof(float));

    char *src  = reinterpret_cast<char*>(fStats.data());
    char *end  = reinterpret_cast<char*>(fStats.data()+1024*1440*6+3);
    char *dest = buf.data();

    while (src<end)
    {
        reverse_copy(src, src+sizeof(float), dest);
        src  += sizeof(float);
        dest += sizeof(float);
    }

    if (!file.AddRow())
        return;

    if (!file.WriteData(buf.data(), 1024*1440*sizeof(float)*6+3))
        return;

    ostringstream str;
    str << "Wrote DRS calibration data (step=" << fStep << ") to '" << fFileName << "'";
    Info(str.str());
#endif
}

bool DataCalib::Close(RUN_TAIL *)
{
    if (fStep==0)
    {
        fOffset.assign(fSum.begin(), fSum.end());
        fNumOffset = fNumEntries;

        // Scale ADC data from 12bit to 2000mV
        GetSampleStats(fStats.data()+3, 2000./4096);
        reinterpret_cast<uint32_t*>(fStats.data())[0] = GetRunId();;
    }
    if (fStep==1)
    {
        fGain.assign(fSum.begin(), fSum.end());
        fNumGain = fNumEntries*fNumOffset;

        // DAC:  0..2.5V == 0..65535
        // V-mV: 1000
        //fNumGain *= 2500*50000;
        //for (int i=0; i<1024*1440; i++)
        //    fGain[i] *= 65536;
        fNumGain *= 1953125;
        for (int i=0; i<1024*1440; i++)
            fGain[i] *= 1024;

        // Scale ADC data from 12bit to 2000mV
        GetSampleStats(fStats.data()+1024*1440*2+3, 2000./4096/fNumOffset);//0.5);
        reinterpret_cast<uint32_t*>(fStats.data())[1] = GetRunId();;
    }
    if (fStep==2)
    {
        fTrgOff.assign(fSum.begin(), fSum.end());
        fNumTrgOff = fNumEntries*fNumOffset;

        // Scale ADC data from 12bit to 2000mV
        GetSampleStats(fStats.data()+1024*1440*4+3, 2000./4096/fNumOffset);//0.5);
        reinterpret_cast<uint32_t*>(fStats.data())[2] = GetRunId();;
    }

    if (fStep>=0 && fStep<=2)
        WriteFits();

    fProcessing = false;

    fDim.Update(fStats);

    fStep++;

    return true;
}

void DataCalib::Apply(int16_t *val, const int16_t *start, uint32_t roi)
{
    CalibData::Apply(val, start, roi,
                     fOffset.data(), fNumOffset,
                     fGain.data(),   fNumGain,
                     fTrgOff.data(), fNumTrgOff);
}

void DataCalib::Apply(float *vec, int16_t *val, const int16_t *start, uint32_t roi)
{
    CalibData::Apply(vec, val, start, roi,
                     fOffset.data(), fNumOffset,
                     fGain.data(),   fNumGain,
                     fTrgOff.data(), fNumTrgOff);
}