QTriggerOptimumFilter.cc

Go to the documentation of this file.

#include "QTriggerOptimumFilter.hh"
#include "TMath.h"
#include <vector>
#include <numeric>
#include "QGlobalDataManager.hh"
#include "QAveragePulseHandle.hh"
#include "QAverageNoiseHandle.hh"
 
 
ClassImp(Diana::QTriggerOptimumFilter);
using std::endl;
using std::cout;
using std::cerr;
using Diana::QVector;
using Diana::QVectorC;
 
Q_BEGIN_NAMESPACE
 
void QTriggerOptimumFilter::GetParametersFromCFG(){
    if(reader==0){
        cout<<"NO READER PASSED TO TRIGGER"<<endl;
        exit(-1);
    }
    SetCFGParameters(reader->GetString(("AveragePulse"+triggerLabel),""),reader->GetString(("AverageNoisePS"+triggerLabel),""),reader->GetInt(("Dataset"+triggerLabel),1,true), reader->GetDouble(("ThresholdInRMS"+triggerLabel).c_str(),5.),reader->GetDouble(("ReceedLevel"+triggerLabel).c_str(),2.),reader->GetString(("AvgPulseOwner"+triggerLabel),"AveragePulses",false),reader->GetString(("AvgNoiseOwner"+triggerLabel),"NoiseAvgPowerSpectrum",false));
}
 
const QVector& QTriggerOptimumFilter::Filter(const QVector& window_stream){
    QVectorC stream_FD;
    filt_stream_window.Clear();
    if((int)window_stream.Size()!=2*windowsize){
        cout<<"FATAL ERROR: Stream not the right size"<<endl;
        cout<<"Size should be "<<2*windowsize<<" Found instead "<<window_stream.Size()<<endl;
    }
    fTransformer_d->TransformToFreq(window_stream,stream_FD);
    stream_FD[0].Set(0,0);
    fTransformer_d->TransformFromFreq(fFilterFD_d.Mult(stream_FD),filt_stream_window);
    return filt_stream_window;
}
 
bool QTriggerOptimumFilter::BuildOptimumFilter(){
    
    Diana::QGlobalDataManager dm;
    
    // Get Average pulse
    /*string apowner = "AveragePulses";
    string anowner = "NoiseAvgPowerSpectrum";
    if(extraLabel.compare("")!=0){
        apowner += "_";
        apowner += extraLabel;
        anowner += "_";
        anowner += extraLabel;
    }*/
 
    QAveragePulseHandle aph (chan);
    aph.SetDataset(dataset);
    QError err = dm.Get(apowner,&aph,average_pulse);
    if (!aph.IsValid()){
        cout<<"Impossible to read average pulse"<<endl;
        cout<<err.ToString()<<"\t"<<err.GetDescription()<<endl;
        return false;
    }
        
    ap = aph.Get();
    ap = ap - ap[0];
    ap_max_norm = ap.GetMax();
    ap/=ap.GetMax();
    //int apmaxIndex = ap.GetMaxIndex();
    
    QAverageNoiseHandle anh(chan);
    anh.SetDataset(dataset);
    err = dm.Get(anowner,&anh,average_noise);
    if (!anh.IsValid()){
        cout<<"Impossible to read average noise"<<endl;
        cout<<err.ToString()<<"\t"<<err.GetDescription()<<endl;
        return false;
    }
    an = anh.Get();
    an[0]=0;
    
    fTransformer = new QRealComplexFFT(ap.Size());
    fTransformer_d= new QRealComplexFFT(2*ap.Size());
    
    smoothing_ap.Resize(ap.Size());
    size_t i = 1;
    size_t transient =0;
    while(transient <= 0){
        transient =  (2*i*ap.Size())/25;
        i++;
    }
    for(size_t i= 0; i<smoothing_ap.Size();i++){
        if(2*i<=transient){
            smoothing_ap[i] = TMath::Cos(2*TMath::Pi()*i/transient+TMath::Pi()) * 0.5 +0.5;
        }else if(2*i>=(2*smoothing_ap.Size() - transient)){
            smoothing_ap[i] = TMath::Cos(2*TMath::Pi()*i/transient) * 0.5 +0.5;
        }else{
            smoothing_ap[i] = 1;
        }
    }
    
    smoothed_ap = ap.Mult(smoothing_ap);
    
    
    fTransformer->TransformToFreq(smoothed_ap,ap_FD);
    
    QVector SNR = ap_FD.GetModulusSquare().Div(an);
    double M = 1./SNR.Sum(ap.Size()-1,1);
    
    fFilter1 = M*ap_FD.Conj().Div(an);
    fFilter1 *= ap.Size();
    fFilter1[0].Set(0,0);
    
    
    fTransformer->TransformFromFreq(fFilter1,fFilterTD);
    
    fFilterTD_d.Resize(2*fFilterTD.Size());
    fFilterTD_d.Initialize(0);
    QVector filtertmp = fFilterTD.Mult(smoothing_ap);
    for(size_t i = 0; i<fFilterTD.Size();i++){
        filtertmp.Append(0);
    }
    fFilterTD_d = filtertmp.Shift((3*fFilterTD.Size())/2);
    fTransformer_d->TransformToFreq(fFilterTD_d,fFilterFD_d);
    
    rms = TMath::Sqrt(M);
    windowsize = (int)ap.Size();
    debounce = windowsize;
    
    
    ap_d.Resize(2*ap.Size());
    ap_d.Initialize(0);
    for(size_t i = 0; i<ap.Size();i++){
        ap_d[i+ap.Size()/2]=ap[i];
    }
    
    
    apfilt = Filter(ap_d);
    int max_index = ap_d.GetMaxIndex();
    double threshold = ap_d[max_index]/receed_lvl;
    int threshold_index = max_index/2;
    i = (size_t)max_index;
    while( (ap_d[i]-threshold)>=0 and (i>=0)){
        threshold_index=i;
        i--;
    }
    if((threshold - ap_d[threshold_index])>(ap_d[threshold_index+1]-threshold)){
        threshold_index+=1;
    }
    threshold_index = max_index - threshold_index;
    IndexDifference = (ap.Size() -max_index +threshold_index);
    optimum_filter_built=true;
    return true;
}
 
