QPulseBasicParams.cc

Go to the documentation of this file.

#include "QPulseBasicParams.hh"
#include <math.h>
#include "TMath.h"
 
using namespace Diana;
 
QPulseBasicParams::QPulseBasicParams(const QVector& pulse, int triggerPosition, double baseline,double baselineRMS, int ADCMax, int ADCMin, int samplingFreq)
 
{
  fSamplingFreq = samplingFreq;
   fBaseline = baseline;
   fBaselineRMS = baselineRMS;
   fTriggerPosition = triggerPosition;
   FindMaximumPosition(pulse);
   ComputeHeight(pulse);
   FindEnd(pulse, 30,3);
   FindStart(pulse,60);
//   ComputeRaise(pulse,10,90);
//   ComputeDecay(pulse,90,10);
   ComputeDecayAndRise(pulse,7,0.1,0.9,0.3);
   ComputeSlowDecay(pulse,7,0.3,0.1);
   CheckSaturation(pulse, ADCMax, ADCMin);
}
 
QPulseBasicParams::QPulseBasicParams(int triggerPosition, double baseline, double baselineRMS, int samplingFreq)
{
  fSamplingFreq = samplingFreq;
   fBaseline = baseline;
   fBaselineRMS = baselineRMS;
   fTriggerPosition = triggerPosition;
}
 
QError QPulseBasicParams::FindMaximumPosition(const QVector& pulse) 
{
    QError err = QERR_SUCCESS;
    //    size_t window = 10;
    int window = fTriggerPosition/2;
    int size = (int) pulse.Size(); 
    int from = fTriggerPosition-window;
    double baseline = pulse.GetMean(TMath::FloorNint(0.8*fTriggerPosition));
    QVector temp = pulse - baseline;
    if(from < 0){
        return QERR_OUT_OF_RANGE;
    }
    if((from+2*window)>=size){
        return QERR_OUT_OF_RANGE;
    }
    //int currMax = -1, prevMax = -1;
    //double prevMaxValue = -1e20, currMaxValue = -1e20;
 
    /*while(from+window < size){
        currMax = pulse.GetMaxIndex(window,from);
        currMaxValue = pulse[currMax];
        if(currMaxValue < prevMaxValue) break;
        else {
            prevMaxValue = currMaxValue;
            prevMax = currMax;
        }
        from += window;
    }*/
    
    fMaxPos = (double)temp.Abs().GetMaxIndex(2*window,from);
    //fMaxPos =(double) prevMax;
    return err;
}
 
 
QError QPulseBasicParams::ComputeHeight(const QVector& pulse) 
{
    QError err = QERR_SUCCESS;
    if(fMaxPos >= 0. && fMaxPos < pulse.Size())
        fHeight = double(pulse[(size_t)lround(fMaxPos)] - fBaseline);
    else fHeight = Q_DOUBLE_DEFAULT;
    return err;
}
 
QError QPulseBasicParams::FindEnd(const QVector& pulse, size_t meansize, double nBaseRMS)
{
    QError err = QERR_SUCCESS;
    size_t fNPoints = pulse.Size();
    double temp=0;
    if( fNPoints > 1000) {
        meansize = meansize*2;
    }
    size_t bufsize=(fNPoints)/2;
    double* binbuf = (double*) malloc(bufsize*sizeof(double));
 
    //fill buffer with mean of meansize subsequent
    for (size_t i=0; i<bufsize; i++) {
        binbuf[i]=0;
        for (size_t j=0; j<meansize; j++) {
            if ((fNPoints-bufsize+i+j)< fNPoints) {
                temp = pulse[(fNPoints-bufsize)+i+j];
                binbuf[i]+= temp;
            }
            else binbuf[i]+=temp;
        }
        binbuf[i]/=meansize;
    }
 
    // compute fit end
    fEnd = fNPoints - 1;
    for (int i = bufsize-1; i>=0; i--) 
        if ((binbuf[i]-fBaseline) > nBaseRMS*fBaselineRMS) {
            fEnd= (fNPoints-bufsize) + i ;
            break;
        }
    free(binbuf);
    return err;
}
 
QError QPulseBasicParams::FindStart(const QVector& pulse, int bufsize)
{
    double binv;
    QError Err = QERR_SUCCESS;
    double binbuf[bufsize];;
    int maxbin = (int)fMaxPos;
    fStart = 0;
    if (maxbin<bufsize){ //if the maximum is in the pretrigger region
        Err= QERR_UNKNOWN_ERR;
        return Err;
    }
    for(int i = maxbin; i > maxbin - bufsize; i--){
        binv = pulse[i] - fBaseline;  
        binbuf[i - (maxbin - bufsize) - 1] = binv;
    }
    if(fBaselineRMS == 0) fBaselineRMS = 1;
    for(int i = maxbin - bufsize; i > 0; i--){
        for(int j = bufsize - 1; j > 0; j--) {
            binbuf[j] = binbuf[j-1];   
        }
        binbuf[0] = pulse[i] - fBaseline;        
        if( binbuf[bufsize - 1] <= binbuf[0] + 2*fBaselineRMS && binbuf[bufsize - 1] >= binbuf[0] - 2*fBaselineRMS  ){
            fStart = i + bufsize -1;
            Err = QERR_SUCCESS;
            break;
        }
    }
    return Err;
}
 
