BOSS 7.1.1
BESIII Offline Software System
Loading...
Searching...
No Matches
MucRecLineFit Class Reference

#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)
 

Detailed Description

Definition at line 15 of file MucRecLineFit.h.

Constructor & Destructor Documentation

◆ MucRecLineFit()

MucRecLineFit::MucRecLineFit ( )

Constructor.

Definition at line 19 of file MucRecLineFit.cxx.

20{
21 // Constructor.
22}

◆ ~MucRecLineFit()

MucRecLineFit::~MucRecLineFit ( )

Destructor.

Definition at line 24 of file MucRecLineFit.cxx.

25{
26 // Destructor.
27}

Member Function Documentation

◆ LineFit() [1/2]

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.

186{
187 double sum, sx, sy, sxx, sxy, syy, det;
188 float chi;
189
190 /* The variable "i" is declared 8 bytes long to avoid triggering an
191 apparent alignment bug in optimized code produced by gcc-2.95.3.
192 The bug doesn't seem to be present in gcc-3.2. */
193 long i;
194
195 chi = 99999999.0;
196
197 /* N must be >= 2 for this guy to work */
198 if (n < 2)
199 {
200 cout << "utiLineFit-W: Too few points for line fit \n" << endl;
201 return -1;
202 }
203
204 /* Initialization */
205 sum = 0.0;
206 sx = 0.0;
207 sy = 0.0;
208 sxx = 0.0;
209 sxy = 0.0;
210 syy = 0.0;
211
212 /* Find sum , sumx ,sumy, sumxx, sumxy */
213 for (i = 0; i < n; i++)
214 {
215 //cout<<"x: "<<x[i]<<" y: "<<y[i]<<" w: "<<w[i]<<endl;
216 sum = sum + w[i];
217 sx = sx + w[i] * x[i];
218 sy = sy + w[i] * y[i];
219 sxx = sxx + w[i] * x[i] * x[i];
220 sxy = sxy + w[i] * x[i] * y[i];
221 syy = syy + w[i] * y[i] * y[i];
222 }
223
224 det = sum*sxx-sx*sx;
225 if (fabs(det) < 1.0e-20)
226 {
227 *a = 1.0e20;
228 *b = x[0];
229 *chisq = 0.0;
230 *siga = 0.0;
231 *sigb = 0.0;
232
233 return 0;
234 }
235
236 /* compute the best fitted parameters A,B */
237 *a = (sum*sxy-sx*sy)/det;
238 *b = (sy*sxx-sxy*sx)/det;
239
240 //cout<<"a: "<<*a<<" b: "<<*b<<endl;
241 /* calculate chi-square */
242 chi = 0.0;
243 for (i=0; i<n; i++)
244 {
245 chi = chi+(w[i])*((y[i])-*a*(x[i])-*b)*
246 ((y[i])-*a*(x[i])-*b);
247 }
248
249 *siga = sqrt(fabs(sum/det));
250 *sigb = sqrt(fabs(sxx/det));
251
252
253 *chisq = chi;
254 return 0;
255}
const Int_t n
Double_t x[10]
double w
double y[1000]
const double b
Definition slope.cxx:9

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

◆ LineFit() [2/2]

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.

79{
80 if(part>2||part<0) {
81 cout<<"BAD part id from MucRecLineFit"<<endl;
82 return -1;
83 }
84
85 int status = -1;
86 float c,d,sigc,sigd;
87
88 if(part != 1 ){
89 for(int i = 0; i < n; i++){
90 float gap0z = (MucGeoGeneral::Instance()->GetGap(part, seg, 0)->GetCenter()).z();
91 //cout<<"in MucRecLineFit change x from "<<x[i]<<" to ";
92 x[i] -= gap0z;
93 //cout<<x[i]<<endl;
94 }
95 status = LineFit(x,y,w,n,a,b,chisq,siga,sigb);
96 //cout<<"in MucRecLineFit sigb = "<<*sigb<<endl;
97 }
98 else{
99 if(orient == 0){ //need change x and y because x has error now.
100
101 for(int i = 0; i < n; i++){
102 float gap0r = (MucGeoGeneral::Instance()->GetGap(part,0 , 0)->GetCenter()).x();
103 y[i] -= gap0r;
104 }
105 status = LineFit(y,x,w,n,&c,&d,chisq,&sigc,&sigd);
106 //need to recalculate a,b,siga,sigb now
107 // y=ax+b; x=cy+d
108 // a = 1/c; b = -d/c
109 // siga^2 = a^2/c^2 * sigc^2
110 // sigb^2 = b^2/c^2 * sigc^2 + b^2/d^2 * sigd^2;
111
112 *a = 1/c;
113 *b = -1 * d / c;
114 *siga = 1/c/c * sigc;
115 //*sigb = sqrt(fabs(1/c/c*sigd*sigd + d*d/c/c/c/c*sigc*sigc));
116 *sigb = sigd; ////////////
117 }
118 else{ // barrel... seg0,4 need not change; seg2,6 change x<->y; other seg need more complicate operation!
119 if(seg == 0||seg == 4){
120 float gap0x = (MucGeoGeneral::Instance()->GetGap(part,seg , 0)->GetCenter()).x();
121 for(int i = 0; i < n; i++){
122 x[i] -= gap0x;
123 }
124 status = LineFit(x,y,w,n,a,b,chisq,siga,sigb);
125 }
126 else if(seg == 2||seg == 6){
127 float gap0y = (MucGeoGeneral::Instance()->GetGap(part,seg , 0)->GetCenter()).y();
128 for(int i = 0; i < n; i++){
129 y[i] -= gap0y;
130 }
131 status = LineFit(y,x,w,n,&c,&d,chisq,&sigc,&sigd);
132 //need to recalculate a,b,siga,sigb now
133 *a = 1/c;
134 *b = -1 * d / c;
135 *siga = 1/c/c * sigc;
136 //*sigb = sqrt(fabs(1/c/c*sigd*sigd + d*d/c/c/c/c*sigc*sigc));
137 *sigb = sigd;
138 }
139 else{
140 for(int i = 0; i < n; i++){
141 float temp = (y[i] + x[i])/sqrt(2.0);
142 y[i] = (y[i] - x[i])/sqrt(2.0);
143 x[i] = temp;
144 }
145
146 if(seg == 1 || seg == 5){
147 for(int i = 0; i < n; i++){
148 float gap0x = (MucGeoGeneral::Instance()->GetGap(part,seg-1 , 0)->GetCenter()).x();
149 x[i] -= gap0x;
150 }
151 status = LineFit(x,y,w,n,a,b,chisq,siga,sigb);
152 }
153 else if(seg == 3 || seg == 7){
154 for(int i = 0; i < n; i++){
155 float gap0y = (MucGeoGeneral::Instance()->GetGap(part,seg-1 , 0)->GetCenter()).y();
156 y[i] -= gap0y;
157 }
158 status = LineFit(y,x,w,n,a,b,chisq,siga,sigb);
159
160 }
161
162 }
163
164
165 }
166
167 }
168
169 return status;
170
171}
HepPoint3D GetCenter() const
Get gap center position in global coordinate.
Definition MucGeoGap.h:137
MucGeoGap * GetGap(const int part, const int seg, const int gap) const
Get a pointer to the gap identified by (part,seg,gap).
static MucGeoGeneral * Instance()
Get a pointer to the single instance of MucGeoGeneral.
int LineFit(float x[], float y[], float w[], int n, float *a, float *b, float *chisq, float *siga, float *sigb)

The documentation for this class was generated from the following files: