MNoiseFrequencyCoherence.cc

Go to the documentation of this file.

#include "MNoiseFrequencyCoherence.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 "QBaseType.hh"
#include <iostream>
 
 
 
REGISTER_MODULE(MNoiseFrequencyCoherence)
 
using namespace Diana;
 
void MNoiseFrequencyCoherence::Init(QEvent& ev)
{
    fOutput = GetString("Output",Q_STRING_DEFAULT);
    std::string window = GetString("WindowType","welch");
    fWindowType = QFFT::StrToWindowType(window);
    if(fWindowType == QFFT::WT_None) {
        Warn("Window type not defined, using rectangular window");
        fWindowType = QFFT::WT_Rect;
    }
    else {
        Info("Window type: %d (%s)", fWindowType,window.c_str());
    }
    if(fWindowType != QFFT::WT_Rect)
        fCoherentGain = GetInt("WindowCoherentGain",2);
    else fCoherentGain = 2;
    fValidityKind = GetString("ValidityKind","dataset",false);
    if(fValidityKind != "dataset" && fValidityKind != "run") {
        Panic("ValidityKind must be set to dataset or run");
    }
 
    fValidityStart = GetInt("ValidityStart",Q_INT_DEFAULT,false);
    fValidityEnd = GetInt("ValidityEnd",Q_INT_DEFAULT,false);
    fPulseLabel = GetString("PulseLabel","DAQ@Pulse",false);
    ev.Require<QHeader>("DAQ","Header");
    ev.RequireByLabel<QPulse>(fPulseLabel);
    ev.Require<QPulseInfo>("DAQ","PulseInfo");
 
    fSRuns.clear();
    fCurrentRun = -1;
    fChannelInfos.clear();
}
 
 
 
void MNoiseFrequencyCoherence::Do(QEvent& ev, const QEventList& neighbours)
{
    if(neighbours.Size() != 1) return;
    const QEvent& neigh = neighbours[0];
 
    const QPulseInfo& pi= ev.Get<QPulseInfo>("DAQ","PulseInfo");
    const QPulseInfo& piB= neigh.Get<QPulseInfo>("DAQ","PulseInfo");
 
    const int chan = pi.GetChannelId();
    const int chanB = piB.GetChannelId();
    if(chan != chanB) return;
    
    const QVector& samples= ev.GetByLabel<QPulse>(fPulseLabel).GetSamples();
    const int size = samples.Size();
    const QVector& samplesB= neigh.GetByLabel<QPulse>(fPulseLabel).GetSamples();
    if(size != (int)samplesB.Size()) return;
 
    // check continuity
    const QHeader& hdr = ev.Get<QHeader>("DAQ","Header");
    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;
    }
 
 
    if (run != fCurrentRun) {
        fCurrentRun = run;
        fSRuns.push_back(fCurrentRun);
    }
 
    if(fChannelInfos.find(chan) == fChannelInfos.end()) {
        ChannelInfo& chInfo = fChannelInfos[chan];
        chInfo.fNumEvents = 0;
        chInfo.fCovariance.Resize(size);
        chInfo.fPowerB.Resize(size);
        chInfo.fPower.Resize(size);
        chInfo.fCovariance.Initialize(0,0);
        chInfo.fPowerB.Initialize(0);
        chInfo.fPower.Initialize(0);
        chInfo.fTransformer = new QRealComplexFFT(size);
        chInfo.fTransformer->SetWindowType(fWindowType, fCoherentGain);
    }
    ChannelInfo& chInfo = fChannelInfos[chan];
 
 
    QVectorC TransformedPulseB(size);
    QVectorC TransformedPulse(size);
    int err = Transform(chInfo.fTransformer,samplesB, TransformedPulseB);
    if(err != 0) {
        Error("QFFTRealComplex::TransformToFreq: pulse B");
        return;
    } 
    TransformedPulseB[0].SetRe(0);
    TransformedPulseB[0].SetIm(0);
    err = Transform(chInfo.fTransformer,samples, TransformedPulse);
    if(err != 0) {
        Error("QFFTRealComplex::TransformToFreq: pulse ");
        return;
    } 
    TransformedPulse[0].SetRe(0);
    TransformedPulse[0].SetIm(0);
 
 
    chInfo.fPowerB += TransformedPulseB.GetModulusSquare();
    chInfo.fPower += TransformedPulse.GetModulusSquare();
    chInfo.fCovariance += TransformedPulse.Mult(TransformedPulseB.Conj());
    chInfo.fNumEvents++;
 
}
 
void MNoiseFrequencyCoherence::Done()
{
    std::map<int,ChannelInfo>::iterator iter = fChannelInfos.begin();
    while (iter != fChannelInfos.end() ) {
        ChannelInfo& chanInfo = iter->second;
        if(chanInfo.fNumEvents == 0) {
            iter++;
            continue;
        }
 
        int channel = iter->first;
        chanInfo.fPower /= chanInfo.fNumEvents;
        chanInfo.fPowerB /= chanInfo.fNumEvents;
        chanInfo.fCovariance /= chanInfo.fNumEvents;
 
        QVector fAvgPower(chanInfo.fPower.Size());
        for(size_t f = 0; f < chanInfo.fPower.Size(); f++) {
            fAvgPower[f] = sqrt(chanInfo.fPower[f]*chanInfo.fPowerB[f]);
        }
 
        QAverageVectorC coherence;
        coherence = chanInfo.fCovariance.Div(fAvgPower);
        coherence[0].Set(0,0);
        coherence.SetSourceRuns(fSRuns); 
        coherence.SetNumEvents(chanInfo.fNumEvents); 
        //std::cout<<"channel: "<<channel<<" coh: "<<coherence<<std::endl;
        // Fill Output
 
        GlobalHandle<QAverageVectorC> handle("Coherence");
        handle.SetChannel(channel);
        if(fValidityKind=="dataset"){
            GlobalHandle<QInt> dHandle("Dataset");
            GlobalData().Get("",&dHandle,"");
            int dataset = dHandle.Get();
            handle.SetDataset(dataset);
        }
        else if (fValidityKind=="run"){
            if (fSRuns.size()==1){
                handle.SetRun(fSRuns[0]);
            }
            else Error("ValidityKind set to run is not permitted when reading multiple runs.");
        }
        else Error("ValidityKind can be set only to run or dataset. Check cfg.");
 
        handle.Set(coherence);
        GlobalData().Set(&handle,fOutput);
        iter++;
 
    }
}
 
int MNoiseFrequencyCoherence::Transform(Diana::QRealComplexFFT* transformer, const Diana::QVector& samples, Diana::QVectorC& out)
{
    const size_t size = samples.Size();
    QVector zeroMeanPulse(size), meanPulse(size);
    double mean = samples.Sum(size)/size;
    meanPulse.Initialize(mean);
    zeroMeanPulse = samples;
    zeroMeanPulse -= meanPulse;
    return transformer->TransformToFreq(zeroMeanPulse, out);
}