QError QPulseBasicParams::ComputeRaise(const QVector& pulse, double percLow,double percHigh) 
{
    QError err = QERR_SUCCESS;                                                                          
    int t1 = (int)fStart;
    int t2 = (int)fMaxPos;
    for(int t = int(fStart)+1; t <  fMaxPos; t++) {
        if((pulse[t] - fBaseline) < percLow*fHeight/100) t1 = t;
        if((pulse[t] - fBaseline) < percHigh*fHeight/100) t2 = t;                                              
    }
    fRaise = t2 - t1;        
    return err;
}
 
QError QPulseBasicParams::ComputeDecay(const QVector& pulse, double percHigh,double percLow)
{
    QError err = QERR_SUCCESS;                                                                          
    int t1 = (int)fMaxPos;
    int t2 = (int)fEnd;
    for(int t = int(fMaxPos)+1; t <  fEnd; t++) {
        if((pulse[t] - fBaseline) >  percHigh*fHeight/100) t1 = t;
        if((pulse[t] - fBaseline) > percLow*fHeight/100) t2 = t;                                              
    }
    fDecay = t2 - t1;        
    return err;
}
 
QError QPulseBasicParams::CheckSaturation(const QVector& pulse, int ADCMax, int ADCMin) 
{
        QError err = QERR_SUCCESS;
        
//      double ADCMax = 65534.5;
//      double ADCMin = 0.5;
        
        ADCMax -= 1;
        ADCMin += 1;
        
        fIsSaturatedHigh=false;
        fIsSaturatedLow=false;
        double max = pulse.GetMax();
        double min = pulse.GetMin();
        
        if(max > ADCMax) fIsSaturatedHigh=true;
        if(min < ADCMin) fIsSaturatedLow=true;
        return err;
}
 
QError QPulseBasicParams::ComputeDecayAndRise(const QVector& pulse, int smoothingSize, float fracMin, float fracMax,float fracMinDecay)
{
    QError err = QERR_SUCCESS;
    double RT, DT;
    QVector smoothed = pulse;
 
    int size = smoothed.Size(); 
 
    int num = smoothingSize;
    for(int i=0;i<size;i++)
    {
        int n=num;
        smoothed[i]*=num;
        for(int k=1;k<=(num-1)/2;k++) {
            if((i-k)>=0)
            {
                smoothed[i]+=pulse[i-k]*(num-2*k);
                n+=num-2*k;
            }
            if((i+k)<size)
            {
                smoothed[i]+=pulse[i+k]*(num-2*k);
                n+=num-2*k;
            }
        }
        smoothed[i]/=n;
        smoothed[i]-=fBaseline;
    }
 
    // Rise
    double max = fHeight;
    int max_pos=(int) fMaxPos;
    int n = 0;
    double x1 = 0,x0 = 0;
    bool times1 = false;
    RT = 0;
    if(max_pos>=0.) {
        for(n = max_pos; n >= 0; n--) {
            if(smoothed[n]<max*fracMax &&times1==false) {
                x1=n+(max*fracMax-smoothed[n])/(smoothed[n+1]-smoothed[n]);
                times1=true;
            }
            if(smoothed[n]<max*fracMin) {
                x0=n+(max*fracMin-smoothed[n])/(smoothed[n+1]-smoothed[n]);
                break;
            }
        }
        if(n!=0) RT=(x1-x0); 
    }
 
 
    // Decay
    times1=false; 
    x1=x0=0;
 
    DT = 0;
    if(max_pos != size-1 && max_pos >= 0.) {
        for(n=max_pos;n<size;n++) {
            if(smoothed[n]<max*fracMax &&times1==false) {
                x1=n-1+((max*fracMax-smoothed[n-1])/(smoothed[n]-smoothed[n-1]));
                times1=true;
            }
            if(smoothed[n]<max*fracMinDecay &&times1==true) {
                x0=n-1+((max*fracMinDecay-smoothed[n-1])/(smoothed[n]-smoothed[n-1]));
                break;
            }
        }
        if(n!=size) DT=(x0-x1);
    }
    fDecay = DT;
    fRaise = RT;
    return err;
}
QError QPulseBasicParams::ComputeSlowDecay(const QVector& pulse, int smoothingSize, float fracMax,float fracMinDecay)
{
    QError err = QERR_SUCCESS;
    double DT;
    QVector smoothed = pulse;
 
    int size = smoothed.Size(); 
 
    int num = smoothingSize;
    for(int i=0;i<size;i++)
    {
        int n=num;
        smoothed[i]*=num;
        for(int k=1;k<=(num-1)/2;k++) {
            if((i-k)>=0)
            {
                smoothed[i]+=pulse[i-k]*(num-2*k);
                n+=num-2*k;
            }
            if((i+k)<size)
            {
                smoothed[i]+=pulse[i+k]*(num-2*k);
                n+=num-2*k;
            }
        }
        smoothed[i]/=n;
        smoothed[i]-=fBaseline;
    }
 
    // Decay
    double max = fHeight;
    int max_pos=(int) fMaxPos;
    int n = 0;
    double x1 = 0,x0 = 0;
    bool times1 = false;
 
    DT = 0;
    if(max_pos != size-1 && max_pos >= 0.) {
        for(n=max_pos;n<size;n++) {
            if(smoothed[n]<max*fracMax &&times1==false) {
                x1=n-1+((max*fracMax-smoothed[n-1])/(smoothed[n]-smoothed[n-1]));
                times1=true;
            }
            if(smoothed[n]<max*fracMinDecay &&times1==true) {
                x0=n-1+((max*fracMinDecay-smoothed[n-1])/(smoothed[n]-smoothed[n-1]));
                break;
            }
        }
        if(n!=size) DT=(x0-x1);
    }
    fSlowDecay = DT;
    return err;
}