MNoiseCrossFrequencyCovariance.cc

Go to the documentation of this file.

#include "MNoiseCrossFrequencyCovariance.hh"
#include "QEvent.hh"
#include "QEventList.hh"
#include "QRealComplexFFT.hh"
#include "QRunData.hh"
#include "QRawEvent.hh"
#include "QGlobalHandle.hh"
#include "QPulseInfo.hh"
#include "QGlobalDataManager.hh"
#include "QRunDataHandle.hh"
#include "QStdVector.hh"
#include "QBaseType.hh"
#include "QVectorView.hh"
#include "QAverageNoiseHandle.hh"
#include <iostream>
 
 
 
REGISTER_MODULE(MNoiseCrossFrequencyCovariance)
 
using namespace Diana;
 
void MNoiseCrossFrequencyCovariance::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;
    fUseMainEvent = GetBool("UseMainEvent",false);
    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);
    fSingleChannel = GetBool("SingleChannel",false,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;
}
 
 
 
void MNoiseCrossFrequencyCovariance::Do(QEvent& ev, const QEventList& neighbours)
{
 
    QEventList allchans;
    if(fUseMainEvent) allchans.Push(&ev);
    for(size_t ch = 0; ch < neighbours.Size(); ch++) {
        allchans.Push(&neighbours[ch]);
    }
    int run = ev.Get<QHeader>("DAQ","Header").GetRun();
    if (run != fCurrentRun) {
        fCurrentRun = run;
        fSRuns.push_back(fCurrentRun);
    }
 
 
    if(fTimeLength == 0) {
        fTimeLength = ev.GetByLabel<QPulse>(fPulseLabel).GetSamples().Size()/2;
        if(fTimeLength %2 != 0) {
            fTimeLength -= 1;
            Info("Window size reduced to %d",fTimeLength);
        }
        fTransformer = new QRealComplexFFT(fTimeLength);
        fTransformer->SetWindowType(fWindowType, fCoherentGain);
    }
 
    QVectorC TransformedPulseIB(fTimeLength);
    QVectorC TransformedPulseIF(fTimeLength);
    QVectorC TransformedPulseJB(fTimeLength);
    QVectorC TransformedPulseJF(fTimeLength);
    for(size_t i = 0;  i < allchans.Size(); i++) {
        const int chanI = allchans[i].Get<QPulseInfo>("DAQ","PulseInfo").GetChannelId();
        const QVector& samplesIorig =  allchans[i].GetByLabel<QPulse>(fPulseLabel).GetSamples();
 
        QVectorConstView samplesIviewB(samplesIorig,0,fTimeLength);
        const QVector& samplesIB = samplesIviewB.GetVector();
        int err = Transform(samplesIB, TransformedPulseIB);
        if(err != 0) {
            Error("QFFTRealComplex::TransformToFreq: pulse I");
            return;
        } 
        TransformedPulseIB[0].SetRe(0);
        TransformedPulseIB[0].SetIm(0);
 
        QVectorConstView samplesIviewF(samplesIorig,fTimeLength,fTimeLength);
        const QVector& samplesIF = samplesIviewF.GetVector();
        err = Transform(samplesIF, TransformedPulseIF);
        if(err != 0) {
            Error("QFFTRealComplex::TransformToFreq: pulse I");
            return;
        } 
        TransformedPulseIF[0].SetRe(0);
        TransformedPulseIF[0].SetIm(0);
 
        for(size_t j = i ;  j < (fSingleChannel ? i+1: allchans.Size()); j++) {
            const int chanJ = allchans[j].Get<QPulseInfo>("DAQ","PulseInfo").GetChannelId();
            const QVector& samplesJorig =  allchans[j].GetByLabel<QPulse>(fPulseLabel).GetSamples();
            QVectorConstView samplesJviewF(samplesJorig,fTimeLength,fTimeLength);
            const QVector& samplesJF = samplesJviewF.GetVector();
            err = Transform(samplesJF, TransformedPulseJF);
            if(err != 0) {
                Error("QFFTRealComplex::TransformToFreq: pulse J");
                return;
            } 
            TransformedPulseJF[0].SetRe(0);
            TransformedPulseJF[0].SetIm(0);
 
            QVectorConstView samplesJviewB(samplesJorig,0,fTimeLength);
            const QVector& samplesJB = samplesJviewB.GetVector();
 
            err = Transform(samplesJB, TransformedPulseJB);
            if(err != 0) {
                Error("QFFTRealComplex::TransformToFreq: pulse J");
                return;
            } 
            TransformedPulseJB[0].SetRe(0);
            TransformedPulseJB[0].SetIm(0);
 
            ChannelInfo& chInfo = fRosetta[chanI][chanJ];
            if(chInfo.fNevents == 0) {
                chInfo.fCovarianceFreq.Resize(fTimeLength,fTimeLength);
                chInfo.fCovarianceFreq.Initialize(0,0);
                if(chanI == chanJ) {
                    fChannelList.push_back(chanJ);
                }
            }
            QMatrixC element,phases;
            element = TransformedPulseIF*TransformedPulseJF.H();
            phases = TransformedPulseIB*TransformedPulseJB.H();
            QMatrixC magnitudes(phases.Magnitude());
            phases.Div(magnitudes);
            phases.Conjugate();
            element.Mult(phases);
            chInfo.fCovarianceFreq += element;
            if(chanI == chanJ) {
                if(chInfo.fNevents == 0) {
                    chInfo.fPower.Resize(fTimeLength);
                    chInfo.fPower.Initialize(0);
                    chInfo.fCoherence.Resize(fTimeLength);
                    chInfo.fCoherence.Initialize(0);
                }
                chInfo.fPower += TransformedPulseIF.GetModulusSquare();
                chInfo.fCoherence += TransformedPulseIB.Mult(TransformedPulseIF.Conj());
            }
            chInfo.fNevents++;
        }
    }
 
}
 
