source: fact/tools/rootmacros/zerosearch.C@ 18350

#include "Region.h"
#include <iostream>

#include "zerosearch.h"

using namespace std;

// searches for zero crossings in a given vector of floats
// for zero crossings in falling edge
// give edge = -1
//
//
// returns pointer to vetctor of ints

vector<Region> *zerosearch(
    vector<float> &input,
    int edge,               // search on rising edge=1, -1:falling
    unsigned int step,      // search in steps of step
    int VerbosityLevel
){
vector<Region> * ZeroPositions = new vector<Region>;
Region currentRegion = {0, 0, 0, 0.0, 0, 0.0, 0.0, 0, 0};

    if (step < 1){
        if (VerbosityLevel > 0)
            cout << "zerosearch - stepsize=" << step
                << " is smaller than 1. returning." << endl;
        return ZeroPositions;
    }
    if (input.size() < step){
        if (VerbosityLevel > 0)
            cout << "zerosearch - input-vector.size()="<< input.size()
            << " is smaller than stepsize=" << step
            << "returning." << endl;
        return ZeroPositions;
    }

    if (edge > 0) // search for zero-x-ings on the rising edge
    {
        for (unsigned int sl = 0 ; sl < (input.size()-step) ; sl+=step)
        {
            if (input[sl] > 0.0) // can't be a rising 0-x-ing
            {
                continue;
            }
            if (input[sl] * input[sl+step] <= 0.0)  // 0-x-ing or both are 0
            {
                currentRegion.begin     = sl;
                currentRegion.end       = sl+step;
                ZeroPositions->push_back(currentRegion);
            }
        }
    }
    else if (edge < 0) // search for zero-x-ings on the falling edge
    {
        for (unsigned int sl = 0 ; sl < (input.size()-step) ; sl+=step)
        {
            if (input[sl] < 0.0) // can't be a falling 0-x-ing
            {
                continue;
            }
            if (input[sl] * input[sl+step] <= 0.0)  // 0-x-ing or both are 0
            {
                currentRegion.begin     = sl;
                currentRegion.end       = sl+step;
                ZeroPositions->push_back(currentRegion);
            }
        }
    }
    else // search for zero-x-ings on both esges
    {
        for (unsigned int sl = 0 ; sl < (input.size()-step) ; sl+=step)
        {
            if (input[sl] * input[sl+step] <= 0.0)  // 0-x-ing or both are 0
            {
                currentRegion.begin     = sl;
                currentRegion.end       = sl+step;
                ZeroPositions->push_back(currentRegion);
            }
        }
    }

    return ZeroPositions;
}

size_t ShiftRegionBy(vector<Region> &src,
    int Shift,
    int VerbosityLevel)
{
    vector<Region>::iterator it;
    // I copy code here ... not good, but I hope nobody is giving VL > 10 ever...
    if (VerbosityLevel > 10){
        for (it = src.begin() ; it < src.end() ; it++){
            cout << " it->begin="<< it->begin;
            it->begin += Shift;
            cout << "--> shifted="<< it->begin << endl;
            cout << " it->end="<< it->end;
            it->end += Shift;
            cout << "--> shifted="<< it->end << endl;
            cout << " it->maxPos="<< it->maxPos;
            it->maxPos += Shift;
            cout << "--> shifted="<< it->maxPos << endl;
        }
    return src.size();
    }

    for (it = src.begin() ; it < src.end() ; it++){
        it->begin += Shift;
        it->end += Shift;
        it->maxPos += Shift;
    }
    return src.size();
}

size_t EnlargeRegion(vector<Region> &src,
    int Left,
    int Right,
    int VerbosityLevel)
{
    vector<Region>::iterator it;
    // I copy code here ... not good, but I hope nobody is giving VL > 10 ever...
    if (VerbosityLevel > 10){
        for (it = src.begin() ; it < src.end() ; it++){
            cout << " it->begin="<< it->begin;
            it->begin -= Left;
            cout << "--> enlarged="<< it->begin << endl;
            cout << " it->end="<< it->end;
            it->end += Right;
            cout << "--> enlarged="<< it->end << endl;
        }
        return src.size();
    }


    for (it = src.begin() ; it < src.end() ; it++){
        it->begin -= Left;
        it->end += Right;
    }
    return src.size();

}

