QTriggerBULLDAQ.cc

Go to the documentation of this file.

#include "QTriggerBULLDAQ.hh"
#include "QReader.hh"
#include "TMath.h"
#include <vector>
#include "QVector.hh"
#include <set>
 
ClassImp(Diana::QTriggerBULLDAQ);
using std::endl;
using std::cout;
using std::cerr;
using Diana::QVector;
 
Q_BEGIN_NAMESPACE
 
void QTriggerBULLDAQ::GetParametersFromCFG(){
    if(reader==0){
        cout<<"NO READER PASSED TO TRIGGER"<<endl;
        exit(-1);
    }
    
    SetCFGParameters(reader->GetDouble(("ThresholdInRMS"+triggerLabel).c_str(),5.),reader->GetDouble(("LowPassFrequency"+triggerLabel).c_str(),50));
}
 
 
 
void QTriggerBULLDAQ::TriggerStream(const vector<double>& stream){
    //if parameters not read from cfg
    if (!read_parameters){
            GetParametersFromCFG();
    }
    
    
    //empty triggers
    triggers.clear();
    triggerHeights.clear();
    filteredStream.clear();
    
    if(samplingFrequency ==0 or windowsize==0)cout<<"Sampling frequecy " << samplingFrequency <<"\tWindow Size = "<<windowsize<<endl;
    
    
    int npoints = TMath::FloorNint(TMath::Sqrt(0.196202+(fCutoff/samplingFrequency)*(fCutoff/samplingFrequency))/(fCutoff/samplingFrequency)); //dsp.stackexchange.com/questions/9966/what-is-the-cut-off-frequency-of-a-moving-average-filter 
    if(npoints<1) npoints=1;
    int triggerDelay = 2;//npoints-1;
    
    
    QVector pulse (stream.size());
    
    //Differentiate (high pass filter)
    for(size_t i = stream.size()-1; i > 0; i--){
        pulse[i] =(double)( stream[i] - stream[i-1] );    
    }
    pulse[0] = pulse[1];
    
    
    //Integrate (pulse height reconstruction)
    for(size_t i =0; i<pulse.Size();i++){
        if(i<pulse.Size()-npoints){
            pulse[i] = pulse.Sum(npoints,i);
        }else {
            pulse[i]=pulse[i-1];
        }
    }
    
    
    //Moving Average (low pass filter) // made on reverse to remove time shift on filters
    for(int i = (int)pulse.Size()-1; i>-1;i--){
        if(i>=npoints){
            pulse[i] = pulse.GetMean(npoints,i-npoints);
        }else{
            pulse[i]=pulse[i+1];
        }
    }
    
     //Find RMS
     if(rms ==-1){
        //int step = TMath::FloorNint((double)windowsize/10);//chosen because of dsp.stackexchange.com/questions/24468/derivative-filter-transfer-function
        int step = TMath::FloorNint(samplingFrequency/10);//chosen because of dsp.stackexchange.com/questions/24468/derivative-filter-transfer-function
            if((size_t)(4*step)>pulse.Size()){
                step = TMath::FloorNint((double)windowsize/5);//chosen because of dsp.stackexchange.com/questions/24468/derivative-filter-transfer-function
        
            }
        double rms_new=-1;   
        std::set<double> rmsss;
        for(size_t i =0; i<pulse.Size()-step;i+=step){
            if(rms_new ==-1){
                rms_new = pulse.GetSubVector(step,i).GetRMS();
                rms = rms_new +1;
            }else{
                rms_new = pulse.GetSubVector(step,i).GetRMS();
            }
            if(rms_new<=rms and rms_new>0){
                rms=rms_new;
            }
            if (rms_new>0)rmsss.insert(rms_new);
        }
        std::set<double>::const_iterator it=rmsss.cbegin();
        std::advance(it,rmsss.size()/2);
        double median=*it;
        if (rmsss.size()%2==0){
            median = (median + *std::prev(it))/2;
        }
    }
    for (size_t i =0 ; i < pulse.Size();i++){
        filteredStream.push_back(pulse[i]);    
    }
    
    size_t first_trig=pulse.Size()+1;
    size_t prev_trig=first_trig;
    double height =0;
    size_t sample;
    for(size_t i=0; i<pulse.Size();i++){
        if(pulse[i]>=triggerLevel*rms){//If sample triggers
            if(first_trig==pulse.Size()+1){//first trigger point in chunk
                first_trig = i;
                prev_trig = i;
            }else if (i==prev_trig+1){//same peak
                prev_trig=i;
            }else if(i>prev_trig+1){//new peak
                first_trig=i;
                prev_trig=i;
            }
        }
        if((i==prev_trig+1 and pulse[i]<triggerLevel*rms) or (i == pulse.Size()-1 and i ==prev_trig)){//as soon as the peak ends or if the peak is ongoing and we are at the end of the stream
            
            if(triggers.size()==0){
                for(size_t j = first_trig; j<=prev_trig;j++){
                    if(j==first_trig){
                        height=pulse[j];
                        sample=j;
                    }else{
                        if(height<pulse[j]){
                            height=pulse[j];
                            sample=j;
                        }
                    }
                }
                
                triggers.push_back((int)sample-triggerDelay);
                triggerHeights.push_back(height);
            }else{
                if((double)((int)first_trig-triggers[triggers.size()-1]) >= (double)debounce){//don't trigger if trigger is too near the previous one
                    for(size_t j =first_trig; j<=prev_trig;j++){
                        if(j==first_trig){
                            height=pulse[j];
                            sample=j;
                        }else{
                            if(height<pulse[j]){
                                height=pulse[j];
                                sample=j;
                            }
                        }
                    }
                    triggers.push_back((int)sample-triggerDelay);
                    triggerHeights.push_back(height);
                }
            }
        }
    }
}
 
 
 
Q_END_NAMESPACE