source: fact/tools/rootmacros/DrsCalibration.C@ 13593

#include "DrsCalibration.h"

#include <iostream>
using namespace std;

float getValue( int slice, int pixel,
        vector<float> &drs_basemean,
        vector<float> &drs_gainmean,
        vector<float> &drs_triggeroffsetmean,
        unsigned int    RegionOfInterest,
        vector<int16_t> AllPixelDataVector,
        vector<int16_t> StartCellVector
){
        const float dconv = 2000/4096.0;

        float vraw, vcal;

        unsigned int pixel_pt;
        unsigned int slice_pt;
        unsigned int cal_pt;
        unsigned int drs_cal_offset;

        // printf("pixel = %d, slice = %d\n", slice, pixel);

        pixel_pt = pixel * RegionOfInterest;
        slice_pt = pixel_pt + slice;
        drs_cal_offset = ( slice + StartCellVector[ pixel ] )%RegionOfInterest;
        cal_pt    = pixel_pt + drs_cal_offset;

        vraw = AllPixelDataVector[ slice_pt ] * dconv;
        vcal = ( vraw - drs_basemean[ cal_pt ] - drs_triggeroffsetmean[ slice_pt ] ) / drs_gainmean[ cal_pt ]*1907.35;

        return( vcal );
}



// FACT raw data is stored in fits files in large vector<int16_t>. these
// containt the data of the full camera for a single event.
// in order to work with this data, one needs to extract the part, which
// belongs to the current pixel and apply a DRS calibration, based on a dedicated
// DRS calibration file.
//
// TODO: The calibration Data as well as the Data and StartCell vectors are a bit
// massive to give this function as parameters.
// better to have these vectors somehow bundled in classes like:
// RawData - class and
// DrsCalibration - class
// those classes would be the result of
// openDataFits() and
// openCalibFits()
//
size_t applyDrsCalibration( vector<float> &destination,
        int pixel, 
  int LeaveOutLeft, 
  int LeaveOutRight,
        vector<float> &drs_basemean,
        vector<float> &drs_gainmean,
        vector<float> &drs_triggeroffsetmean,
        unsigned int RegionOfInterest,
        vector<int16_t> AllPixelDataVector,
        vector<int16_t> StartCellVector,
        int verbosityLevel
){
        // In order to minimize mem free and alloc actions
        // the destination vector is only resized, if it is way too big
        // or too small
        size_t DestinationLength = RegionOfInterest-LeaveOutRight-LeaveOutLeft;
        if ( destination.size() < DestinationLength || destination.size() > 2 * DestinationLength ){
                destination.resize( DestinationLength );
        }

        //      We do not entirely know how the calibration constants, which are saved in a filename.drs.fits file 
        // were calculated, so it is not fully clear how they should be applied to the raw data for calibration.
        // apparently the calibration constants were transformed to the unit mV, which means we have to do the same to
        // the raw data prior to apply the calibration
        //
        // on the FAD board, there is a 12bit ADC, with a 2.0V range, so the factor between ADC units and mV is
        // ADC2mV = 2000/4096. = 0.48828125 (numerically exact)
        //
        // from the schematic of the FAD we learned, that the voltage at the ADC 
        // should be 1907.35 mV when the calibration DAC is set to 50000. 
        // 
        // One would further assume that the calibration constants are calculated like this:

        // The DRS Offset of each bin in each channel is the mean value in this very bin,
        // obtained from so called DRS pedestal data
        // Its value is about -1820 ADC units or -910mV

        // In order to obtain the DRS Gain of each bin of each channel
        // again data is takes, with the calibration DAC set to 50000
        // This is called DRS calibration data.
        // We assume the DRS Offset is already subtracted from the DRS calibration data
        // so the typical value is assumed to be ~3600 ADC units oder ~1800mV 
        // As mentioned before, the value *should* be 1907.35 mV 
        // So one might assume that the Gain is a number, which actually converts ~3600 ADC units into 1907.35mV for each bin of each channel.
        // So that the calibration procedure looks like this
        // TrueValue = (RawValue - Offset) * Gain
        // The numerical value of Gain would differ slightly from the theoretical value of 2000/4096.
        // But this is apparently not the case. 
        // The Gain, as it is stored in the DRS calibration file of FACT++ has numerical values 
        // around +1800.
        // So it seems one should calibrate like this:
        // TrueValue = (RawValue - Offset) / Gain * 1907.35

        // When these calibrations are done, one ends up with a quite nice calibrated voltage.
        // But it turns out that, if one returns the first measurement, and calculates the mean voltages 
        // in each bin of the now *logical* DRS pipeline, the mean voltage is not zero, but slightly varies 
        // So one can store these systematical deviations from zero in the logical pipeline as well, and subtract them.
        // The remaining question is, when to subtract them.
        // I assume, in the process of measuring this third calibration constant, the first two 
        // calibrations are already applied to the raw data. 

        // So the calculation of the calibrated volatage from some raw voltage works like this:
        // assume the raw voltage is the s'th sample in channel c. While the Trigger made the DRS stopp in its t'th cell.
        // note, that the DRS pipeline is always 1024 bins long. This is constant of the DRS4 chip. 

        // TrueValue[c][s] = ( RawValue[c][s] - Offset[c][ (c+t)%1024 ] ) / Gain[c][ (c+t)%1024 ] * 1907.35 - TriggerOffset[c][s] 
        
        const float dconv = 2000/4096.0;
        float vraw;

        if ( RegionOfInterest != AllPixelDataVector.size()/1440 ){
                if (verbosityLevel > 0){
                        cout << "RegionOfInterest != AllPixelDataVector.size()/1440" << endl;
                        cout << "RegionOfInterest: " << RegionOfInterest << endl;
                        cout << "AllPixelDataVector.size(): " << AllPixelDataVector.size() << endl;
                        cout << "aborting" << endl;
                }
                return 0;
        }

        // the vector drs_triggeroffsetmean is not 1440 * 1024 entries long
        // but has hopefully the length 1440 * RegionOfInterest (or longer)
        if ( drs_triggeroffsetmean.size() < 1440*RegionOfInterest ){
                if (verbosityLevel > 0){
                        cout << "Error: drs_triggeroffsetmean.size() < 1440*RegionOfInterest" << endl;
                        cout << "drs_triggeroffsetmean.size():" << drs_triggeroffsetmean.size() << endl;
                        cout << "RegionOfInterest" << RegionOfInterest << endl;
                        cout << "aborting" << endl;
                }
                return 0;
        }

        int DataPos, OffsetPos, TriggerOffsetPos;
        
        for ( unsigned int sl = LeaveOutLeft; sl < RegionOfInterest-LeaveOutRight ; sl++){

                DataPos = pixel * RegionOfInterest + sl;
                // Offset and Gain vector *should look the same
                OffsetPos = pixel * drs_basemean.size()/1440 + (sl + StartCellVector[pixel])%(drs_basemean.size()/1440);

                TriggerOffsetPos = pixel * drs_triggeroffsetmean.size()/1440 + sl;

                vraw = AllPixelDataVector[ DataPos ] * dconv;
                vraw -= drs_basemean[ OffsetPos ];
                vraw -= drs_triggeroffsetmean[ TriggerOffsetPos ];
                vraw /= drs_gainmean[ OffsetPos ];
                vraw *= 1907.35;

//              slice_pt = pixel_pt + sl;
//              drs_cal_offset = ( sl + StartCellVector[ pixel ] ) % RegionOfInterest;
//              cal_pt    = pixel_pt + drs_cal_offset;
//              vraw = AllPixelDataVector[ slice_pt ] * dconv;
//              vcal = ( vraw - drs_basemean[ cal_pt ] - drs_triggeroffsetmean[ slice_pt ] ) / drs_gainmean[ cal_pt ]*1907.35;
//              destination.push_back(vcal);

                destination[sl-LeaveOutLeft] = vraw;
        }

        return destination.size();
}