MDownSampling.cc

Go to the documentation of this file.

#include "MDownSampling.hh"
#include "QEventAssembler.hh"
#include "QRawEvent.hh"
#include "QBaseType.hh"
#include "QRunDataHandle.hh"
#include "QPulse.hh"
#include "QVector.hh"
#include <cmath>
#include "QFir.hh"
 
 
 
 
REGISTER_MODULE(MDownSampling)
 
using namespace Diana;
using namespace std;
 
QDriver::Action MDownSampling::Init(QEventAssembler& eva)
{
    
    
    fCutOffFrequency = GetDouble("CutOffFrequency",300000);
    fAttenuation = GetDouble("AttenuationDB",60);
    fPulseLabel = GetString("PulseLabel","DAQ@Pulse",false);
    fOverWritePulse = GetBool("OverwritePulse",true);
    fDownSampling = GetDouble("DownSampling",1.,false);
    if(fDownSampling < 1.) Panic("Downsampling factor (%f) must be >= 1",fDownSampling);
    
    
    QEvent& ev = eva.GetEvent();
    ev.Require<QHeader>("DAQ","Header");
    ev.Require<QPulseInfo>("DAQ","PulseInfo");
    ev.RequireByLabel<QPulse>(fPulseLabel);
    
    if(!fOverWritePulse) {
        ev.Add<QPulse>(fPulseLabel.name);
        ev.Add<QPulseInfo>("PulseInfo");
    }
 
    Info("Cutoff set to %f",fCutOffFrequency);
    return ACT_NEXTEV;
}
 
QDriver::Action MDownSampling::Do(QEventAssembler& eva)
{
    QEvent& ev = eva.GetEvent();
    QPulseInfo& pi = eva.Get<QPulseInfo>("DAQ","PulseInfo");
    const int channel = pi.GetChannelId();
 
    if(fFilters.find(channel) == fFilters.end()) {
        const QHeader& header = ev.Get<QHeader>("DAQ","Header");
        QRunDataHandle rHandle(header.GetRun());
        GlobalData().Get("DAQ",&rHandle,"");
        const QRunData& runData = rHandle.Get();
        const QChannelRunData& channelRunData = runData.GetChannelRunData(channel);
        const double samplFreq = channelRunData.fSamplingFrequency;
        const double ADC2mV = channelRunData.fADC2mV;
        double fracFreq = fCutOffFrequency/samplFreq;
 
        int M = fAttenuation * 0.025/fracFreq;
        if(M%2 == 0) M++;
        QCFirData configuration;
        configuration.fM = M;
        configuration.fCutoff = fracFreq;
        configuration.fAttDB = fAttenuation;
        configuration.fMethod = QCFirData::KayserBessel;
        QFir* fir = new  QFir(configuration);
        size_t trigPos = pi.GetMasterSample().GetSampleIndex();
        if(fir->GetFilterReduction() > trigPos/2) {
            Panic("Filter reduction (%d) is greater than triggerPosition/2 (%d), try a lower attenuation or a higher cutoff frequency (Size of filter scales linearly with them)",fir->GetFilterReduction(),trigPos/2);
        }
 
        if(!fir) Panic("Cannot create FIR for channel %d", channel);
 
 
        ChannelInfo info;
        info.fir = fir;
        info.numberOfSamples = channelRunData.fNumberOfSamples - fir->GetFilterReduction();
 
        //undersampling
        info.R = fDownSampling;//samplFreq/(fCutOffFrequency*2)/4;//??????????
        //if(info.R == 0) Panic("Undersampling failed. Cutoff frequency (%f) must be lower than (sampling frequency/8)= %f",fCutOffFrequency,samplFreq/8);
        
        //info.numberOfSamples /= (info.R);
        info.numberOfSamples /= fDownSampling;
        //if(info.numberOfSamples %2 != 0) info.numberOfSamples++;
 
        fFilters[channel] = info;
        GlobalHandle<QCFirData>  cfh("FirData");
        cfh.SetChannel(channel);
        cfh.SetRun(header.GetRun());
        cfh.Set(configuration);
        GlobalData().Set(&cfh,"CurrentWriter");
 
        if(fOverWritePulse) {
            // overwrite rundata
            void* rdv = (void*)&channelRunData;
            QChannelRunData* rd = (QChannelRunData*)rdv;
            rd->fNumberOfSamples = info.numberOfSamples;
            rd->fADC2mV = ADC2mV;
            rd->fSamplingFrequency /= info.R;
            GlobalData().Set(&rHandle,"CurrentWriter");
 
        }
 
    }
 
 
    QPulse& oldPulse =  eva.Get<QPulse>(fPulseLabel.owner.c_str(),fPulseLabel.name.c_str());
    const QVector& input = oldPulse.GetSamples();
 
    const ChannelInfo& info = fFilters[channel];
    fOutput.Clear();
    info.fir->Filter(input,fOutput);
 
    //cout<<"Input Event: "<<input.Size()<<"\t Output Event: "<<fOutput.Size()<<"\t FIR Reduction: "<<info.fir->GetFilterReduction()<<endl;
    int samps = 0;
    for(size_t i = 0; i < info.numberOfSamples; i++) {
        if(i*info.R<fOutput.Size()){
            samps++;
            fOutput[i] = fOutput[i*info.R];
        }
    }
 
    fOutput.Resize(samps);
 
    int ms = pi.GetMasterSample().GetSampleIndex()-info.fir->GetFilterReduction();
    
    ms /= info.R;
 
    if(fOverWritePulse) {
        pi.GetMasterSample().SetSampleIndex(ms);
        oldPulse.SetSamplesADC(fOutput);
    } else {
        ev.Get<QPulseInfo>("PulseInfo") = pi;
        ev.Get<QPulseInfo>("PulseInfo").GetMasterSample().SetSampleIndex(ms);
        QPulse& filteredPulse = ev.Get<QPulse>(fPulseLabel.name.c_str());
        filteredPulse.SetSamplesADC(fOutput);
    }
 
    //MV FIXME: add other sampleInfo
 
 
    return ACT_NEXTEV;
}
 
QDriver::Action MDownSampling::Done()
{
    return ACT_NEXTEV;
}