QPulseTimeConstants.cc

Go to the documentation of this file.

#include "QPulseTimeConstants.hh"
#include "QVector.hh"
#include "QInterval.hh"
 
using namespace std;
using namespace Diana;
 
QPulseTimeConstants::QPulseTimeConstants(const std::vector<QInterval>& riseIntervals,const std::vector<QInterval>& decayIntervals)
{
    if(riseIntervals.empty() && decayIntervals.empty()) {
        QError err = QERR_OUT_OF_RANGE;
        err.SetDescription(__FILE__,__LINE__,"Empty rise and decay intervals");
        DianaThrow(err);
    }
    for(size_t v = 0; v < riseIntervals.size(); v++) {
        if(riseIntervals[v].GetMin() >= riseIntervals[v].GetMax() || riseIntervals[v].GetMax() > 1) {
            QError err = QERR_OUT_OF_RANGE;
            err.SetDescription(__FILE__,__LINE__,"Rise intervals must be < 1, and Min < Max");
            DianaThrow(err);
        }
        fRise.push_back(make_pair(riseIntervals[v].GetMin(),v));
        fRise.push_back(make_pair(riseIntervals[v].GetMax(),v));
    }
 
    for(size_t v = 0; v < decayIntervals.size(); v++) {
        if(decayIntervals[v].GetMin() >= decayIntervals[v].GetMax() || decayIntervals[v].GetMax() > 1) {
            QError err = QERR_OUT_OF_RANGE;
            err.SetDescription(__FILE__,__LINE__,"Decay intervals must be < 1, and Min < Max");
            DianaThrow(err);
        }
        fDecay.push_back(make_pair(decayIntervals[v].GetMin(),v));
        fDecay.push_back(make_pair(decayIntervals[v].GetMax(),v));
    }
    sort(fRise.begin(),fRise.end(),PairComparatorAsc);
    sort(fDecay.begin(),fDecay.end(),PairComparatorAsc);
        
}
 
 
QPulseTimeConstants::~QPulseTimeConstants()
{
}
 
QError QPulseTimeConstants::Compute(const QVector& input, const int maxpos, const double baseLine,vector<double>& rise, vector<double>& decay)
{
    int auxMaxpos=maxpos;
    if(auxMaxpos < 1 || auxMaxpos >= (int)input.Size()) {
 
        //Compute max pos here
        size_t window = 10;
        size_t size = (int)input.Size();
        size_t from = 0;
        auxMaxpos=-1;
        int currMax = -1;
        double prevMaxValue = -1e15, currMaxValue = -1e15;
        while(from+window < size)
        {
          currMax = input.GetMaxIndex(window,from);
          currMaxValue = input[currMax];
          if(currMaxValue > prevMaxValue)
          {
              prevMaxValue = currMaxValue;
              auxMaxpos = currMax;
          }
          from += window;
        }
    }
    double base = baseLine;
    if(base < 0) {
        base = input[0];
        if(auxMaxpos > 0) base = input.GetMean(auxMaxpos*3./4);
    }
    double max = input[auxMaxpos]-base; 
 
    if(!fRise.empty()) {
        rise.resize(fRise.size()/2,Q_DOUBLE_DEFAULT);
        int riseCounter = fRise.size()-1;
        for(int i = auxMaxpos; i >= 0; i--) {
            if(input[i]-base < fRise[riseCounter].first*max) {
                double time = i + (fRise[riseCounter].first*max - input[i]+base)/(input[i+1] - input[i]);
                if(rise[fRise[riseCounter].second] == Q_DOUBLE_DEFAULT) rise[fRise[riseCounter].second] = time+input.Size(); 
                else rise[fRise[riseCounter].second] -= time+input.Size();
                riseCounter--;
                i++;
            }
            if(riseCounter < 0) break;
        }
    }
    for(size_t v = 0; v < rise.size();v++) {
        if(rise[v] >= input.Size()) rise[v] = Q_DOUBLE_DEFAULT;
    }
 
    if(!fDecay.empty()) {
        decay.resize(fDecay.size()/2,Q_DOUBLE_DEFAULT);
        int decayCounter = fDecay.size()-1;
        for(int i = auxMaxpos; i < (int)input.Size()-1; i++) {
            if(input[i]-base < fDecay[decayCounter].first*max) {
                double time = i - (fDecay[decayCounter].first*max - input[i]+base)/(input[i-1] - input[i]);
                if(decay[fDecay[decayCounter].second] == Q_DOUBLE_DEFAULT) decay[fDecay[decayCounter].second] = -time; 
                else decay[fDecay[decayCounter].second] += time;
                decayCounter--;
                i--;
            }
            if(decayCounter < 0) break;
        }
    }
 
    for(size_t v = 0; v < decay.size();v++) {
        if(decay[v] < 0) decay[v] = Q_DOUBLE_DEFAULT;
    }
 
 
 
    return QERR_SUCCESS;
}
 
bool QPulseTimeConstants::PairComparatorAsc ( pair<double, int> l,pair<double, int> r) 
{ 
    return l.first < r.first; 
}
 
bool QPulseTimeConstants::PairComparatorDesc ( pair<double, int> l,pair<double, int> r) 
{ 
    return l.first > r.first; 
}