Lutz.cxx

Go to the documentation of this file.

///
/// \file Lutz.cxx 
/// \brief "Breakpoint" track fitter, G. Lutz, NIM A273 (1988) 349-361
/// \author messier@indiana.edu
///

#include "Lutz.h"
#include "BreakPointFitter/func/BPException.h"

#include <iostream>

namespace bpfit {
  Lutz::Lutz(unsigned int n, unsigned int N) :
    fn(n),
    fN(N),
    fz(n),
    fXI(n),
    fwx(n),
    fZ(N),
    fwSJ(N),
    fZeta(n),
    fD(N+2,n),
    fAinv(N+2),
    fa(0),
    fb(0),
    fAlpha(fN),
    fx(n)
  {
    if(N >= n)
      throw BPException(__FILE__, __LINE__, kNn);

    for(unsigned int i = 0; i < n; i++)
      fZeta[i].resize(N);
  }
  
  //......................................................................
  
  void Lutz::SetMeasurement(unsigned int i, double z, double x, double sigx) 
  {
    if(i >= fn)
      throw BPException(__FILE__, __LINE__, kINDEX);
    if(z != z || x != x || sigx != sigx)
      throw BPException(__FILE__, __LINE__, kNAN);
    if(sigx <= 0.0)
      throw BPException(__FILE__, __LINE__, kSIGX, sigx);

    fz[i]  = z;
    fXI[i] = x;
    fwx[i] = 1.0/(sigx*sigx);
  }
  
  //......................................................................
  
  void Lutz::SetScatteringPlane(unsigned int j, double Z, double sigSJ) 
  {
    if(j >= fN)
      throw BPException(__FILE__, __LINE__, kINDEX);
    if(Z != Z || sigSJ != sigSJ)
      throw BPException(__FILE__, __LINE__, kNAN);
    if(sigSJ <= 0.0)
      throw BPException(__FILE__, __LINE__, kSIGSJ, sigSJ);
    
    fZ[j]   = Z;
    fwSJ[j] = 1.0/(sigSJ*sigSJ);
  }
  
  //......................................................................
  
  void Lutz::CalcZeta() 
  {
    unsigned int i, j;
    for (i=0; i<fn; ++i) {
      for (j=0; j<fN; ++j) {
        fZeta[i][j] = fz[i]-fZ[j];
        if (fZeta[i][j]<0.0) fZeta[i][j] = 0.0;
      }
    }
  }
  
  //......................................................................
  //
  // Equations 31.
  //
  void Lutz::CalcD() 
  {
    unsigned int i, j;
    for (i=0; i<fn; ++i) {
      fD[0][i] = fwx[i];
      fD[1][i] = fwx[i]*fz[i];
      for (j=0; j<fN; ++j) {
        fD[j+2][i] = fwx[i]*fZeta[i][j];
      }
    }
  }
  
  //......................................................................
  //
  // Equations 29.
  //
  bool Lutz::CalcAinv()
  {
    unsigned int i, J, K;
    
    double sum00 = 0;
    double sum01 = 0;
    double sum11 = 0;
    double fwxfz;
    for (i=0; i<fn; ++i) {
      fwxfz  = fwx[i]*fz[i]; 
      sum00 += fwx[i];
      sum01 += fwxfz;
      sum11 += fwxfz*fz[i];
    }
    fAinv[0][0] = sum00;
    fAinv[0][1] = sum01;
    fAinv[1][0] = sum01;
    fAinv[1][1] = sum11;
    
    double sum0J2;
    double sum1J2;
    double fwxZeta;
    for (J=0; J<fN; ++J) {
      sum0J2 = 0;
      sum1J2 = 0;
      for (i=0; i<fn; ++i) {
        fwxZeta = fwx[i]*fZeta[i][J];
        sum0J2 += fwxZeta;
        sum1J2 += fz[i]*fwxZeta;
      }
      fAinv[0][J+2] = sum0J2;
      fAinv[J+2][0] = sum0J2;
      fAinv[1][J+2] = sum1J2;
      fAinv[J+2][1] = sum1J2;
    }
    
    double sum;
    for (J=0; J<fN; ++J) {
      for (K=J; K<fN; ++K) {
        sum = 0;
        for (i=0; i<fn; ++i) sum += fwx[i]*fZeta[i][J]*fZeta[i][K];
        if (J==K) sum += fwSJ[J];
        fAinv[J+2][K+2] = sum;
        fAinv[K+2][J+2] = sum;
      }
    }
    
    double det = 0.0;

    /// Inverting the matrix
    /// Sometimes this operation fails, producing error: matrix is singular
    /// Thus we need to put it into try/catch block to avoid the hard exit.
    try {
      fAinv.InvertFast(&det);
    }
    catch(...) {
      std::cerr << __FILE__ << ":" << __LINE__
                << " Matrix inversion failed" << std::endl;
      return false;
    } // end of inverting the matrix
    
    if (det<=0.0) {
      std::cerr << __FILE__ << ":" << __LINE__ 
                << " Matrix inversion failed" << std::endl;
      return false;
    }
    return true;
  }
  
