MucRecLineFit.cxx

Go to the documentation of this file.

//$id$
//
//$log$
 
/*
 *    2003/12/24   Zhengyun You     Peking University
 *
 *    2004/09/12   Zhengyun You     Peking University
 *                 transplanted to Gaudi framework
 */
 
using namespace std;
 
#include "MucRecEvent/MucRecLineFit.h"
#include "MucGeomSvc/MucGeoGeneral.h"
#include <iostream>
#include <math.h>
 
MucRecLineFit::MucRecLineFit() {
  // Constructor.
}
 
MucRecLineFit::~MucRecLineFit() {
  // Destructor.
}
 
/* --------------------------------------------------------------
 
 utiLineFit - a least squares straight line fitting program
 
 DESCRIPTION:
 Performs weighted least squares fit to line (Y  =  A*X  +  B )
 using linear regression.
 
 INPUT ARGUMENTS:
  X(N),Y(N),W(N) - Input values and weights (N=1,2,3...)
  N              - Number of pairs of data points,
  X              - Array of data points in X-axis,
  Y              - Array of data points in Y-axis,
  W              - Array of weights for data points,
 
 OUTPUT ARGUMENTS:
  B              - Y intercept of best fitted straight line,
  SIGB           - Standard deviation of B,
  A              - Slope of fitted straight line,
  SIGA           - Standard deviation of A,
  CHISQ          - Chi-square
  LSWLF          - = 0  return without error
                   = -1 got some fitting problems
 
 AUTHOR:
 Jawluen Tang, Physics department, UT-Austin
 J. T. Mitchell - adapted for PHENIX use. Converted to C.
 
 USAGE:
 utiLineFit(X,Y,W,N,&A,&B,&CHISQ,&SIGA,&SIGB)
 The arrays must start counting at element 1.
 
--------------------------------------------------------------- */
 
// add (part,seg,orient) to choose better coordinator.
 
int MucRecLineFit::LineFit( float x[], float y[], float w[], int part, int seg, int orient,
                            int n, float* a, float* b, float* chisq, float* siga,
                            float* sigb ) {
  if ( part > 2 || part < 0 )
  {
    cout << "BAD part id from MucRecLineFit" << endl;
    return -1;
  }
 
  int   status = -1;
  float c, d, sigc, sigd;
 
  if ( part != 1 )
  {
    for ( int i = 0; i < n; i++ )
    {
      float gap0z = ( MucGeoGeneral::Instance()->GetGap( part, seg, 0 )->GetCenter() ).z();
      // cout<<"in MucRecLineFit change x from "<<x[i]<<" to ";
      x[i] -= gap0z;
      // cout<<x[i]<<endl;
    }
    status = LineFit( x, y, w, n, a, b, chisq, siga, sigb );
    // cout<<"in MucRecLineFit sigb = "<<*sigb<<endl;
  }
  else
  {
    if ( orient == 0 )
    { // need change x and y because x has error now.
 
      for ( int i = 0; i < n; i++ )
      {
        float gap0r = ( MucGeoGeneral::Instance()->GetGap( part, 0, 0 )->GetCenter() ).x();
        y[i] -= gap0r;
      }
      status = LineFit( y, x, w, n, &c, &d, chisq, &sigc, &sigd );
      // need to recalculate a,b,siga,sigb now
      //  y=ax+b;  x=cy+d
      //  a = 1/c; b = -d/c
      //  siga^2 = a^2/c^2 * sigc^2
      //  sigb^2 = b^2/c^2 * sigc^2 + b^2/d^2 * sigd^2;
 
      *a    = 1 / c;
      *b    = -1 * d / c;
      *siga = 1 / c / c * sigc;
      //*sigb = sqrt(fabs(1/c/c*sigd*sigd + d*d/c/c/c/c*sigc*sigc));
      *sigb = sigd; ////////////
    }
    else
    { // barrel... seg0,4 need not change; seg2,6 change x<->y; other seg need more complicate
      // operation!
      if ( seg == 0 || seg == 4 )
      {
        float gap0x = ( MucGeoGeneral::Instance()->GetGap( part, seg, 0 )->GetCenter() ).x();
        for ( int i = 0; i < n; i++ ) { x[i] -= gap0x; }
        status = LineFit( x, y, w, n, a, b, chisq, siga, sigb );
      }
      else if ( seg == 2 || seg == 6 )
      {
        float gap0y = ( MucGeoGeneral::Instance()->GetGap( part, seg, 0 )->GetCenter() ).y();
        for ( int i = 0; i < n; i++ ) { y[i] -= gap0y; }
        status = LineFit( y, x, w, n, &c, &d, chisq, &sigc, &sigd );
        // need to recalculate a,b,siga,sigb now
        *a    = 1 / c;
        *b    = -1 * d / c;
        *siga = 1 / c / c * sigc;
        //*sigb = sqrt(fabs(1/c/c*sigd*sigd + d*d/c/c/c/c*sigc*sigc));
        *sigb = sigd;
      }
      else
      {
        for ( int i = 0; i < n; i++ )
        {
          float temp = ( y[i] + x[i] ) / sqrt( 2.0 );
          y[i]       = ( y[i] - x[i] ) / sqrt( 2.0 );
          x[i]       = temp;
        }
 
        if ( seg == 1 || seg == 5 )
        {
          for ( int i = 0; i < n; i++ )
          {
            float gap0x =
                ( MucGeoGeneral::Instance()->GetGap( part, seg - 1, 0 )->GetCenter() ).x();
            x[i] -= gap0x;
          }
          status = LineFit( x, y, w, n, a, b, chisq, siga, sigb );
        }
        else if ( seg == 3 || seg == 7 )
        {
          for ( int i = 0; i < n; i++ )
          {
            float gap0y =
                ( MucGeoGeneral::Instance()->GetGap( part, seg - 1, 0 )->GetCenter() ).y();
            y[i] -= gap0y;
          }
          status = LineFit( y, x, w, n, a, b, chisq, siga, sigb );
        }
      }
    }
  }
 
  return status;
}
 
int MucRecLineFit::LineFit( float x[], float y[], float w[], int n, float* a, float* b,
                            float* chisq, float* siga, float* sigb )
 
{
  double sum, sx, sy, sxx, sxy, syy, det;
  float  chi;
 
  /* The variable "i" is declared 8 bytes long to avoid triggering an
     apparent alignment bug in optimized code produced by gcc-2.95.3.
     The bug doesn't seem to be present in gcc-3.2. */
  long i;
 
  chi = 99999999.0;
 
  /* N must be >= 2 for this guy to work */
  if ( n < 2 )
  {
    cout << "utiLineFit-W: Too few points for line fit \n" << endl;
    return -1;
  }
 
  /* Initialization */
  sum = 0.0;
  sx  = 0.0;
  sy  = 0.0;
  sxx = 0.0;
  sxy = 0.0;
  syy = 0.0;
 
  /* Find sum , sumx ,sumy, sumxx, sumxy */
  for ( i = 0; i < n; i++ )
  {
    // cout<<"x: "<<x[i]<<"  y: "<<y[i]<<"  w: "<<w[i]<<endl;
    sum = sum + w[i];
    sx  = sx + w[i] * x[i];
    sy  = sy + w[i] * y[i];
    sxx = sxx + w[i] * x[i] * x[i];
    sxy = sxy + w[i] * x[i] * y[i];
    syy = syy + w[i] * y[i] * y[i];
  }
 
  det = sum * sxx - sx * sx;
  if ( fabs( det ) < 1.0e-20 )
  {
    *a     = 1.0e20;
    *b     = x[0];
    *chisq = 0.0;
    *siga  = 0.0;
    *sigb  = 0.0;
 
    return 0;
  }
 
  /* compute the best fitted parameters A,B */
  *a = ( sum * sxy - sx * sy ) / det;
  *b = ( sy * sxx - sxy * sx ) / det;
 
  // cout<<"a: "<<*a<<"  b: "<<*b<<endl;
  /* calculate chi-square */
  chi = 0.0;
  for ( i = 0; i < n; i++ )
  {
    chi =
        chi + ( w[i] ) * ( ( y[i] ) - *a * ( x[i] ) - *b ) * ( ( y[i] ) - *a * ( x[i] ) - *b );
  }
 
  *siga = sqrt( fabs( sum / det ) );
  *sigb = sqrt( fabs( sxx / det ) );
 
  *chisq = chi;
  return 0;
}