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

#include "DataCalib.h"

#include "FAD.h"
#include "FitsFile.h"
#include "DimDescriptionService.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, 2);
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   = 1;
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, 2);
    fTrgOff.assign(1440*1024, 0);

    fNumOffset = 1;
    fNumGain   = 1;
    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++] = 0.5; // 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;
}

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;

        // 0.5: Scale from ADC to Millivolt
        GetSampleStats(fStats.data()+3, 0.5);
        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;

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

        // 0.5: Scale from ADC to Millivolt
        GetSampleStats(fStats.data()+1024*1440*4+3, 1./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);
}