MDecorrelator.cc

Go to the documentation of this file.

#include "MDecorrelator.hh"
#include "QEvent.hh"
#include "QEventList.hh"
#include "QPulse.hh"
#include "QVectorView.hh"
#include "QPulseInfo.hh"
#include "QRealComplexFFT.hh"
#include "QCountPulsesData.hh"
#include "QBaselineData.hh"
#include "QBool.hh"
#include "QBaseType.hh"
#include "QStdVector.hh"
#include "QHeader.hh"
#include "QMultiChannelDecorrelatorOF.hh"
#include "QAveragePulseHandle.hh"
#include "QRunDataHandle.hh"
#include <algorithm>
 
using namespace Diana;
using namespace std;
 
 
REGISTER_MODULE(MDecorrelator)
 
 
void MDecorrelator::Init(QEvent& ev)
{
    ev.Add<QPulse>("Pulse").SetWrite(GetBool("SaveSamples",true));
    ev.Add<QBool>("Decorrelated");
    ev.Require<QPulseInfo>("DAQ","PulseInfo");
    ev.Require<QHeader>("DAQ","Header");
    ev.Require<QPulse>("DAQ","Pulse");
    fFirstCall = true;
    fUseOF = GetBool("UseOF",false);
    if(fUseOF) {
        ev.Add<QDouble>("OFAmplitude");
        ev.Add<QDouble>("OFMaximumPosition");
    }
}
 
