MucRecLineFit Class Reference

#include <MucRecLineFit.h>

List of all members.

Public Member Functions

 MucRecLineFit ()
 Constructor.
 ~MucRecLineFit ()
 Destructor.
int LineFit (float x[], float y[], float w[], int n, float *a, float *b, float *chisq, float *siga, float *sigb)
int 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)


Detailed Description

Definition at line 15 of file MucRecLineFit.h.


Constructor & Destructor Documentation

MucRecLineFit::MucRecLineFit (  ) 

Constructor.

Definition at line 19 of file MucRecLineFit.cxx.

00020 {
00021     // Constructor.
00022 }

MucRecLineFit::~MucRecLineFit (  ) 

Destructor.

Definition at line 24 of file MucRecLineFit.cxx.

00025 {
00026   // Destructor.
00027 }


Member Function Documentation

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 
)

Definition at line 67 of file MucRecLineFit.cxx.

References MucGeoGap::GetCenter(), MucGeoGeneral::GetGap(), genRecEmupikp::i, MucGeoGeneral::Instance(), LineFit(), subSeperate::temp, and x.

00079 {
00080   if(part>2||part<0)  {
00081     cout<<"BAD part id from MucRecLineFit"<<endl;
00082     return -1;
00083   }
00084 
00085   int status = -1;
00086   float c,d,sigc,sigd;
00087   
00088   if(part != 1 ){
00089     for(int i = 0; i < n; i++){
00090       float gap0z = (MucGeoGeneral::Instance()->GetGap(part, seg, 0)->GetCenter()).z();
00091       //cout<<"in MucRecLineFit change x from "<<x[i]<<" to ";
00092       x[i] -= gap0z;
00093       //cout<<x[i]<<endl;
00094     }
00095     status = LineFit(x,y,w,n,a,b,chisq,siga,sigb);
00096     //cout<<"in MucRecLineFit sigb = "<<*sigb<<endl;
00097   }
00098   else{
00099     if(orient == 0){  //need change x and y because x has error now.
00100 
00101       for(int i = 0; i < n; i++){
00102         float gap0r = (MucGeoGeneral::Instance()->GetGap(part,0 , 0)->GetCenter()).x();
00103         y[i] -= gap0r;
00104       }
00105       status = LineFit(y,x,w,n,&c,&d,chisq,&sigc,&sigd);
00106       //need to recalculate a,b,siga,sigb now
00107       // y=ax+b;  x=cy+d
00108       // a = 1/c; b = -d/c
00109       // siga^2 = a^2/c^2 * sigc^2
00110       // sigb^2 = b^2/c^2 * sigc^2 + b^2/d^2 * sigd^2;
00111 
00112       *a = 1/c;
00113       *b = -1 * d / c;
00114       *siga = 1/c/c * sigc;
00115       //*sigb = sqrt(fabs(1/c/c*sigd*sigd + d*d/c/c/c/c*sigc*sigc));
00116       *sigb = sigd;  
00117     }
00118     else{  // barrel... seg0,4 need not change; seg2,6 change x<->y; other seg need more complicate operation!
00119       if(seg == 0||seg == 4){
00120         float gap0x = (MucGeoGeneral::Instance()->GetGap(part,seg , 0)->GetCenter()).x();
00121         for(int i = 0; i < n; i++){
00122           x[i] -= gap0x;
00123         }
00124         status = LineFit(x,y,w,n,a,b,chisq,siga,sigb);
00125       }
00126       else if(seg == 2||seg == 6){
00127         float gap0y = (MucGeoGeneral::Instance()->GetGap(part,seg , 0)->GetCenter()).y();
00128         for(int i = 0; i < n; i++){
00129           y[i] -= gap0y;
00130         }   
00131         status = LineFit(y,x,w,n,&c,&d,chisq,&sigc,&sigd);
00132         //need to recalculate a,b,siga,sigb now
00133         *a = 1/c; 
00134         *b = -1 * d / c;
00135         *siga = 1/c/c * sigc;
00136         //*sigb = sqrt(fabs(1/c/c*sigd*sigd + d*d/c/c/c/c*sigc*sigc));
00137         *sigb = sigd;
00138       }
00139       else{
00140         for(int i = 0; i < n; i++){
00141           float temp = (y[i] + x[i])/sqrt(2.0);
00142           y[i] = (y[i] - x[i])/sqrt(2.0);
00143           x[i] = temp;
00144         }
00145 
00146         if(seg == 1 || seg == 5){
00147           for(int i = 0; i < n; i++){
00148             float gap0x = (MucGeoGeneral::Instance()->GetGap(part,seg-1 , 0)->GetCenter()).x();
00149             x[i] -= gap0x;
00150           }
00151           status = LineFit(x,y,w,n,a,b,chisq,siga,sigb);
00152         }
00153        else if(seg == 3 || seg == 7){
00154          for(int i = 0; i < n; i++){
00155            float gap0y = (MucGeoGeneral::Instance()->GetGap(part,seg-1 , 0)->GetCenter()).y();
00156            y[i] -= gap0y;
00157          }
00158          status = LineFit(y,x,w,n,a,b,chisq,siga,sigb);
00159 
00160        }
00161 
00162       }
00163 
00164 
00165     }
00166     
00167   }
00168 
00169   return status;
00170 
00171 }