  //......................................................................
  //
  // Equation 35
  //
  void Lutz::Calcq() 
  {
    unsigned int j;
    //
    // Solve according to equation 35. Don't bother with the second term
    // which is by construction zero.
    //
    TVectorT<double> q(fXI);
    q *= fD;
    q *= fAinv;
    
    //
    // Unpack the resutls
    //
    fa = q[0];
    fb = q[1];
    for (j=0; j<fN; ++j) fAlpha[j] = q[2+j];
  }
  
  //......................................................................
  
  void Lutz::Calcx() 
  {
    unsigned int i, j;
    for (i=0; i<fn; ++i) {
      fx[i] = fa + fb*fz[i];
      for (j=0; j<fN; ++j) {
        if (fZeta[i][j]>0.0) fx[i] += fAlpha[j]*fZeta[i][j];
      }
    }
  }
  
  //......................................................................
  
  double Lutz::X(double z) const
  {
    unsigned int j;
    double zeta;
    double x = fa + fb*z;
    for (j=0; j<fN; ++j) {
      zeta = z-fZ[j];
      if (zeta>0.0) x += fAlpha[j]*zeta;
    }
    return x;
  }
  
  //......................................................................
  
  void Lutz::Fit(double *a, double* b, double* alpha, double* chi2) 
  {
    bool ok;
    unsigned int j;
    //
    // Make the fit
    //
    this->CalcZeta();
    this->CalcD();
    ok = this->CalcAinv();
    if (ok) {
      this->Calcq();
      this->Calcx();    
    }
    else this->NullResult();
    //
    // Unpack the results
    //
    if (a!=0) *a = fa;
    if (b!=0) *b = fb;
    if (alpha!=0) {
      for (j=0; j<fN; ++j) alpha[j] = fAlpha[j];
    }
    if (chi2!=0) {
      *chi2 = this->Chi2XI() + this->Chi2Beta();
    }
  }
  
  //......................................................................
  
  double Lutz::Chi2XIi(unsigned int i) const
  {
    double chi;
    
    chi  = fXI[i]-fx[i];
    chi *= chi;
    chi *= fwx[i];
    
    return chi;
  }
  
  //......................................................................
  
  double Lutz::Chi2BetaJ(unsigned int j) const
  {
    double chi;
    
    chi = 0.0-fAlpha[j];
    chi *= chi;
    chi *= fwSJ[j];
    
    return chi;
  }
  
  //......................................................................
  
  double Lutz::Chi2XI() const
  {
    unsigned int i;
    double sum = 0.0;
    for (i=0; i<fn; ++i) sum += this->Chi2XIi(i);
    return sum;
  }
  
  //......................................................................
  
  double Lutz::Chi2Beta() const
  {
    unsigned int j;
    double sum = 0.0;
    for (j=0; j<fN; ++j) sum += this->Chi2BetaJ(j);
    return sum;
  }

  //....................................................................

  void Lutz::NullResult()
  {
    unsigned int i, j;
    fa = 0.0;
    fb = 0.0;
    for (i=0; i<fn; ++i) fx[i]     = 0.0;
    for (j=0; j<fN; ++j) fAlpha[j] = 0.0;
  }
}
  
////////////////////////////////////////////////////////////////////////