MOptimumFilter.cc

Go to the documentation of this file.

#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 Diana;
 
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);
   
    double defaultChi2Threshold1 = 0.1;
    fChi2Threshold1 = GetDouble("Chi2Threshold1",defaultChi2Threshold1,false);
    if(fChi2Threshold1 < 0 || fChi2Threshold1 > 1) {
        Error("Chi2Threshold1 value %f must be between 0 and 1, resetting to default %f",fChi2Threshold1, defaultChi2Threshold1);
        fChi2Threshold1 = defaultChi2Threshold1;
    }
 
 
    fMaxJitter = GetInt("MaxJitter",-1,false);
    fParametersOutput = GetString("ParametersOutput","",false);
    fDifferentiationOn = GetBool("DifferentiationOn",true);
    fTomV = GetBool("ConvertTomV",true,false);
 
    ev.Add<QDouble>("Amplitude");
    ev.Add<QDouble>("Integral");
    ev.Add<QDouble>("ChiSquare");
    ev.Add<QDouble>("HighSNRChiSquare");
    ev.Add<QDouble>("Jitter");
    ev.Add<QDouble>("ChiLeft");
    ev.Add<QDouble>("ChiRight");
 
    fSaveSamples =  GetBool("SaveSamples",false,false);
    fShiftFilteredSamples = GetBool("ShiftFilteredSamples",true,false);
    
    ev.Add<QVector>("FilteredSamples").SetWrite(fSaveSamples);
    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 = 1;
    if(fTomV) 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 [samples], Cutoff = %.0f [Hz]",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;
    }
    if(fShiftFilteredSamples)
        ev.Get<QVector>("FilteredSamples") = chanInfo.of->GetFilteredShifted();
    else
        ev.Get<QVector>("FilteredSamples") = chanInfo.of->GetFiltered();
 
    double a1, chi2, integral, jitter,chileft,chiright;
    if(!fInterpolationOn || pulseInfo.GetIsNoise()) {
        err = chanInfo.of->Get(jitter,chi2, a1, integral,chileft,chiright);
    } else {
        err = chanInfo.of->GetInterpolated(jitter,chi2, a1,integral,chileft,chiright);
    }
    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<QDouble>("ChiLeft") = chileft;
    ev.Get<QDouble>("ChiRight") = chiright;
 
    double highSNRChi21 = Q_DOUBLE_DEFAULT;
    err = chanInfo.of->GetHighSNRChi2(fChi2Threshold1,highSNRChi21);
    if(err != QERR_SUCCESS) {
        Debug("Event number %d %s",header.GetEventNumber(),err.GetDescription().c_str());
        return;
    }
    ev.Get<QDouble>("HighSNRChiSquare") = highSNRChi21;
 
 
 
 
    return;
}  
 
 
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.
     */
}