/home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Muc/MucRecEvent/MucRecEvent-00-02-52/src/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 <iostream>
#include <math.h>
#include "MucGeomSvc/MucGeomSvc.h"
#include "MucRecEvent/MucRecLineFit.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;
}   

---