QBiComponentOptimumFilter.cc

Go to the documentation of this file.

#include "QBiComponentOptimumFilter.hh"
#include "QError.hh"
#include "QMatrix.hh"
#include "QRealComplexFFT.hh"
#define EXTRA_JITTER 1
 
// EXTRA_JITTERS:
// 1: minimum number needed for interpolation
// 2: more point to see if we have a minimum at the boundary (enable nice parabola)
 
using std::vector;
using namespace Diana;
 
QBiComponentOptimumFilter::QBiComponentOptimumFilter(const QVector& ap1, const QVector& ap2, const QVector& an, size_t maxJitter, bool diff) : 
fIsDiff(diff),
fMaxJitter(maxJitter),
fTransformer(0)
{
    fCheck = QERR_UNKNOWN_ERR;
    fSize = ap1.Size();
    if(ap2.Size() != fSize || an.Size() != fSize) {
        QError err(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,"Size of input QVectors mismatches");
        DianaThrow(err);
    }
    fTransformer = new QRealComplexFFT(fSize);
 
    QVector ap1n = NormalizeVector(ap1);
    QVector ap2n = NormalizeVector(ap2);
    fAveragePulse1 = ap1n;
    fAveragePulse2 = ap2n;
 
    fAvgPulseIntegral1= ap1n.Sum(ap1n.Size(),0);
    fAvgPulseIntegral2= ap2n.Sum(ap2n.Size(),0);
 
    if(fIsDiff) {
        ap1n.Differentiate();
        ap2n.Differentiate();
 
        QVector andi = an;
        for(size_t i = 1 ; i < fSize;i++) {
            if(i<=fSize/2) andi[i]*=2.*(1.-cos(2.*M_PI/(fSize)*i));
            else andi[i]*=2.*(1.-cos(2.*M_PI/(fSize)*(fSize-i)));
        }
        fAverageNoise = andi;
    } else {
        fAverageNoise = an;
    }
    fAverageNoise[0] = 0;
    BuildFilter(ap1n,ap2n);
}
 
QBiComponentOptimumFilter::~QBiComponentOptimumFilter()
{
    if(fTransformer) {
        delete fTransformer;
        fTransformer = 0;
    }
}
 
QError QBiComponentOptimumFilter::Filter(const Diana::QVector& p)
{
    fCheck = QERR_SUCCESS;
    if(p.Size() != fSize) {
        QError err(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,"Size of filtered QVector is not correct");
        fCheck = err;
        return err;
    }
    fFiltered1.Clear(); 
    fFiltered2.Clear(); 
    QVectorC pf,pf1,pf2;
    if(fIsDiff) {
        QVector pd = p;
        pd.Differentiate();
        fTransformer->TransformToFreq(pd,pf);
    } else {
        fTransformer->TransformToFreq(p,pf);
    }
    pf1 = fFilter1.Mult(pf);
    pf2 = fFilter2.Mult(pf);
    fTransformer->TransformFromFreq(pf1,fFiltered1);
    fTransformer->TransformFromFreq(pf2,fFiltered2);
 
    fCheck = ComputeAmplitude(pf);
    
    return fCheck;
}
void QBiComponentOptimumFilter::GetFiltered(Diana::QVector& f1, Diana::QVector& f2) const 
{ 
    if(fCheck == QERR_SIZE_NOT_MATCH) DianaThrow(fCheck);
    f1 = fFiltered1; 
    f2 = fFiltered2; 
}
const std::vector<double>& QBiComponentOptimumFilter::GetChiSquareWithJitter() const 
{ 
    if(fCheck == QERR_SIZE_NOT_MATCH) DianaThrow(fCheck);
    return fChiSquareWithJitter;
}
 
int QBiComponentOptimumFilter::GetJitterAtMinimum() const 
{ 
    if(fCheck != QERR_SUCCESS) DianaThrow(fCheck);
    return fJitters[fOptimumJitter] ;
}
 
double QBiComponentOptimumFilter::GetChiSquareAtMinimium() const 
{ 
    if(fCheck != QERR_SUCCESS) DianaThrow(fCheck);
    return fChiSquareWithJitter[fOptimumJitter];
}
 
void QBiComponentOptimumFilter::GetAmplitude(double& a1, double& a2, double& integral) const 
{ 
    if(fCheck != QERR_SUCCESS) DianaThrow(fCheck);
    a1 = fAmplitude1; 
    a2 = fAmplitude2; 
    integral = fAvgPulseIntegral1*fAmplitude1+fAvgPulseIntegral2*fAmplitude2;
}
 
QError QBiComponentOptimumFilter::GetInterpolated(double& intJitter, double& inta1, double& inta2,double& integral) const
{
    if(fCheck != QERR_SUCCESS) DianaThrow(fCheck);
    double xa = fJitters[fOptimumJitter-1];
    double xb = fJitters[fOptimumJitter];
    double xc = fJitters[fOptimumJitter+1]; 
    double ya = fChiSquareWithJitter[fOptimumJitter-1];
    double yb = fChiSquareWithJitter[fOptimumJitter];
    double yc = fChiSquareWithJitter[fOptimumJitter+1];
    if(yb>ya || yb>yc) {
        QError err = QERR_OUT_OF_RANGE;
        err.SetDescription("Chisquare not at minimum, interpolation not possible");     
        return err;
    }
 
        
 
    double detA = xa*xa*(xb-xc)-xb*xb*(xa-xc)+xc*xc*(xa-xb);
    double detA1 = ya*(xb-xc)-yb*(xa-xc)+yc*(xa-xb);
    double detA2 = xa*xa*(yb-yc)-xb*xb*(ya-yc)+xc*xc*(ya-yb);
    //double detA3 = xa*xa*(xb*yc-yb*xc)-xb*xb*(xa*yc-ya*xc)+xc*xc*(xa*yb-ya*xb);
 
    double p0 = detA1/detA;
    double p1 = detA2/detA;
//    double p2 = detA3/detA;
 
    intJitter = -p1/2./p0;
    int low = floor(intJitter); 
    int high = ceil (intJitter); 
    double x =  intJitter-low;
    if(low < 0) low += fSize;
    if(high < 0) high += fSize;
    double a1Low = fFiltered1[low];
    double a1High = fFiltered1[high];
    double a2Low = fFiltered2[low];
    double a2High = fFiltered2[high];
    inta1 = (a1High-a1Low)*x+a1Low;
    inta2 = (a2High-a2Low)*x+a2Low;
    integral = fAvgPulseIntegral1*inta1+fAvgPulseIntegral2*inta2;
    return QERR_SUCCESS;
}
 