#include <limits>
size_t findAbsMaxInRegions(
    vector<Region> &regions,
    vector<float> &data,
    int VerbosityLevel)
{
    vector<Region>::iterator reg;
    for (reg = regions.begin() ; reg < regions.end() ; reg++){
        reg->maxVal=-numeric_limits<float>::max();
        // 1st check if both ends of the region are in the data-vector
        if (reg->begin < 0) reg->begin=0;
        if ((unsigned int)reg->end > data.size()-1) reg->end=data.size()-1;

        // TODO or TOTHINKOF:
        // Region might overlap ... I don't care .. I just search
        // in a later step Maxima at the same Positions can be removed...
        // of course this means, some searches were needless ...
        //cout << "searching from " << reg->begin << " to " << reg->end << endl;
        for (unsigned int pos=reg->begin; pos <= (unsigned int)reg->end; pos++){
            if (data[pos] > reg->maxVal)
            {
                reg->maxVal = data[pos];
                reg->maxPos = pos;
            }
        }
        if (VerbosityLevel > 3) {
            cout << "findAbsMaxInRegions - found Max at:"
                << reg->maxPos << " with Value:" << reg->maxVal << endl;
        }

        if (reg->maxPos <= 0)
        {
            reg = regions.erase( reg ) ;
            --reg;
            continue;
        }
    }
    return regions.size();
}

#include <algorithm>
bool compMaxPosOfRegions (Region a, Region b) {
  return (a.maxPos == b.maxPos);
}

size_t removeEqualMaxima(
    vector<Region> &regions,
    int VerbosityLevel)
{
    if (VerbosityLevel > 2){
        cout << "removeEqualMaxima:" << endl;
        cout << "region before contains:" << endl;
        for (unsigned int i=0; i<regions.size(); i++){
            cout << i <<"\t";
            cout << regions[i].begin <<"\t";
            cout << regions[i].end <<"\t";
            cout << regions[i].maxPos <<"\t";
            cout << regions[i].maxVal <<"\t";
            cout << endl;
        }
    }

    vector<Region>::iterator it;
    it = unique (regions.begin(), regions.end() , compMaxPosOfRegions);
    regions.resize( it - regions.begin() );

    if (VerbosityLevel > 1){
        cout << "region now contains:" << endl;
        for (unsigned int i=0; i<regions.size(); i++){
            cout << i <<"\t";
            cout << regions[i].begin <<"\t";
            cout << regions[i].end <<"\t";
            cout << regions[i].maxPos <<"\t";
            cout << regions[i].maxVal <<"\t";
            cout << endl;
        }
    }
    return regions.size();
}

size_t removeRegionOnFallingEdge(
    vector<Region> &regions,
    unsigned int FallingEdgeWidth,
    int VerbosityLevel)
{
    //throw exceptions
    if ( FallingEdgeWidth < 1){
        if (VerbosityLevel > 0){
            cout << " removeRegionOnFallingEdge: FallingEdgeWidth  < 1" << endl;
            cout << " FallingEdgeWidth =" << FallingEdgeWidth << endl;
            cout << "returning." << endl;
            return regions.size();
        }
    }
    //throw exceptions
    if (regions.size() < 1){
        if (VerbosityLevel > 3)
            cout << "removeRegionOnFallingEdge: source vector empty." << endl;
        return 0;
    }

 vector<Region>::reverse_iterator it = regions.rbegin();
  while( it != regions.rend()-1 )
  {
  if (VerbosityLevel >10)
        cout << it->maxPos << "\t" << (it+1)->maxPos << "\t" << it->maxPos - (it+1)->maxPos << endl;

    if ( it->maxPos - (it+1)->maxPos < (int)FallingEdgeWidth) {
      if (VerbosityLevel > 3)
        cout << "erasing Region @ " << it->maxPos << endl;
       regions.erase( --it.base() ) ;
      ++it;
    }else
      ++it;
  }

  return regions.size();
}


