#include <MucRecLineFit.h>
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) |
Definition at line 15 of file MucRecLineFit.h.
MucRecLineFit::MucRecLineFit | ( | ) |
MucRecLineFit::~MucRecLineFit | ( | ) |
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 }