MOptimumFilter.cc

#include "MOptimumFilter.hh"
#include "QEvent.hh"
#include "QEventList.hh"
#include "QRawEvent.hh"
#include "QBaseModule.hh"
#include "QSampleInfo.hh"
#include "QRunDataHandle.hh"
#include "QBaseType.hh"
#include "QBaselineData.hh"
#include "QAveragePulseHandle.hh"
#include "QAverageNoiseHandle.hh"
#include "QCOFParametersHandle.hh"
#include "QOptimumFilter.hh"

REGISTER_MODULE(MOptimumFilter)

    using namespace Cuore;

void MOptimumFilter::Init(QEvent& ev)
{
    fAvgNoiseInput = GetString("AvgNoiseInput","");
    fAvgPulseInput= GetString("AvgPulseInput","");
    fAvgPulseOwner = GetString("AvgPulseOwner","AveragePulses",false);
    fAvgNoiseOwner = GetString("AvgNoiseOwner","NoiseAvgPowerSpectrum",false);

    std::string searchMode = GetString("AmplitudeMode","AbsoluteMaximum");
    if(searchMode == "AbsoluteMaximum") fJitterMode = QOptimumFilter::J_ABSMAX;
    else if(searchMode == "NoSearch") fJitterMode = QOptimumFilter::J_NOJITTER;
    else if(searchMode == "LocalMaximum") fJitterMode = QOptimumFilter::J_LOCMAX;
    else if(searchMode == "FixedDistance") fJitterMode = QOptimumFilter::J_FIXED;
    else Panic("AmplitudeMode %s not available",searchMode.c_str());

    fInterpolationOn = GetBool("InterpolationOn",true,true);
    fMaxJitter = GetInt("MaxJitter",-1,false);
    fParametersOutput = GetString("ParametersOutput","",false);
    fDifferentiationOn = GetBool("DifferentiationOn",true);

    ev.Add<QDouble>("Amplitude");
    ev.Add<QDouble>("Integral");
    ev.Add<QDouble>("ChiSquare");
    ev.Add<QDouble>("Jitter");
    ev.Add<QVector>("FilteredSamples").SetWrite(GetBool("SaveSamples",false,false));

    fPulseLabel = GetString("PulseLabel","DAQ@Pulse",false);
    ev.RequireByLabel<QPulse>(fPulseLabel);
    ev.Require<QPulseInfo>("DAQ","PulseInfo");
    ev.Require<QHeader>("DAQ","Header");

}