size_t removeRegionWithMaxOnEdge(
    vector<Region> &regions,
    unsigned int EdgeWidth,
    int VerbosityLevel)
{
    if (EdgeWidth < 1){
        if (VerbosityLevel > 0){
            cout << "removeRegionWithMaximaOnEdge: EdgeWidth < 1" << endl;
            cout << "EdgeWidth=" << EdgeWidth << endl;
            cout << "returning." << endl;
        return regions.size();
        }
    }

    vector<Region>::iterator it = regions.begin();
    while( it != regions.end() )
    {
//      cout << it->begin << "\t";
//      cout << it->maxPos << "\t";
//      cout << it->end << "\t";
//      cout << it->maxVal << endl;

        if (it->maxPos < (int)(it->begin+EdgeWidth)) {
            if (VerbosityLevel > 3)
                cout << "erasing Region(max@left edge) " << it->maxPos << endl;
            it = regions.erase( it ) ;
            //++it;
        }else if (it->maxPos > (int)(it->end-EdgeWidth)) {
            if (VerbosityLevel > 3)
                cout << "erasing Region(max@right edge) " << it->maxPos << endl;
            it = regions.erase( it ) ;
            //++it;
        }else
            ++it;

    }

    return regions.size();
}

size_t removeMaximaBelow(
    vector<Region> &regions,
    float threshold,
    int VerbosityLevel)
{
//      if (threshold < 0){
//              if (VerbosityLevel > 0)
//                      cout << "removeMaximaBelow: threshold < 0" << endl;
//                      cout << "threshold=" << threshold << endl;
//                      cout << "returning." << endl;
//              return regions.size();
//      }

    vector<Region>::iterator it = regions.begin();
    while( it != regions.end() )
    {
        if (it->maxVal < threshold ) {
            if (VerbosityLevel > 3){
                cout << "removing max " << it->maxVal << "\t";
                cout << " @ " << it->maxPos << "\t";
                cout << endl;
            }
            it = regions.erase( it ) ;
        }else
            ++it;
    }

    return regions.size();
}

size_t removeMaximaAbove(
    vector<Region> &regions,
    float threshold,
    int VerbosityLevel)
{
    if (threshold < 0){
        if (VerbosityLevel > 0)
            cout << "removeMaximaAbove: threshold < 0" << endl;
            cout << "threshold=" << threshold << endl;
            cout << "returning." << endl;
        return regions.size();
    }

    vector<Region>::iterator it = regions.begin();
    while( it != regions.end() )
    {
        if (it->maxVal > threshold ) {
            if (VerbosityLevel > 3){
                cout << "removing max " << it->maxVal << "\t";
                cout << " @ " << it->maxPos << "\t";
                cout << endl;
            }
            it = regions.erase( it ) ;
        }else
            ++it;
    }

    return regions.size();
}

Region FindAbsMax(
    vector<Region> &regions,
    int VerbosityLevel)
{
    Region AbsMax;
    AbsMax.begin = -1;
    AbsMax.maxPos = -1;
    AbsMax.end = -1;
    AbsMax.maxVal=-numeric_limits<float>::max();
    if (regions.size() < 1)
        return AbsMax;

    for (vector<Region>::iterator it = regions.begin(); it < regions.end(); ++it)
    {
        if (it->maxVal > AbsMax.maxVal ) {
            AbsMax = *it;
        }
    }
    if (VerbosityLevel > 5){
        AbsMax.begin +=0;
    }
    return AbsMax;
}

////////////////
////////////// old Version of the code
/*
    for (unsigned int sl = step; sl < input.size(); sl+=step){
        if (input[sl] * input[sl-] < 0 || input[sl]==0.0 ){ // sign change --> zero crossing OR really zero

            if ( (input[sl-1] - input[sl]) * edge < 0){ // this is the crossing the user wanted

                // check if we go lower than a certain limit in the next few slices
                for ( int lala=0; lala<post; lala++){
                    if (input[sl+lala] > 1.5) {
                        zeroPositions->push_back(sl);
                        sl += lala;
                        break;
                    }
                }

            } else if ( (input[sl-1] - input[sl]) * edge > 0){ // this is the zero x-ing the user did not want

                // do nothing

            } else { // sl and sl-1 are equal .. don't know waht do to...

            }

        }

    } // end of loop over slices - between pre and post

*/

//size_t calcAmplitudeWeightedTime(
//        vector<Region>  &regions,
//        vector<float>   &data,