void QTriggerOptimumFilter::TriggerStream(const vector<double>& stream){
    if (!optimum_filter_built){
        //reader->Debug("Building OF");
        cout<<"Building OF"<<endl;
        if(!read_parameters){
            GetParametersFromCFG();
        }
        optimum_filter_built = BuildOptimumFilter();
    }
    //empty triggers
    triggers.clear();
    triggerHeights.clear();
    filteredStream.clear();
    
    if(samplingFrequency ==0 or windowsize==0)cout<<"Sampling frequecy "<<samplingFrequency<<"Hz"<< "\t Window Size = "<<windowsize<<endl;
    
    
 
    QVector filtered_win;
    QVector win(2 * windowsize);
    win.Initialize(0);
    for(size_t i = 0; i<(size_t)windowsize and i<stream.size() ;i++){
        win[i] = stream[i]-stream[0];
    }
    
    win=win.Shift(windowsize/2);
    //if((int)stream.size()<2*windowsize){
    //}
    //for(size_t i = 0; i<(size_t)win.Size();i++){
    //    cout<<win[i]<<",";
    //}
    //cout<<endl<<endl;
    filtered_win = Filter(win);
    //for(size_t i = 0 ; i<filtered_win.Size();i++){ 
    //    cout<<filtered_win[i]<<",";
    //}
    //cout<<endl<<endl;
//    for(size_t i = 0 ; i<(size_t)((double)windowsize*1.5);i++){
    for(size_t i = 0; i<(size_t)windowsize;i++){
        filteredStream.push_back(filtered_win[(windowsize/2)+i]);
    }
 
    for(size_t i =(size_t)(windowsize*0.5);i<stream.size();i+=(size_t)(windowsize)){
        if(i+2*windowsize<stream.size()){
            win = vector<double> (stream.begin()+i, stream.begin()+i+2*windowsize);
                win = win - stream[i];
        }else{
            win.Resize(2*windowsize);
            win.Initialize(0);
            for(size_t k = 0; k<stream.size()-i;k++){
                win[k]=stream[i+k]-stream[i];
            }
            win=win.Shift(windowsize/2);
        }
        
        filtered_win = Filter(win);
        
        for(size_t j = 0;j<(size_t)windowsize;j++){
            filteredStream.push_back(filtered_win[j+(size_t)((double)windowsize*0.5)]);
        }
    }
    
    
    filteredStream.resize(stream.size());
    
    vector<int> sample_idx(filteredStream.size());
    vector<bool> enabled(filteredStream.size(),true);
    vector<int> mapping(filteredStream.size(),true);
    
    std::iota(sample_idx.begin(),sample_idx.end(),0); //Initializing
    sort( sample_idx.begin(),sample_idx.end(), [&](int i,int j){return filteredStream[i]>filteredStream[j];} );
    
    std::iota(mapping.begin(),mapping.end(),0); //Initializing
    sort(mapping.begin(),mapping.end(), [&](int i,int j){return sample_idx[i]<sample_idx[j];});
    
    for(size_t i = 0; i<enabled.size();i++){
        while(i <enabled.size()){
            if(enabled[i]==true) break;
            i++;
        }
        
        if(i==enabled.size()) break;
        if(filteredStream[sample_idx[i]] <= rms*triggerLevel) break;
        enabled[i] = false;
        
        if((sample_idx[i]-IndexDifference)>=0){
            triggers.push_back(sample_idx[i]-IndexDifference);
            triggerHeights.push_back(filteredStream[sample_idx[i]]);
        }else{
            triggers.push_back(sample_idx[i]);
            triggerHeights.push_back(filteredStream[sample_idx[i]]);
        }
        int min_deb = debounce*4/5;
        int max_deb = debounce/5;
        while(sample_idx[i]-min_deb<0){
            min_deb--;
        }
        while(sample_idx[i]+max_deb>=(int)enabled.size()){
            max_deb--;
        }
        for(int j = -min_deb; j<max_deb;j++){
            enabled[mapping[sample_idx[i]+j]] = false;
        }
        
    }
}
 
 
 
Q_END_NAMESPACE