void MDecorrelator::Do(QEvent& ev, const QEventList& neighbours)
{
    const QPulseInfo& pi = ev.Get<QPulseInfo>("DAQ","PulseInfo");
    const int channel = pi.GetChannelId();
    const int run = ev.Get<QHeader>("DAQ","Header").GetRun();
    if(fFirstCall) {
        std::string corrInput = GetString("CovarianceInput",Q_STRING_DEFAULT);
 
        GlobalHandle<QInt> dataset("Dataset");
        GlobalData().Get("",&dataset,"");
        GlobalHandle<QChannelCovariance> CovGetter("Covariance");
        CovGetter.SetDataset(dataset.Get());
        GlobalData().Get("NoiseCrossPowerSpectrum",&CovGetter,corrInput);
        if(!CovGetter.IsValid()) {
            CovGetter.SetDataset(Q_INT_DEFAULT);
            CovGetter.SetRun(run);
            GlobalData().Get("NoiseCrossPowerSpectrum",&CovGetter,corrInput);
        }
 
        fCovariance = CovGetter.Get();
        fFirstCall = false;
    }
 
 
    if(fDecorrelators.find(channel) == fDecorrelators.end()) {
        QMultiChannelDecorrelator* newdeco = 0;
        bool pDeco=false,pOF=false,pOFMultiDim=false;
        if(fUseOF) {
            pOF = true;
            std::string apInput = GetString("AveragePulseInput",Q_STRING_DEFAULT);
            std::string apOwner = GetString("AveragePulseOwner","AveragePulses",false);
            std::string relativesInput = GetString("RelativeChannelsInput","",false);
            std::vector<QVector> averagePulses;
            std::vector<int> relatives;
 
            GlobalHandle<QStdVector<int> > relativesHandle("RelativeChannels");
            relativesHandle.SetChannel(channel);
            GlobalData().Get("User",&relativesHandle,relativesInput);
            if(relativesHandle.IsValid() && !relativesHandle.Get().empty()) {
                relatives = relativesHandle.Get();
            }
            if(!relatives.empty()) {
                pOF = false;
                pOFMultiDim = true;
            }
            relatives.insert(relatives.begin(),channel);
            for(size_t i =0; i < relatives.size(); i++) {
                const int rchan = relatives[i];
                QVector ap;
                if(apOwner == "Old") {
                    GlobalHandle<QVector> apHandle("Average");
                    apHandle.SetChannel(rchan);
                    GlobalData().Get(apOwner,&apHandle,apInput);
                    ap = apHandle.Get();
                } else {
                    QAveragePulseHandle apHandle(rchan);
                    GlobalHandle<QInt> dataset("Dataset");
                    GlobalData().Get("",&dataset,"");
                    apHandle.SetDataset(dataset.Get());
                    GlobalData().Get(apOwner,&apHandle,apInput);
                    ap = apHandle.Get();
                }
                averagePulses.push_back(ap);
            }
            newdeco = new QMultiChannelDecorrelatorOF(relatives,&fCovariance,averagePulses);
 
        } 
        std::string sideInput = GetString("SideChannelsInput",Q_STRING_DEFAULT);
        GlobalHandle<QStdVector<int> > sidesHandle("CorrelatedChannels");
        sidesHandle.SetChannel(channel);
        GlobalData().Get("User",&sidesHandle,sideInput);
        if(sidesHandle.IsValid()) {
            pDeco = true;
            if(!pOF && !pOFMultiDim) newdeco = new QMultiChannelDecorrelator(channel,&fCovariance);
            newdeco->SetSideChannels(sidesHandle.Get());
        }
 
       if(GetBool("ConvertTomV",false,false)) {
            QRunDataHandle rHandle(run);
            GlobalData().Get("",&rHandle,"");
            double ADC2mV = rHandle.Get().GetChannelRunData(channel).fADC2mV;
            if(newdeco) {
                newdeco->SetScale(ADC2mV);
                // forse re-initialization
                newdeco->SetSideChannels(newdeco->GetSideChannels());
            }
        }
        string msg;
        if(newdeco) {
            
            stringstream msgstr;
            msgstr<<", OrigRMS="<<newdeco->GetOriginalResolution();
            msgstr<<", FiltRMS="<<newdeco->GetFilteredResolution();
            msg = msgstr.str();
        }
 
        fDecorrelators[channel].fDecorrelator = newdeco;
        fDecorrelators[channel].fUsed = false;
 
        Info("Channel %04d: Deco %d, OF1D %d, OFMD %d%s",channel,pDeco,pOF,pOFMultiDim,msg.c_str());
 
    }
 
    QMultiChannelDecorrelator* deco = fDecorrelators[channel].fDecorrelator;
    QVector decoPulse;
    QPulse& decorrelated_pulse = ev.Get<QPulse>("Pulse");
    const QPulse& pulse = ev.Get<QPulse>("DAQ","Pulse");
    
 
    ev.Get<QBool>("Decorrelated") = false;
 
    if(deco != 0) {
        std::vector<const QVector*> sideSamples;
        const std::vector<int>& sideChannels = deco->GetSideChannels();
        for(size_t ich = 0; ich < sideChannels.size(); ich++) {
            const QVector& sideSample = GetNeighbourSamples(sideChannels[ich],neighbours);
            if(sideSample.Size() > 0) {
                sideSamples.push_back(&sideSample);
            }
        }
        if(sideSamples.size() != sideChannels.size()) {
            decoPulse = pulse.GetSamples();
            decorrelated_pulse.SetSamplesADC(decoPulse);
            return;
        }
        std::vector<const QVector*> samples;
        samples.push_back(&pulse.GetSamples());
        const std::vector<int>& channels = deco->GetChannels();
        for(size_t ich = 1; ich < channels.size(); ich++) {
            const QVector& sample = GetNeighbourSamples(channels[ich],neighbours);
            if(sample.Size() > 0) {
                samples.push_back(&sample);
            }
        }
        if(samples.size() != channels.size()) {
            decoPulse = pulse.GetSamples();
            decorrelated_pulse.SetSamplesADC(decoPulse);
            return;
        }
 
        deco->Filter(samples,sideSamples,decoPulse);
        if(fUseOF) {
            double amp = Q_DOUBLE_DEFAULT;
            double maxpos = Q_DOUBLE_DEFAULT;
            if(pi.GetIsNoise()) {
                amp = decoPulse[pi.GetMasterSample().GetSampleIndex()];
                maxpos = pi.GetMasterSample().GetSampleIndex();
            } else {
                amp = decoPulse.GetMax();
                maxpos = decoPulse.GetMaxIndex();
            }
            ev.Get<QDouble>("OFAmplitude") = amp;
            ev.Get<QDouble>("OFMaximumPosition") = maxpos;
        }
        fDecorrelators[channel].fUsed = true;
        ev.Get<QBool>("Decorrelated") = true;
        decorrelated_pulse.SetSamplesADC(decoPulse);
    } else {
        decoPulse = pulse.GetSamples();
        decorrelated_pulse.SetSamplesADC(decoPulse);
        return;
    }
 
    //decorrelated_pulse.SetSamplesADC(decoPulse);
    return;
}   
 
void MDecorrelator::Done()
{
    std::map<int,ChannelInfo>::const_iterator iter = fDecorrelators.begin();
    while(iter != fDecorrelators.end()) {
        if(iter->second.fDecorrelator && !iter->second.fUsed) Warn("Channel %d was never decorrelated because the expected neighbours were not present",iter->first);
        iter++;
    }
}
 
const QVector& MDecorrelator::GetNeighbourSamples(const int channel, const QEventList& neighbours)
{
    for(size_t i = 0; i < neighbours.Size(); i++) {
        const QPulseInfo& pi = neighbours[i].Get<QPulseInfo>("DAQ","PulseInfo");
        if(pi.GetChannelId() == channel) {
            const QPulse& pulse = neighbours[i].Get<QPulse>("DAQ","Pulse");
            return pulse.GetSamples(); 
        }
    }
    return fEmpty;
}