//        ){

//}

size_t calcCFDPositions(
        vector<Region>  &regions,
        vector<float>   &cfd_data
        )
{
    vector<Region>::iterator reg;
    for (reg = regions.begin() ; reg < regions.end() ; reg++)
    {
        int     neg_zero_crossing = -1;
        int     min_in_cfd        = -1;
        int     pos_zero_crossing = reg->maxPos;
        float   minVal            = 10.0;

        for (int i = pos_zero_crossing; i > pos_zero_crossing - 40; i--)
        {
            if (i <= 1 || i >= 1439) break;

            if ( cfd_data[i] < minVal)
            {
                minVal      = cfd_data[i];
                min_in_cfd  = i;
            }

            if (cfd_data[i] > 0 && cfd_data[i+1] <= 0 && cfd_data[i+2] < 0 )
            {
                neg_zero_crossing = i+1;
            }
        }

        reg->cfdMinPos          = min_in_cfd;
        reg->cfdNegZerocrossing = neg_zero_crossing;
    }
    return regions.size();
}



size_t findTimeOfHalfMaxLeft(
    vector<Region>  &regions,
    vector<float>   &data,
    float           baseline,
    int             beginRisingEdge,
    int             endRisingEdge,
    int             VerbosityLevel
        )
{
//    float   pulse_size      = 0;     //
    float   thr             = 0;
    //int begin = beginRisingEdge;
    int     end             = 20;
    int     counter         = 0;
    // local vars for line fitting
    float   xmean           = 0.0;
    float   ymean           = 0.0;
    float   cov             = 0.0;          // empiric covariance between x and y
    float   var             = 0.0;          // empiric variance of x
    // result of line fitting
    float   slope           = 0.0;
    float   intercept       = 0.0;
    // result of this function -- the position of the half rising edge
    float   pos_of_thr_Xing = 0.0;
    int     result          = 0;


    vector<Region>::iterator reg;
    for (reg = regions.begin() ; reg < regions.end() ; reg++)
    {
        // check if linear rising edge is completely in pipeline
        // if not --> delete
        if ( reg->maxPos - end < 0 )
        {
            // delete this region from vector<Region>
            reg = regions.erase( reg ) ;
            --reg;
            continue;
        }

        // I think 5 slices left of the max there begins the rising edge
        // and I think the rising edge is about 6 slices long.
        // but these numbers are input as parameters
        // beginRisingEdge &
        // endRisingEdge

        // I fit a linear function to the falling edge
        // ALARM check for out of range values...

        float average_of_max = -1;
        int stepwidth = 2;

        counter = 0;
        for ( int i = reg->maxPos -stepwidth; i <= reg->maxPos + stepwidth; i++)
        {
            average_of_max += data[i];
            counter++;
        }
        average_of_max /= counter;

        thr = average_of_max / 2.0;

        float   max_slope       = 0;
        int     max_slope_pos   = -1;
        int     half_height_pos = -1;

        bool passed = false;
        for (int slice=reg->maxPos; slice > reg->maxPos - end; --slice)
        {
            if (slice <= 0) break;
//            if ( data[slice] < 0.9*data[reg->maxPos] && !passed)
//            {
//                endRisingEdge = slice;
//                beginRisingEdge = endRisingEdge -3;
//                passed = true;
//            }
//            else if ( data[slice] < 0) //0.4*data[reg->maxPos] )
//            {
//                beginRisingEdge = slice;
//                reg->lengthOfRisingEdge=endRisingEdge - beginRisingEdge;
//                break;
//            }

            //calculate max slope on leading edge
            if (data[slice] - data[slice-4] > max_slope && slice-4 > 0)
            {
                max_slope   = data[slice] - data[slice-4];
                max_slope_pos      = slice-2;
            }

            if ( data[slice] <= thr && !passed){
                half_height_pos  = slice;
                passed  = true;
            }
        }
        reg->maxSlope           = max_slope/4;
        reg->maxSlopePos        = max_slope_pos;
        reg->halfRisingEdgePos  = half_height_pos;
        reg->halfHeight         = data[half_height_pos];

        result                  = half_height_pos;

        endRisingEdge   = result +2;
        beginRisingEdge = result -2;

        counter = 0;
        //calculate mean of x and y coordinate
        for (int slice= beginRisingEdge; slice < endRisingEdge; ++slice)
        {
            if (slice <= 0) break;

            xmean += slice;
            ymean += data[slice];
            counter++;
        }
        xmean /= counter;
        ymean /= counter;

        if (VerbosityLevel > 2) cout << "## xmean: " << xmean << endl;
        if (VerbosityLevel > 2) cout << "## ymean: " << ymean << endl;

        cov = 0.0;
        var = 0.0;

        //calculate covariance and variance
        for (int slice=beginRisingEdge; slice < endRisingEdge; ++slice)
        {
            if (slice <= 0) break;

            cov += (data[slice] - ymean) * (slice - xmean);
            var += (slice - xmean) * (slice - xmean);
        }
        if (VerbosityLevel > 2) cout << "## cov: " << cov << endl;
        if (VerbosityLevel > 2) cout << "## var: " << var << endl;

        slope       = cov / var;
        intercept   = ymean - slope * xmean;

        // The maximum was probably found using smoothed data,
        // so I read the value at the position of the maximum from the data
        // again. So I rely on a well defined maxPos.
//        if (VerbosityLevel > 4) cout << "## baseline: " << baseline << endl;
//        pulse_size = data[reg->maxPos] - baseline;
//        if (VerbosityLevel > 4) cout << "## pulse_size: " << pulse_size << endl;
//        thr = pulse_size / 2. + baseline;
//        if (VerbosityLevel > 4) cout << "## thr: " << thr << endl;

        // now calculate, where the fittet line crosses the threshold
//        pos_of_thr_Xing = (thr - intercept) / slope;
//        result = (int)(pos_of_thr_Xing + 0.5);


        if (VerbosityLevel > 2)
        {
            cout << "findTimeOfHalfMaxLeft() is still in debugging phase:" << endl;
            cout << "please edit the code in oder to suppress this output:" << endl;
            cout << "threshold: " << thr << endl;
            cout << "slope: " << slope << " [mV/timeslice]" << endl;
            cout << "intercept: " << intercept << "[mV]" << endl;
            cout << "position where line crosses threshold " << pos_of_thr_Xing << endl;
            cout << "converted to slice number " << result << endl;

        }

        reg->distanceEdgeToMax      = reg->maxPos - result;
        reg->interceptRisingEdge    = reg->maxPos - intercept;
        reg->slopeOfRisingEdge      = slope;

        if (reg->halfRisingEdgePos <= 0 || reg->halfRisingEdgePos >= 1024)
        {
            reg = regions.erase( reg ) ;
            --reg;
            continue;
        }
    }
    return regions.size();
}

