source: fact/tools/rootmacros/discriminator.C@ 17428

#include <stdio.h>
#include <iostream>

using namespace std;

// this function reads a (hopefully smoothed, and baseline corrected) DRS pipeline and produces 
// The DRs pipiline is given as a vector<float>. 
//
// the Output is a vector of Over-Threshold-Regions, which are described by 
//      * the index of its start
//      * the index of its Maximum
//      * the index of its end
//      * the height of its Maximum
//
// it needs a threshold and the length of the falling edge 

#include "discriminator.h"

vector<Region> * discriminator(
        vector<float>& input,           // vector of floats, the discriminator acts on
        float thr,                                      // threshold
        bool clean,                                     // choose here, if the discriminator should already discard Regions, which are assumend to sit on the predessesors falling edge
        int fallingEdge,                        // number of slices, after the maximum, which should be discarded...
        bool debug
){ //----- start of function discrimnator--------

        // this vector will contain the Region, found by the simple 
        // discriminator 
        vector<Region> * uncleaned = new vector<Region>;

        // The former vector will be cleaned, by checking the distance of the 
        // position of the Maxima .... Maxima, sitting on 
        // falling edges of other pulses, should be discarded.
        vector<Region> * Cleaned = new vector<Region>;
        
        Region currentRegion;
        currentRegion.begin = 0;
        currentRegion.end = 0;
        currentRegion.maxPos = 0;
        currentRegion.maxVal = 0.;
        bool over_thr_found = false;

        for ( unsigned int sl = 0; sl < input.size() ; sl++ ){
                if ( input[sl] > thr ) {
                        // mark as -start-
                        if ( !over_thr_found ) {
                                over_thr_found = true;
                                currentRegion.begin = sl;
                        }

                        // Find the -Maximum- in this Region
                        if ( input[sl] > currentRegion.maxVal) {
                                currentRegion.maxVal = input[sl];
                                currentRegion.maxPos = sl;
                        }

                } else if ( input[sl] < thr ) {

                        // mark the -end- 
                        if ( over_thr_found ) {
                                over_thr_found = false;
                                currentRegion.end = sl;
                                // store the region in a vector.        
                                uncleaned->push_back(currentRegion);
                                currentRegion.begin = 0;
                                currentRegion.end = 0;
                                currentRegion.maxPos = 0;
                                currentRegion.maxVal = 0.;
                        }
                }               
        } // end of for loop over all slices

                if (debug){  // output the vector
                        cout << "------------Discriminator Debug Output:"
                                        <<"----------------" << endl;
                        cout << "\t uncleaned Over Threshold Regions found: "
                                        << uncleaned->size() << endl;
                        for (unsigned int p=0 ; p < uncleaned->size() ; p++ ){
                                cout << p << ":\t"; 
                                cout << uncleaned->at(p).begin << "\t";
                                cout << uncleaned->at(p).end << "\t";
                                cout << uncleaned->at(p).maxPos << "\t";
                                cout << uncleaned->at(p).maxVal << endl;
                        }
                }

        if (clean){
                cleanRegionsOnFallingEdge ( *Cleaned, *uncleaned, fallingEdge, debug);
                delete uncleaned;
                return Cleaned;
        } else
                return uncleaned;




} // end of function - discriminator
///////////////////////////////////////////////////////////////////////////////////////




////////// Cleaning Funcs ////////////

int     cleanRegionsOnFallingEdge ( 
                vector<Region> &dest,
                vector<Region> &src,
                int fallingEdgeLen, 
                bool debug
){
        // This function scans the src-vector for Maxima, sitting on its predesessors falling edge.
        // the out vector is cleaned from these regions. 
        // note: the cleaned regions are appended to the dest-vector.
        //
        //
        // note further: in case pulses pile up to form a broad Region over threshold
        // the maximum of this region is caused by pile up, even though is is far away from the next region.
        // TODO
        // write another cleaning function - or implement into this function a check, 
        // if the maximum the first of all local maxima in the region.
        // maybe the search for the maximum should be oursourced and not be done in the first 
        // discrimnator loop in line 51ff...

        if (dest.size() > 0 && debug)
                cout << "discriminator::cleanRegionsOnFallingEdge: destination vector.size()=" << dest.size() << endl;

        // if nothing to do
        if (src.size() == 0){
                if (debug)
                        cout << "discriminator::cleanRegionsOnFallingEdge: src vector empty" << endl;
                return -1;                      
        }

        // local copy of source
        vector<Region> localSrc(src);
        // last Region in (local)-src
        Region last;
        while (!localSrc.empty()){
                last = localSrc.back();
                localSrc.pop_back();
                                //              last.begin = result->back().begin;
                                //              last.end = result->back().end;
                                //              last.maxPos = result->back().maxPos;
                                //              last.maxVal = result->back().maxVal;    

                // if localSrc is _now_ empty, then there was no predecessor
                if (localSrc.empty()){
                        dest.push_back(last);
                        break;
                }
                if (last.maxPos - localSrc.back().maxPos > fallingEdgeLen)
                        dest.push_back(last);
        }

                if (debug){  // output the vector
                        cout << "------------Discriminator::CleanRegionsOnFallingEdge -  Debug Output:----------------" << endl; 
                        cout << "\t Cleaned Regions found: " << dest.size() << endl; 
                        for (unsigned int p=0 ; p < dest.size() ; p++ ){
                                cout << p << ":\t"; 
                                cout << dest.at(p).begin << "\t";
                                cout << dest.at(p).end << "\t";
                                cout << dest.at(p).maxPos << "\t";
                                cout << dest.at(p).maxVal << endl;
                        }
                }

return 0;
}
 
//////  alternative loop, but it didn't work ////////////////
/*
        vector<vector<DiscOut>::iterator> toBeDeleted;

        for (vector<DiscOut>::iterator it=result->begin(); it != result->end()-1; ++it ){
                if ( (*it).begin==0 ){
                                toBeDeleted.push_back(it);
                }
                if ( ((*(it+1)).maxPos - (*it).maxPos) < fallingEdge ){
                        toBeDeleted.push_back(it);
                        }
        }               

        for (int i=0; i<toBeDeleted.size(); i++) {
                result->erase(toBeDeleted[i]);
                if (debug)
                        cout << "erased peak at:" << (*toBeDeleted[i]).maxPos << endl;
        }
*/