CFD.hpp

//  CFD.hpp
//  
//  Author : Jeromy Tompkins
//  Date   : 8/14/2013 

#ifndef CFD_H
#define CFD_H

#include <deque>
#include <cmath>
#include "TrapFilter.hpp"
#include "Solver.h"
#include "TObject.h"

namespace TrAnal
{


template<class T>
struct CFDResult
{
    TrIterator<T> trig; 
    double fraction; 

    // ROOT dictionary generation 
    ClassDef(CFDResult,0);
};



template<class T> 
CFDResult<T> CFD(const TrRange<T>& range,  
        int rise_len, 
        int gap_len,   
        int delay,     
        int scale_factor,
        const Solver& solver=LinearSolver()  )

{ 
    // TrapFilters are defined from their lower bound. First locate the 
    // the location of the lead filter
    TrIterator<T> lead_start = range.begin()+delay;
    // Verify that it is in range
    if (lead_start>=range.end() ) return CFDResult<T>();

    // Construct the two trap filters
    TrapFilter<T> lead_filter (lead_start, rise_len, gap_len);
    TrapFilter<T> dely_filter (range.begin(), rise_len, gap_len);


    // initialize the queue
    std::deque<double> cfd_res;
    for (unsigned int i=0; i<solver.npoints() && lead_filter<range.end(); ++i) {
        
        // Compute the CFD value
        cfd_res.push_back((*lead_filter) - (*dely_filter)/pow(2.0,scale_factor+1));

        ++lead_filter;
        ++dely_filter;
    }

    // Initialize the counter for number of negative points
    unsigned int nneg=0;
    while (lead_filter<range.end()) {
        // Compute new CFD values
        cfd_res.push_back((*lead_filter) - (*dely_filter)/pow(2.0,scale_factor+1));
    
        // Have we crossed over zero? 
        if (cfd_res.back() < 0) ++nneg;
    
        // Do we have enough points to determine the zero crossing?
        // if so, stop!
        if (nneg>=solver.npoints()-1) break;

        ++lead_filter;
        ++dely_filter;

        cfd_res.pop_front();
    }

    // If not enough data exists to find a zero crossing, return a null result
    if (nneg!=solver.npoints()-1) return CFDResult<T>();

    // This will return the fraction index of the zero crossing with respect to the
    // front of the queue
    double trig = solver(cfd_res);
    int floortrig = static_cast<int>(::floor(trig));

    // trig location happened at ceiltrig from end
    int n_from_end=static_cast<int>(cfd_res.size())-floortrig;
    TrIterator<T> trigit = lead_filter.max_extent()-n_from_end;
    double frac = trig - floortrig; 

    // Construct the result and move return it
    CFDResult<T> res;
    res.trig = trigit;
    res.fraction = frac;

    return res;

} // end CFD function

} // end namespace 

#endif