MPeakDetector.cc

Go to the documentation of this file.

#include "MPeakDetector.hh"
#include "QEvent.hh"
#include "QBaseType.hh"
#include "QVector.hh"
#include "QVectorI.hh"
#include "QPulse.hh"
#include "QPulseInfo.hh"
#include "QStdVector.hh"
#include "QCountPulsesData.hh"
#include <algorithm>
#include "QFir.hh"
#include <cmath>
 
 
REGISTER_MODULE(MPeakDetector)
 
using namespace Diana;
using namespace std;
 
void MPeakDetector::Init(QEvent& ev)
{
    fPulseLabel = GetString("PulseLabel","DAQ@Pulse",false);
    fThreshold = GetDouble("ThresholdNumberOfSigma", 5);
    if(fThreshold <= 2) Panic("ThresholdNumberOfSigma must be > 2");
    fSubSampling = GetInt("SubSampling",5,false);
    if(fSubSampling < 2) Panic("SubSampling must be > 1");
 
    ev.RequireByLabel<QPulse>(fPulseLabel);
    ev.RequireByLabel<QPulseInfo>("DAQ@PulseInfo");
 
    ev.Add<QCountPulsesData>("CountPulsesData");
 
    if(fSubSampling > 1) {
        QCFirData configuration;
        configuration.fCutoff = 2./(fSubSampling);
        configuration.fAttDB = 40;
        configuration.fM = configuration.fAttDB*0.025/configuration.fCutoff;
        if(configuration.fM%2 == 0) configuration.fM++;
        configuration.fMethod = QCFirData::KayserBessel;
        fFir = new QFir(configuration);
 
    } else {
        fFir = 0;
    }
}
 
void MPeakDetector::Do(QEvent& ev)
{
    const QVector&   samples  = ev.GetByLabel<QPulse>(fPulseLabel).GetSamples();
 
    // number of samples to evaluate the RMS
    const int minSamples = int(ev.Get<QPulseInfo>("DAQ","PulseInfo").GetMasterSample().GetSampleIndex()*3./4);
 
    // find peaks in window
    vector<int> peakPositions = GetPeaks(samples,minSamples,1);
    // find peaks in window subsampled by a factor fSubSampling
    vector<int> peakPositionsSub = GetPeaks(samples,minSamples,fSubSampling);
   
    // merge peaks into save vector
    std::vector<int> save = peakPositions;
    for(size_t p = 0; p < peakPositions.size(); p++) {
        int peakPos = peakPositions[p];
        for(vector<int>::iterator pr = peakPositionsSub.begin(); pr != peakPositionsSub.end(); ) {
            int peakPosSub = *pr;
            int delta = abs(peakPosSub - peakPos);
            if(delta <= 3*fSubSampling) {
                pr = peakPositionsSub.erase(pr);
            } else {
                pr++;
            }
        }
    }
    for(size_t pr = 0; pr < peakPositionsSub.size(); pr++) {
        save.push_back(peakPositionsSub[pr]);
    }
    std::sort(save.begin(),save.end());
 
    // save results to event
    ev.Get<QCountPulsesData>("CountPulsesData").fNumberOfPulses = save.size();
    QVector peakpos(save.size());
    for(size_t i =0; i < save.size(); i++) peakpos[i] = save[i];
    ev.Get<QCountPulsesData>("CountPulsesData").fTimeIntervals = peakpos;
 
}
 
void MPeakDetector::Done()
{
}
 
vector<int> MPeakDetector::GetPeaks(const QVector& samples, int minSamples, const size_t red)
{
    QVector samplesReduced;
    const int numSamples = samples.Size();
    int offset = 0;
 
    if(red > 1) {
        offset = fFir->GetFilterReduction();
        QVector filtered;
        fFir->Filter(samples,filtered);
        samplesReduced.Resize(samples.Size()/red-offset);
        samplesReduced.Initialize(0);
        size_t j = 0;
        for(int i = 0; i < numSamples; i++) {
            if( i % red == 0 && j < samplesReduced.Size()) {
                samplesReduced[j]=filtered[i];
                j++;
            }
        }
    } else {
        samplesReduced = samples;
    }
 
    minSamples /= red;
 
    QVector derivative = samplesReduced.Derivative3P();
 
//    derivative[0]=derivative[derivative.Size()-1]=0.;
//    derivative[1]=derivative[derivative.Size()-2]=0.;
    size_t i = 0;
    double rmsMin = 1e20;
    while(i < derivative.Size()-minSamples) {
        double rms = derivative.GetRMS(minSamples,i);
        if(rms < rmsMin) rmsMin = rms;
        i += minSamples/2;
    }
    double thrrms = fThreshold*rmsMin;
 
    vector<int> peakPositions;
    i = 0;
    while ( i < derivative.Size() ) {
        double prevDerivative = derivative[i];
        int sign = 1;
        if(prevDerivative < 0) sign = -1;
 
        int maxPos = i;
        if(prevDerivative*sign > thrrms) {
            i++;
            // look for maximum
            while(i < derivative.Size()) {
                if(derivative[i]*sign < -2*rmsMin) {
                    maxPos = i;
                    break;
                }
                i++;
            }
            // come back to baseline
            while (i < derivative.Size()) {
                if(derivative[i]*sign > (fThreshold-1)*rmsMin) {
                    break;
                }
                i++;
            }
            peakPositions.push_back(maxPos*red+offset);
        } else {
            i++;
        }
 
    }
 
    return peakPositions;
 
}