MApplyCalibration.cc

Go to the documentation of this file.

#include "MApplyCalibration.hh"
#include "QBaseType.hh"
#include "QEvent.hh"
#include "QRawEvent.hh"
#include "QMatrix.hh"
#include "QString.hh"
#include "TF1.h"
#include "QCalibrationParameters.hh"
#include "QCalibrationParametersHandle.hh"
 
REGISTER_MODULE(MApplyCalibration)
 
using std::string;
using std::cout;
using std::endl;
using namespace Diana;
 
MApplyCalibration::~MApplyCalibration()
{
    // called once in a sequence
}
 
void MApplyCalibration::Init(QEvent& ev)
{
    // called before event loop
    fInput = GetString("ParametersInput", "");
 
    // get dataset
    fDataset = GetInt("Dataset",-1,false);
    if(fDataset < 0) {
        GlobalHandle<QInt> dHandle("Dataset");
        GlobalData().Get("",&dHandle,"");
        fDataset = dHandle.Get();
    }
 
    fAmplitudeLabel = GetString("AmplitudeLabel","CorrectAmplitudes@Amplitude");
    ev.Add<QDouble>("Energy");
    ev.Add<QDouble>("EnergyUncertainty");
    fBlackList.clear();
 
    ev.Require<QPulseInfo>("DAQ","PulseInfo");
    ev.Require<QHeader>("DAQ","Header");
    ev.RequireByLabel<QDouble>(fAmplitudeLabel);
}
 
void MApplyCalibration::Do(QEvent& ev)
{
    // called event by event
    const QPulseInfo& pulseInfo = ev.Get<QPulseInfo>("DAQ","PulseInfo");
    const QHeader& header = ev.Get<QHeader>("DAQ","Header");
    int channel = pulseInfo.GetChannelId();
    int run = header.GetRun();
 
    const double amp = ev.GetByLabel<QDouble>(fAmplitudeLabel);
    double energy = -1;
 
    // get calibration parameters for this channel
    QCalibrationParametersHandle paramhandle(channel,fDataset);
    GlobalData().Get(GetString("CalCoefficientsOwner","CalCoefficients",false),&paramhandle,fInput);
 
    if (!paramhandle.IsValid()) {
        if(fBlackList.count(channel)) return;
        fBlackList.insert(channel);
        Error("Channel %d, will not be calibrated: %s",channel,paramhandle.GetError().GetDescription().c_str());
        return;
    }
 
    // check that the amplitude to be calibrated is the same on which the calibration parameters were calculated (stored in DB)
    // this check is performed only when the input of the calibration parameters is DB
    if (fInput == "DB" && fChannels.find(channel) == fChannels.end()) {
        fChannels.insert(channel);
        std::string sLabelstr;
        GlobalHandle<QString> sLabel("StabilizationLabel");
        sLabel.SetChannel(channel);
        sLabel.SetRun(run);
        GlobalData().Get(fAmplitudeLabel.GetOwner(),&sLabel,GetString("StabilizationLabelInput","CurrentReader",false));
        if (sLabel.IsValid()) {
            sLabelstr = sLabel.Get();
        }
        QGlobalLabel stabLabel = sLabelstr;;
 
        if (stabLabel.fOwner != paramhandle.GetStabLabel())
        {
            Error("The amplitude to be calibrated (%s) is not the same on which the calibration parameters were calculated (%s). Calibrating %s anyway.",stabLabel.fOwner.c_str(),paramhandle.GetStabLabel().c_str(),fAmplitudeLabel.GetStringLabel().c_str());
        }
 
        // writes a GlobalData (QString) with the name of the calibration module
        GlobalHandle<QString> clabel("CalibrationLabel");
        clabel.SetChannel(channel);
        clabel.SetDataset(fDataset);
        QGlobalLabel label;
        label.fOwner= paramhandle.GetLabel().fOwner;
        label.fName= paramhandle.GetLabel().fName;
        clabel.Set(label.GetStringLabel());
        GlobalData().Set(&clabel,"CurrentWriter");
 
    }
 
    // Calculate energy from calibration function
    const QCalibrationParameters& calParams = paramhandle.Get();
    energy = calParams.fFunction.Eval(amp);
    // eliminate 'nan's
    if (!(energy >=0) && !(energy<=0)) return;
 
    ev.Get<QDouble>("Energy") = energy;
 
    // Calculate estimated error on energy by evaluating function at limits of parameters
    double *params = calParams.fFunction.GetParameters();
    const double *paramErrors = calParams.fFunction.GetParErrors();
    int nParams = calParams.fFunction.GetNpar();
 
    TF1 calFunction;
    calParams.fFunction.Copy(calFunction);
 
    double *limitParams = new double[nParams];
    for (int i = 0; i < nParams; i++)
        limitParams[i] = params[i] + paramErrors[i];
    double energyUpper = calFunction.EvalPar(&amp, limitParams);
    for (int i = 0; i < nParams; i++)
        limitParams[i] = params[i] - paramErrors[i];
    double energyLower = calFunction.EvalPar(&amp, limitParams);
 
    ev.Get<QDouble>("EnergyUncertainty") = (energyUpper - energyLower) / 2;
}
 
void MApplyCalibration::Done()
{
    // called at the end of the event loop
}