Butterworth.cpp

/***************************************************************************** 
 * Project: RooFit                                                           * 
 *                                                                           * 
 * This code was autogenerated by RooClassFactory                            * 
 *****************************************************************************/ 

// Butterworth transfer function in T-domain 

#include "Riostream.h" 
#include "RooAbsReal.h" 
#include "RooAbsCategory.h" 

#include "TMath.h" 

#include "Butterworth.h" 


ClassImp(Butterworth) 

Butterworth::Butterworth( const char *name, const char *title, 
                          RooAbsReal& _t,
                          RooAbsReal& _omegaC,
                          UInt_t _n ) :
    RooAbsReal( name, title ), 

    t( "t", "t", this, _t ),
    omegaC( "omegaC", "omegaC", this, _omegaC ),
    n( _n )
{ 

} 


Butterworth::Butterworth( const Butterworth& other, const char* name ) :  
    RooAbsReal( other, name ), 

    t( "t", this, other. t ),
    omegaC( "omegaC", this, other.omegaC ),
    n( other.n )
{ 

} 


TComplex Butterworth::Sk( UInt_t k ) const
{
    // Returns k-th pole s_k of Butterworth transfer
    // function in S-domain. Note that omega_c
    // is not taken into account here
    Double_t arg = TMath::Pi() * ( 2 * k + n - 1 ) / ( 2 * n ); 

    return TComplex( TMath::Cos( arg ), TMath::Sin( arg ) );
}


TComplex Butterworth::J( UInt_t k ) const
{
    // Returns (s - s_k) * H(s), where
    // H(s) - BW transfer function
    // s_k  - k-th pole of H(s)
    TComplex res = TComplex( 1. );
    for( UInt_t m = 1; m <= n; m++ )
    {
        if( m == k ) continue;//skip s_m = s_k

        res /= (Sk( k ) - Sk( m ));
    }
    
    return res;
}


TComplex Butterworth::H( Double_t x ) const
{
    // Returns h(t) - BW transfer function in t-domain.
    // Note how omega_c is included here. 
    TComplex res = 0.; 

    if( x >= 0. )
    {
        //         n
        // h(t) = SUM[ (s - s_k) * H(s_k) * exp(s_k * t) ]
        //        k=1
        for( UInt_t k = 1; k <= n; k++ )
        {
            res += TComplex::Exp( omegaC * Sk( k ) * x ) * J( k );
        }
    }

    return (Double_t)omegaC * res;
}


Double_t Butterworth::evaluate() const 
{ 
    // It is easy to show that h(t) is a real function,
    // i.e. h(t) = Re( h(t) ).
    return H( t ).Re();
}