MNoiseCrossPowerSpectrum.cc

Go to the documentation of this file.

#include "MNoiseCrossPowerSpectrum.hh"
#include "QEvent.hh"
#include "QEventList.hh"
#include "QRealComplexFFT.hh"
#include "QCountPulsesData.hh"
#include "QBaselineData.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 "QStdVector.hh"
#include <iostream>
 
using namespace Diana;
 
REGISTER_MODULE(MNoiseCrossPowerSpectrum)
 
 
void MNoiseCrossPowerSpectrum::Init(QEvent& ev)
{
    // called before event loop
    fTimeLength = 0;
    fTransformer = 0;
    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);
 
    std::vector<std::string> pulseLabels;
 
    pulseLabels.push_back("DAQ@Pulse");
    fPulseLabels = GetVectorString("PulseLabels",pulseLabels,false);
    if(fPulseLabels.empty()) Panic("PulseLabels is empty");
 
    ev.Require<QHeader>("DAQ","Header");
    ev.Require<QPulseInfo>("DAQ","PulseInfo");
    for(size_t i = 0; i < fPulseLabels.size(); i++) {
        ev.RequireByLabel<QPulse>(fPulseLabels[i]);
    }
 
 
    fSRuns.clear();
    fCurrentRun = -1;
 
}
 
void MNoiseCrossPowerSpectrum::Do(QEvent& ev, const QEventList& neighbours)
{
    // called event by event
    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>(fPulseLabels[0]).GetSamples().Size();
        fTransformer = new QRealComplexFFT(fTimeLength);
        fTransformer->SetWindowType(fWindowType, fCoherentGain);
        fCov.SetNFreq(fTimeLength);
    }
 
    QVectorC TransformedPulseI(fTimeLength);
    QVectorC TransformedPulseJ(fTimeLength);
    QVector zeroMeanPulse(fTimeLength), meanPulse(fTimeLength);
    for(size_t i = 0;  i < allchans.Size(); i++) {
        int chanI = allchans[i].Get<QPulseInfo>("DAQ","PulseInfo").GetChannelId();
        for(size_t li = 0; li < fPulseLabels.size(); li++) {
            if(li>0) chanI += 1000;
            const QVector& samplesI =  allchans[i].GetByLabel<QPulse>(fPulseLabels[li]).GetSamples();
            double mean = samplesI.Sum(samplesI.Size())/samplesI.Size();
            meanPulse.Initialize(mean);
            zeroMeanPulse = samplesI;
            zeroMeanPulse -= meanPulse;
            int err = fTransformer->TransformToFreq(zeroMeanPulse, TransformedPulseI);
            if(err != 0) {
                Error("QFFTRealComplex::TransformToFreq: pulse I");
                return;
            } 
            TransformedPulseI[0].SetRe(0);
            TransformedPulseI[0].SetIm(0);
            for(size_t j = i ;  j < allchans.Size(); j++) {
                int chanJ = allchans[j].Get<QPulseInfo>("DAQ","PulseInfo").GetChannelId();
                for(size_t lj = 0; lj < fPulseLabels.size(); lj++) {
                    if(lj>0) chanJ += 1000;
                    const QVector& samplesJ =  allchans[j].GetByLabel<QPulse>(fPulseLabels[lj]).GetSamples();
                    mean = samplesJ.Sum(samplesJ.Size())/samplesJ.Size();
                    meanPulse.Initialize(mean);
                    zeroMeanPulse = samplesJ;
                    zeroMeanPulse -= meanPulse;
 
                    int err = fTransformer->TransformToFreq(zeroMeanPulse, TransformedPulseJ);
                    if(err != 0) {
                        Error("QFFTRealComplex::TransformToFreq: pulse J");
                        return;
                    } 
                    TransformedPulseJ[0].SetRe(0);
                    TransformedPulseJ[0].SetIm(0);
 
                    ChannelInfo& chInfo = fRosetta[chanI][chanJ];
                    if(chInfo.fCovariance.Size() == 0) {
                        chInfo.fNevents = 0;
                        chInfo.fCovariance.Resize(fTimeLength);
                        chInfo.fCovariance.Initialize(0,0);
                        if(chanI == chanJ) {
                            fChannelList.push_back(chanJ);
                        }
                    }
                    ChannelInfo& chInfoT = fRosetta[chanJ][chanI];
                    if(chInfoT.fCovariance.Size() == 0) {
                        chInfoT.fNevents = 0;
                        chInfoT.fCovariance.Resize(fTimeLength);
                        chInfoT.fCovariance.Initialize(0,0);
                    }
                    QVectorC element = TransformedPulseI.Mult(TransformedPulseJ.Conj());
                    chInfo.fCovariance += element;
                    chInfo.fNevents++;
 
                    chInfoT.fCovariance += element.Conj();
                    chInfoT.fNevents++;
                }
            }
        }
    }
 
}
 
void MNoiseCrossPowerSpectrum::Done()
{
    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());
 
    fCov.SetChannels(fChannelList);
 
 
    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;
            QVectorC vec = chInfo.fCovariance/chInfo.fNevents;
            fCov.SetCovariance(chI,chJ,vec);
            if(chI == chJ) {
                fCov.SetNEvents(chI,chJ,chInfo.fNevents/2);
            } else {
                fCov.SetNEvents(chI,chJ,chInfo.fNevents);
            }
            iterJ++;
        }
        iterI++;
    }
 
    GlobalHandle<QChannelCovariance> CovSaver("Covariance");
    CovSaver.Set(fCov);
    GlobalHandle<QStdVector<std::string> > LabelsSaver("PulseLabels");
    LabelsSaver.Set(fPulseLabels);
    GlobalHandle<QInt> StepSaver("PulseStep");
    StepSaver.Set(1000);
    if(fValidityKind=="dataset"){
        GlobalHandle<QInt> dHandle("Dataset");
        GlobalData().Get("",&dHandle,"");
        int dataset = dHandle.Get();
        CovSaver.SetDataset(dataset);
        LabelsSaver.SetDataset(dataset);
        StepSaver.SetDataset(dataset);
    }
    else if (fValidityKind=="run"){
        if (fSRuns.size()==1){
            CovSaver.SetRun(fSRuns[0]);
            CovSaver.SetRun(fSRuns[0]);
            LabelsSaver.SetRun(fSRuns[0]);
            StepSaver.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.");
 
 
 
 
    GlobalData().Set(&CovSaver,fOutput);
    GlobalData().Set(&StepSaver,fOutput);
    GlobalData().Set(&LabelsSaver,fOutput);
    if(fTransformer != 0) delete fTransformer;
}