QVector QBiComponentOptimumFilter::NormalizeVector(const QVector& vec) const
{
    QVector ret = vec;
    QVector baseline(ret.Size());
    baseline.Initialize(ret[0]);
    ret -= baseline;
    ret /= ret.GetMax();
    return ret;
}
 
void QBiComponentOptimumFilter::BuildFilter(const QVector& ap1, const QVector& ap2)
{
    const size_t size = ap1.Size();
 
    QVectorC ap1f,ap2f;
    fTransformer->TransformToFreq(ap1,ap1f);
    fTransformer->TransformToFreq(ap2,ap2f);
    
 
    Diana::QMatrix M(2,2);
    M(0,0) = (ap1f.GetModulusSquare().Div(fAverageNoise)).Sum(size-1,1);
    M(0,1) = ((ap1f.Conj().Mult(ap2f)).Re().Div(fAverageNoise)).Sum(size-1,1);
    M(1,0) = M(0,1);
    M(1,1) = (ap2f.GetModulusSquare().Div(fAverageNoise)).Sum(size-1,1);    
    M.Invert();
 
    QVectorC anc(size);
    anc.Initialize(0,0);
    anc.SetRe(fAverageNoise);
 
    fFilter1 = (M(0,0)*ap1f.Conj()+M(0,1)*ap2f.Conj()).Div(anc);
    fFilter2 = (M(1,0)*ap1f.Conj()+M(1,1)*ap2f.Conj()).Div(anc);
 
    fFilteredRMS1 = ((fFilter1.Mult(anc)).GetModulusSquare().Div(fAverageNoise)).Sum(size-1,1);
    fFilteredRMS1 = sqrt(fFilteredRMS1);
    fFilteredRMS2 = ((fFilter2.Mult(anc)).GetModulusSquare().Div(fAverageNoise)).Sum(size-1,1);
    fFilteredRMS2 = sqrt(fFilteredRMS2);
 
    fFilter1 *= size;
    fFilter2 *= size;
    QComplex null(0,0);
    fFilter1[0] = null;
    fFilter2[0] = null;
 
    int maxJitter = fMaxJitter;
    int minJitter = -maxJitter;
    for(int j = minJitter-EXTRA_JITTER; j <= maxJitter+EXTRA_JITTER; j++) {
        QVector ap1Shifted = ap1;
        ap1Shifted.Shift(j);
        QVector ap2Shifted = ap2;
        ap2Shifted.Shift(j);
        QVectorC ap1ShiftedF,ap2ShiftedF;
        fTransformer->TransformToFreq(ap1Shifted,ap1ShiftedF);
        fTransformer->TransformToFreq(ap2Shifted,ap2ShiftedF);
        fAveragePulse1Shifted.push_back(ap1ShiftedF);
        fAveragePulse2Shifted.push_back(ap2ShiftedF);
        fJitters.push_back(j);
    }
}
 
QError QBiComponentOptimumFilter::ComputeAmplitude(const QVectorC& transformed)
{
    QError err = QERR_SUCCESS;
    fChiSquareWithJitter = ChiSquare(transformed,fFiltered1,fFiltered2);
 
    fOptimumJitter = EXTRA_JITTER;
    double chi  = fChiSquareWithJitter[fOptimumJitter];
    for(int i = EXTRA_JITTER+1; i < int(fChiSquareWithJitter.size())-EXTRA_JITTER; i++)  {
        if(fChiSquareWithJitter[i] < chi) {
            fOptimumJitter = i;
            chi = fChiSquareWithJitter[i];
        }
    }
    int index = fJitters[fOptimumJitter];
    if(index < 0) index += fSize;
    fAmplitude1 = fFiltered1[index];
    fAmplitude2 = fFiltered2[index];
    return err;
}
 
 
 
vector<double> QBiComponentOptimumFilter::ChiSquare(const QVectorC& p, const QVector& ap1Filtered, const QVector& ap2Filtered) const
{
    vector<double> chiSquare;
    for(size_t j = 0; j < fJitters.size(); j++) {
        double a1 = 0,a2= 0;
        int jitter = fJitters[j];
        if(jitter < 0) jitter += fSize;
        a1 = ap1Filtered[jitter];
        a2 = ap2Filtered[jitter];
        const Diana::QVectorC& ap1f = fAveragePulse1Shifted[j];
        const Diana::QVectorC& ap2f = fAveragePulse2Shifted[j];
        double chi2 = (p-a1*ap1f-a2*ap2f).GetModulusSquare().Div(fAverageNoise).Sum(fSize-1,1);
        chi2 /= (fSize-2);
        chiSquare.push_back(chi2);
    }
    return chiSquare;
}
 
 
QVector QBiComponentOptimumFilter::GetFitFunction(const double jitter, const double a1, const double a2) const
{
        QVector output = fAveragePulse1*a1+fAveragePulse2*a2;
        output.ShiftReal(jitter);
        return output;
}