size_t removeRegionWithToFlatSlope(
    vector<Region> &regions,
    float minSlope,
    int VerbosityLevel)
{
    vector<Region>::iterator it = regions.begin();
    while( it != regions.end() )
    {
        if (it->slopeOfRisingEdge < minSlope ) {
            if (VerbosityLevel > 3){
                cout << "removing max " << it->maxVal << "\t";
                cout << " @ " << it->maxPos << "\t";
                cout << endl;
            }
            it = regions.erase( it ) ;
        }else
            ++it;
    }

    return regions.size();
}

size_t GetMaxPositions(
    vector<Region>  &regions,
    vector<int>     &positions,
    int VerbosityLevel)
{
    positions.clear();
    vector<Region>::iterator it = regions.begin();
    while( it != regions.end() )
    {
        positions.push_back(it->maxPos);
        if (VerbosityLevel > 3){
            cout << "getting maxPos@ " << it->maxPos << "\t";
            cout << endl;
        }
        ++it;
    }

    return regions.size();
}

size_t GetEdgePositions(
    vector<Region>  &regions,
    vector<int>     &positions,
    int VerbosityLevel)
{
    positions.clear();
    vector<Region>::iterator it = regions.begin();
    while( it != regions.end() )
    {
        positions.push_back(it->halfRisingEdgePos);
        if (VerbosityLevel > 3){
            cout << "getting maxPos@ " << it->halfRisingEdgePos << "\t";
            cout << endl;
        }
        ++it;
    }

    return regions.size();
}