//void MOptimumFilter::Do(QEvent& ev, const QEventList& neighbours)
void MOptimumFilter::Do(QEvent& ev)
{
    const QPulseInfo& pulseInfo = ev.Get<QPulseInfo>("DAQ","PulseInfo");
    const QHeader& header = ev.Get<QHeader>("DAQ","Header");
    const int chan = pulseInfo.GetChannelId();

    QRunDataHandle rHandle(header.GetRun());
    GlobalData().Get("",&rHandle,"");

    const QRunData& runData = rHandle.Get();
    const QChannelRunData& channelRunData = runData.GetChannelRunData(chan);
    double samplingFreq = channelRunData.fSamplingFrequency;
    double ADC2mV = channelRunData.fADC2mV;

    if(fMap.find(chan) == fMap.end()) {

        GlobalHandle<QInt> dHandle("Dataset");
        GlobalData().Get("",&dHandle,"");

        fMap[chan].BlackSheep = false;

        QAverageNoiseHandle AVG_Noise_Handle(chan);
        AVG_Noise_Handle.SetDataset(dHandle.Get());
        GlobalData().Get(fAvgNoiseOwner,&AVG_Noise_Handle,fAvgNoiseInput);
        // try with an of run
        if(!AVG_Noise_Handle.IsValid()) {
            AVG_Noise_Handle.SetRun(header.GetRun());
            AVG_Noise_Handle.SetDataset(Q_INT_DEFAULT);
            GlobalData().Get(fAvgNoiseOwner,&AVG_Noise_Handle,fAvgNoiseInput);
        }
        if(!AVG_Noise_Handle.IsValid()) {
            Warn("%s. Channel %d will not be processed",AVG_Noise_Handle.GetError().GetDescription().c_str(),chan);
            fMap[chan].BlackSheep = true;
        }
        QAveragePulseHandle AVG_Pulse_Handle(chan);
        AVG_Pulse_Handle.SetDataset(dHandle.Get());
        GlobalData().Get(fAvgPulseOwner,&AVG_Pulse_Handle,fAvgPulseInput);

        if(!AVG_Pulse_Handle.IsValid()) {
            AVG_Pulse_Handle.SetRun(header.GetRun());
            AVG_Pulse_Handle.SetDataset(Q_INT_DEFAULT);
            GlobalData().Get(fAvgPulseOwner,&AVG_Pulse_Handle,fAvgPulseInput);
        }

        if(!fMap[chan].BlackSheep && !AVG_Pulse_Handle.IsValid()) {
            Warn("%s. Channel %d will not be processed",AVG_Pulse_Handle.GetError().GetDescription().c_str(),chan);
            fMap[chan].BlackSheep = true;
        }
        if(!fMap[chan].BlackSheep) {

            fMap[chan].ADC2Time = 1000./samplingFreq; //ms

            const QVector& avg_noise = AVG_Noise_Handle.Get();
            const QVector& avg_pulse = AVG_Pulse_Handle.Get();

            fMap[chan].of = new QOptimumFilter(avg_pulse,  avg_noise,fMaxJitter , fDifferentiationOn,false );
            Info("Channel %d: Maximum Jitter = +/- %d, Cutoff = %.0f",chan, fMap[chan].of->GetMaxJitter(),  fMap[chan].of->GetCutOffFrequency()*samplingFreq);

            QCOFData cofData;
            cofData.fResolutionmV = fMap[chan].of->GetFilteredNoiseRMS()*ADC2mV;

            QCOFParametersHandle cofhandle;
            cofhandle.SetChannel(chan);         
            cofhandle.SetDataset(dHandle.Get());
            cofhandle.Set(cofData);

            // if fStoreParametersInDB is true call COFhandle 
            if (!fParametersOutput.empty()) {
                cofhandle.SetAPLabel(AVG_Pulse_Handle.GetLabel());  
                cofhandle.SetAPVersion(AVG_Pulse_Handle.GetVersion());  
                if(AVG_Pulse_Handle.GetDataset()!=Q_INT_DEFAULT){
                    cofhandle.SetAPValidityKind("data_set");
                }
                else {
                    cofhandle.SetAPValidityKind("run");
                }
                cofhandle.SetNPSLabel(AVG_Noise_Handle.GetLabel());
                cofhandle.SetNPSVersion(AVG_Noise_Handle.GetVersion());  
                if(AVG_Noise_Handle.GetDataset()!=Q_INT_DEFAULT){
                    cofhandle.SetNPSValidityKind("data_set");  
                }
                else {
                    cofhandle.SetNPSValidityKind("run");
                }
            }
            GlobalData().Set(&cofhandle,fParametersOutput);

        }
    }

    ChannelInfo& chanInfo = fMap[chan];
    if (chanInfo.BlackSheep) return;

    const QPulse& rawPulse = ev.GetByLabel<QPulse>(fPulseLabel);
    const QVector& Samples = rawPulse.GetSamples();
    QError err = chanInfo.of->Filter(Samples);
    if(err!=QERR_SUCCESS) {
        Debug("Event number %d %s",header.GetEventNumber(),err.GetDescription().c_str());
        return;
    }

    if(pulseInfo.GetIsNoise()) {
        err = chanInfo.of->SetJitter(QOptimumFilter::J_NOJITTER);
    }  else if(fJitterMode == QOptimumFilter::J_LOCMAX) {
        err = chanInfo.of->SetJitter(QOptimumFilter::J_LOCMAX,pulseInfo.GetMasterSample().GetSampleIndex());
    }  else {
        err = chanInfo.of->SetJitter(fJitterMode);

    }
    if(err!=QERR_SUCCESS) {
        Debug("Event number %d %s",header.GetEventNumber(),err.GetDescription().c_str());
        return;
    }

    double a1, chi2, integral, jitter;
    if(!fInterpolationOn || pulseInfo.GetIsNoise()) {
        err = chanInfo.of->Get(jitter,chi2, a1, integral);
    } else {
        err = chanInfo.of->GetInterpolated(jitter,chi2, a1,integral);
    }
    if(err!=QERR_SUCCESS) {
        Debug("Event number %d %s",header.GetEventNumber(),err.GetDescription().c_str());
        return;
    }

    ev.Get<QDouble>("Amplitude") = a1*ADC2mV;
    ev.Get<QDouble>("Integral") = integral*ADC2mV/samplingFreq;

    ev.Get<QDouble>("ChiSquare") = chi2;
    ev.Get<QDouble>("Jitter") = jitter*fMap[chan].ADC2Time;

    ev.Get<QVector>("FilteredSamples") = chanInfo.of->GetFilteredShifted();
}  


void MOptimumFilter::Done()
{
    /* This method is called at the end of the event loop. 
     * Here you can:
     *
     * 1) Operate on the sequence execution parameters(see QBaseModule.hh).
     *
     * 2) Read/Write global data.
     */
}