void MNoiseCrossFrequencyCovariance::Done()
{
    if(fTransformer != 0) delete fTransformer;
    std::stringstream str;
    str<<"Channel List: ";
    for(size_t i=0;i<fChannelList.size();i++) {      
        str<<fChannelList[i]<<" ";
    }
    str<<std::endl;
    Debug(str.str().c_str());
 
 
    std::map<int,std::map<int,ChannelInfo> >::const_iterator iterI = fRosetta.begin();
    while(iterI != fRosetta.end() ) {
        std::map<int,ChannelInfo>::const_iterator iterJ = iterI->second.begin();
        while(iterJ != iterI->second.end() ) {
            int chI =iterI->first;
            int chJ =iterJ->first;
            const ChannelInfo& chInfo=iterJ->second;
            QMatrixC cov = chInfo.fCovarianceFreq/chInfo.fNevents;
            GlobalHandle<QMatrixC> covHandle("CovarianceFreq");
            
            int dataset = -1;
            if(fValidityKind=="dataset") {
                GlobalHandle<QInt> dHandle("Dataset");
                GlobalData().Get("",&dHandle,"");
                dataset = dHandle.Get();
                covHandle.SetDataset(dataset);
            } else if (fValidityKind=="run") {
                if (fSRuns.size()==1){
                    covHandle.SetRun(fSRuns[0]);
                }
                else Panic("ValidityKind set to run is not permitted when reading multiple runs.");
            }
            else Panic("ValidityKind can be set only to run or dataset. Check cfg.");
            covHandle.SetChannel(chI);
            covHandle.SetChannel2(chJ);
            covHandle.Set(cov);
            GlobalData().Set(&covHandle,fOutput);
 
            if(chI == chJ) {
                QAverageNoiseHandle nfh(chI,"NoiseAvgPS");
                nfh.Set(chInfo.fPower / chInfo.fNevents);
                GlobalHandle<QVectorC> cohHandle("Coherence");
                cohHandle.SetChannel(chI);
                cohHandle.Set(chInfo.fCoherence / chInfo.fNevents);
                if(fValidityKind =="run") {
                    nfh.SetRun(fSRuns[0]);
                    cohHandle.SetRun(fSRuns[0]);
                } else {
                    nfh.SetDataset(dataset);
                    cohHandle.SetDataset(dataset);
                }
                GlobalData().Set(&nfh, fOutput);
                GlobalData().Set(&cohHandle,fOutput);
            }
            iterJ++;
        }
        iterI++;
    }
 
}
 
int MNoiseCrossFrequencyCovariance::Transform(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 fTransformer->TransformToFreq(zeroMeanPulse, out);
}