MLockinDenoising.cc

Go to the documentation of this file.

#include "MLockinDenoising.hh"
#include "QEvent.hh"
#include "QEventList.hh"
#include "QRunData.hh"
#include "QRawEvent.hh"
#include "QGlobalHandle.hh"
#include "QPulseInfo.hh"
#include "QGlobalDataManager.hh"
#include "QRunDataHandle.hh"
#include "QAverageVectorC.hh"
#include "QRealComplexFFT.hh"
#include "QBool.hh"
#include "QBaseType.hh"
#include <iostream>
 
 
 
 
REGISTER_MODULE(MLockinDenoising)
 
using namespace Diana;
 
void MLockinDenoising::Init(QEvent& ev)
{
    fPulseLabel = GetString("PulseLabel","DAQ@Pulse",false);
    ev.Require<QHeader>("DAQ","Header");
    ev.RequireByLabel<QPulse>(fPulseLabel);
    ev.Require<QPulseInfo>("DAQ","PulseInfo");
    fNoiseCoherenceOwner = GetString("NoiseCoherenceOwner","NoiseFrequencyCoherence",false);
    fNoiseCoherenceInput = GetString("NoiseCoherenceInput","");
    ev.Add<QPulse>("Pulse");
    ev.Add<QBool>("IsDenoised");
 
 
}
 
void MLockinDenoising::Do(QEvent& ev, const QEventList& neighbours)
{
    const QHeader& hdr = ev.Get<QHeader>("DAQ","Header");
    const QPulseInfo& pi= ev.Get<QPulseInfo>("DAQ","PulseInfo");
    const QPulse& pulse = ev.GetByLabel<QPulse>(fPulseLabel);
    ev.Get<QPulse>("Pulse") = pulse;
    ev.Get<QBool>("IsDenoised") = false;
 
    const int chan = pi.GetChannelId();
 
 
 
    if(fMap.find(chan) == fMap.end()) {
        fMap[chan].BlackSheep = false;
        GlobalHandle<QInt> dHandle("Dataset");
        GlobalData().Get("",&dHandle,"");
        GlobalHandle<QAverageVectorC> Coherence_Handle("Coherence");
        Coherence_Handle.SetDataset(dHandle.Get());
        Coherence_Handle.SetChannel(chan);
        GlobalData().Get(fNoiseCoherenceOwner,&Coherence_Handle,fNoiseCoherenceInput);
        // try with an of run
        if(!Coherence_Handle.IsValid()) {
            Coherence_Handle.SetRun(hdr.GetRun());
            Coherence_Handle.SetDataset(Q_INT_DEFAULT);
            GlobalData().Get(fNoiseCoherenceOwner,&Coherence_Handle,fNoiseCoherenceInput);
        }
        if(!Coherence_Handle.IsValid()) {
            Warn("%s. Channel %d will not be processed",Coherence_Handle.GetError().GetDescription().c_str(),chan);
            fMap[chan].BlackSheep = true;
        } else {
            fMap[chan].fCoherence = Coherence_Handle.Get();
        }
    }
 
    ChannelInfo& chanInfo = fMap[chan];
    if (chanInfo.BlackSheep) return;
 
    if(neighbours.Size() != 1) return;
 
    const QEvent& neigh = neighbours[0];
 
    const QPulseInfo& piB= neigh.Get<QPulseInfo>("DAQ","PulseInfo");
    const int chanB = piB.GetChannelId();
    if(chan != chanB) return;
    
    const QVector& samples= pulse.GetSamples();
    const int size = samples.Size();
    const QVector& samplesB= neigh.GetByLabel<QPulse>(fPulseLabel).GetSamples();
    if(size != (int)samplesB.Size()) return;
 
    // check continuity
    int run = hdr.GetRun();
    QRunDataHandle rHandle(run);
    GlobalData().Get("",&rHandle,"");
    const QRunData& runData = rHandle.Get();
    const QChannelRunData& chanRunData = runData.GetChannelRunData(chan);
    const double samplPeriod = 1. / chanRunData.fSamplingFrequency*1e9;
 
    const QHeader& hdrB = neigh.Get<QHeader>("DAQ","Header");
    Long64_t lastSampleB = size-piB.GetMasterSample().GetSampleIndex()-1;
    lastSampleB +=  hdrB.GetTime().GetFromStartRunNs()/samplPeriod;
    Long64_t firstSample = hdr.GetTime().GetFromStartRunNs()/samplPeriod-pi.GetMasterSample().GetSampleIndex();
    if(firstSample-lastSampleB != 1) {
 //       Warn("Samples on event %d not contiguous (lastSampleB %lld, firstSample %lld",hdr.GetEventNumber(),lastSampleB,firstSample);
 //       return;
    }
 
 
    // Lockin denoising
    QRealComplexFFT transformer(size);
    QVectorC TransformedPulseB(size);
    QVectorC TransformedPulse(size);
 
    int err = transformer.TransformToFreq(samplesB, TransformedPulseB);
    if(err != 0) {
        Error("QFFTRealComplex::TransformToFreq: pulse B");
        return;
    } 
 
    err = transformer.TransformToFreq(samples, TransformedPulse);
    if(err != 0) {
        Error("QFFTRealComplex::TransformToFreq: pulse ");
        return;
    } 
 
    for(int i =1; i < size; i++) {
        TransformedPulse[i] -= TransformedPulseB[i]*chanInfo.fCoherence[i];
    }
    
 
    QVector output;
    transformer.TransformFromFreq(TransformedPulse,output);
 
    ev.Get<QPulse>("Pulse").Clear();
    ev.Get<QPulse>("Pulse").SetSamplesADC(output);
    ev.Get<QPulse>("Pulse").Validate();
    ev.Get<QBool>("IsDenoised") = true;
 
    
 
}
 
void MLockinDenoising::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.
     */
}