int MucRecLineFit::LineFit ( float  x[],
float  y[],
float  w[],
int  n,
float *  a,
float *  b,
float *  chisq,
float *  siga,
float *  sigb 
)

Definition at line 176 of file MucRecLineFit.cxx.

References genRecEmupikp::i.

Referenced by MucRec2DRoad::Fit(), LineFit(), MucRec2DRoad::SimpleFit(), and MucRec2DRoad::SimpleFitNoCurrentGap().

00186 {
00187   double sum, sx, sy, sxx, sxy, syy, det;
00188   float chi;
00189   
00190   /* The variable "i" is declared 8 bytes long to avoid triggering an
00191      apparent alignment bug in optimized code produced by gcc-2.95.3.
00192      The bug doesn't seem to be present in gcc-3.2. */
00193   long i; 
00194 
00195   chi = 99999999.0;
00196 
00197   /* N must be >= 2 for this guy to work */
00198   if (n < 2) 
00199     {
00200       cout << "utiLineFit-W: Too few points for line fit \n" << endl;
00201       return -1;
00202     }
00203 
00204   /* Initialization */
00205   sum = 0.0;
00206   sx  = 0.0;
00207   sy  = 0.0;
00208   sxx = 0.0;
00209   sxy = 0.0;
00210   syy = 0.0;
00211 
00212   /* Find sum , sumx ,sumy, sumxx, sumxy */
00213   for (i = 0; i < n; i++)
00214     {
00215       //cout<<"x: "<<x[i]<<"  y: "<<y[i]<<"  w: "<<w[i]<<endl;
00216       sum = sum + w[i];
00217       sx  = sx + w[i] * x[i];
00218       sy  = sy + w[i] * y[i];
00219       sxx = sxx + w[i] * x[i] * x[i];
00220       sxy = sxy + w[i] * x[i] * y[i];
00221       syy = syy + w[i] * y[i] * y[i];
00222     }
00223 
00224   det = sum*sxx-sx*sx;
00225   if (fabs(det) < 1.0e-20) 
00226     {
00227       *a = 1.0e20;
00228       *b = x[0];
00229       *chisq = 0.0;
00230       *siga  = 0.0;
00231       *sigb  = 0.0;
00232 
00233       return 0;
00234     }
00235 
00236   /* compute the best fitted parameters A,B */
00237   *a = (sum*sxy-sx*sy)/det;
00238   *b = (sy*sxx-sxy*sx)/det;
00239     
00240   //cout<<"a: "<<*a<<"  b: "<<*b<<endl;
00241   /* calculate chi-square */
00242   chi = 0.0;
00243   for (i=0; i<n; i++) 
00244     {
00245       chi = chi+(w[i])*((y[i])-*a*(x[i])-*b)*
00246         ((y[i])-*a*(x[i])-*b);
00247     }
00248   
00249   *siga = sqrt(fabs(sum/det));
00250   *sigb = sqrt(fabs(sxx/det));
00251   
00252   
00253   *chisq = chi;
00254   return 0;
00255 }   


Generated on Tue Nov 29 23:20:30 2016 for BOSS_7.0.2 by  doxygen 1.4.7