File: | libraries/TRACKING/DHelicalFit.cc |
Location: | line 1165, column 8 |
Description: | Value stored to 'dzdphi' during its initialization is never read |
1 | // Set of routines for fitting tracks in both the cdc and fdc assuming a |
2 | // uniform magnetic field. The track therefore is assumed to follow a helical |
3 | // path. |
4 | |
5 | #include <iostream> |
6 | #include <algorithm> |
7 | using namespace std; |
8 | |
9 | #include "DANA/DApplication.h" |
10 | #include <TDecompLU.h> |
11 | #include <math.h> |
12 | |
13 | #include "DHelicalFit.h" |
14 | #define qBr2p0.003 0.003 // conversion for converting q*B*r to GeV/c |
15 | |
16 | #define ONE_THIRD0.33333333333333333 0.33333333333333333 |
17 | #define SQRT31.73205080756887719 1.73205080756887719 |
18 | #define EPS1e-8 1e-8 |
19 | #ifndef M_TWO_PI6.28318530717958647692 |
20 | #define M_TWO_PI6.28318530717958647692 6.28318530717958647692 |
21 | #endif |
22 | |
23 | // The following is for sorting hits by z |
24 | class DHFHitLessThanZ{ |
25 | public: |
26 | bool operator()(DHFHit_t* const &a, DHFHit_t* const &b) const { |
27 | return a->z < b->z; |
28 | } |
29 | }; |
30 | |
31 | inline bool DHFHitLessThanZ_C(DHFHit_t* const &a, DHFHit_t* const &b) { |
32 | return a->z < b->z; |
33 | } |
34 | inline bool RiemannFit_hit_cmp(DHFHit_t *a,DHFHit_t *b){ |
35 | return (a->z>b->z); |
36 | } |
37 | |
38 | |
39 | |
40 | //----------------- |
41 | // DHelicalFit (Constructor) |
42 | //----------------- |
43 | DHelicalFit::DHelicalFit(void) |
44 | { |
45 | x0 = y0 = r0 = tanl = z_vertex = p_trans = phi= theta = q = p = dzdphi =0; |
46 | chisq = 0; |
47 | ndof=0; |
48 | chisq_source = NOFIT; |
49 | bfield = NULL__null; |
50 | |
51 | // For Riemann Fit |
52 | CovR_=NULL__null; |
53 | CovRPhi_=NULL__null; |
54 | c_origin=0.; |
55 | normal.SetXYZ(0.,0.,0.); |
56 | } |
57 | |
58 | //----------------- |
59 | // DHelicalFit (Constructor) |
60 | //----------------- |
61 | DHelicalFit::DHelicalFit(const DHelicalFit &fit) |
62 | { |
63 | Copy(fit); |
64 | } |
65 | //----------------- |
66 | // Copy |
67 | //----------------- |
68 | void DHelicalFit::Copy(const DHelicalFit &fit) |
69 | { |
70 | x0 = fit.x0; |
71 | y0 = fit.y0; |
72 | r0 = fit.r0; |
73 | q = fit.q; |
74 | p = fit.p; |
75 | tanl=fit.tanl; |
76 | p_trans = fit.p_trans; |
77 | phi = fit.phi; |
78 | theta = fit.theta; |
79 | z_vertex = fit.z_vertex; |
80 | chisq = fit.chisq; |
81 | ndof=fit.ndof; |
82 | dzdphi = fit.dzdphi; |
83 | chisq_source = fit.chisq_source; |
84 | bfield = fit.GetMagneticFieldMap(); |
85 | Bz_avg = fit.GetBzAvg(); |
86 | z_mean = fit.GetZMean(); |
87 | phi_mean = fit.GetPhiMean(); |
88 | |
89 | const vector<DHFHit_t*> myhits = fit.GetHits(); |
90 | for(unsigned int i=0; i<myhits.size(); i++){ |
91 | DHFHit_t *a = new DHFHit_t; |
92 | *a = *myhits[i]; |
93 | hits.push_back(a); |
94 | } |
95 | if (fit.CovR_!=NULL__null){ |
96 | CovR_ = new DMatrix(hits.size(),hits.size()); |
97 | for (unsigned int i=0;i<hits.size();i++) |
98 | for (unsigned int j=0;j<hits.size();j++) |
99 | CovR_->operator()(i,j)=fit.CovR_->operator()(i,j); |
100 | } |
101 | else CovR_=NULL__null; |
102 | if (fit.CovRPhi_!=NULL__null){ |
103 | CovRPhi_ = new DMatrix(hits.size(),hits.size()); |
104 | for (unsigned int i=0;i<hits.size();i++) |
105 | for (unsigned int j=0;j<hits.size();j++) |
106 | CovRPhi_->operator()(i,j)=fit.CovRPhi_->operator()(i,j); |
107 | } |
108 | else CovRPhi_=NULL__null; |
109 | } |
110 | |
111 | //----------------- |
112 | // operator= (Assignment operator) |
113 | //----------------- |
114 | DHelicalFit& DHelicalFit::operator=(const DHelicalFit& fit) |
115 | { |
116 | if(this == &fit)return *this; |
117 | Copy(fit); |
118 | |
119 | return *this; |
120 | } |
121 | |
122 | //----------------- |
123 | // DHelicalFit (Destructor) |
124 | //----------------- |
125 | DHelicalFit::~DHelicalFit() |
126 | { |
127 | Clear(); |
128 | } |
129 | //----------------- |
130 | // AddHit -- add an FDC hit to the list |
131 | //----------------- |
132 | jerror_t DHelicalFit::AddHit(const DFDCPseudo *fdchit){ |
133 | DHFHit_t *hit = new DHFHit_t; |
134 | hit->x = fdchit->xy.X(); |
135 | hit->y = fdchit->xy.Y(); |
136 | hit->z = fdchit->wire->origin.z(); |
137 | hit->covx=fdchit->covxx; |
138 | hit->covy=fdchit->covyy; |
139 | hit->covxy=fdchit->covxy; |
140 | hit->chisq = 0.0; |
141 | hits.push_back(hit); |
142 | |
143 | return NOERROR; |
144 | |
145 | } |
146 | |
147 | |
148 | //----------------- |
149 | // AddHit |
150 | //----------------- |
151 | jerror_t DHelicalFit::AddHit(float r, float phi, float z) |
152 | { |
153 | /// Add a hit to the list of hits using cylindrical coordinates |
154 | |
155 | return AddHitXYZ(r*cos(phi), r*sin(phi), z); |
156 | } |
157 | |
158 | //----------------- |
159 | // AddHitXYZ |
160 | //----------------- |
161 | jerror_t DHelicalFit::AddHitXYZ(float x, float y, float z) |
162 | { |
163 | /// Add a hit to the list of hits using Cartesian coordinates |
164 | DHFHit_t *hit = new DHFHit_t; |
165 | hit->x = x; |
166 | hit->y = y; |
167 | hit->z = z; |
168 | hit->covx=1.; |
169 | hit->covy=1.; |
170 | hit->covxy=0.; |
171 | hit->chisq = 0.0; |
172 | hits.push_back(hit); |
173 | |
174 | return NOERROR; |
175 | } |
176 | |
177 | // Add a hit to the list of hits using Cartesian coordinates |
178 | jerror_t DHelicalFit::AddHitXYZ(float x,float y, float z,float covx, |
179 | float covy, float covxy){ |
180 | DHFHit_t *hit = new DHFHit_t; |
181 | hit->x = x; |
182 | hit->y = y; |
183 | hit->z = z; |
184 | hit->covx=covx; |
185 | hit->covy=covy; |
186 | hit->covxy=covxy; |
187 | hits.push_back(hit); |
188 | |
189 | return NOERROR; |
190 | } |
191 | |
192 | |
193 | //----------------- |
194 | // PruneHit |
195 | //----------------- |
196 | jerror_t DHelicalFit::PruneHit(int idx) |
197 | { |
198 | /// Remove the hit specified by idx from the list |
199 | /// of hits. The value of idx can be anywhere from |
200 | /// 0 to GetNhits()-1. |
201 | if(idx<0 || idx>=(int)hits.size())return VALUE_OUT_OF_RANGE; |
202 | |
203 | delete hits[idx]; |
204 | hits.erase(hits.begin() + idx); |
205 | |
206 | return NOERROR; |
207 | } |
208 | |
209 | //----------------- |
210 | // Clear |
211 | //----------------- |
212 | jerror_t DHelicalFit::Clear(void) |
213 | { |
214 | /// Remove covariance matrices |
215 | if (CovR_!=NULL__null) delete CovR_; |
216 | if (CovRPhi_!=NULL__null) delete CovRPhi_; |
217 | |
218 | /// Remove all hits |
219 | for(unsigned int i=0; i<hits.size(); i++)delete hits[i]; |
220 | hits.clear(); |
221 | |
222 | return NOERROR; |
223 | } |
224 | |
225 | //----------------- |
226 | // PrintChiSqVector |
227 | //----------------- |
228 | jerror_t DHelicalFit::PrintChiSqVector(void) const |
229 | { |
230 | /// Dump the latest chi-squared vector to the screen. |
231 | /// This prints the individual hits' chi-squared |
232 | /// contributions in the order in which the hits were |
233 | /// added. See PruneHits() for more detail. |
234 | |
235 | cout<<"Chisq vector from DHelicalFit: (source="; |
236 | switch(chisq_source){ |
237 | case NOFIT: cout<<"NOFIT"; break; |
238 | case CIRCLE: cout<<"CIRCLE"; break; |
239 | case TRACK: cout<<"TRACK"; break; |
240 | }; |
241 | cout<<")"<<endl; |
242 | cout<<"----------------------------"<<endl; |
243 | |
244 | for(unsigned int i=0;i<hits.size();i++){ |
245 | cout<<i<<" "<<hits[i]->chisq<<endl; |
246 | } |
247 | cout<<"Total: "<<chisq<<endl<<endl; |
248 | |
249 | return NOERROR; |
250 | } |
251 | |
252 | //----------------- |
253 | // FitCircle |
254 | //----------------- |
255 | jerror_t DHelicalFit::FitCircle(void) |
256 | { |
257 | /// Fit the current set of hits to a circle |
258 | /// |
259 | /// This takes the hits which have been added thus far via one |
260 | /// of the AddHit() methods and fits them to a circle. |
261 | /// The alogorithm used here calculates the parameters directly using |
262 | /// a technique very much like linear regression. The key assumptions |
263 | /// are: |
264 | /// 1. The magnetic field is uniform and along z so that the projection |
265 | /// of the track onto the X/Y plane will fall on a circle |
266 | /// (this also implies no multiple-scattering or energy loss) |
267 | /// 2. The vertex is at the target center (i.e. 0,0 in the coordinate |
268 | /// system of the points passed to us. |
269 | /// |
270 | /// IMPORTANT: The value of phi which results from this assumes |
271 | /// the particle was positively charged. If the particle was |
272 | /// negatively charged, then phi will be 180 degrees off. To |
273 | /// correct this, one needs z-coordinate information to determine |
274 | /// the sign of the charge. |
275 | /// |
276 | /// ALSO IMPORTANT: This assumes a charge of +1 for the particle. If |
277 | /// the particle actually has a charge of +2, then the resulting |
278 | /// value of p_trans will be half of what it should be. |
279 | |
280 | float alpha=0.0, beta=0.0, gamma=0.0, deltax=0.0, deltay=0.0; |
281 | chisq_source = NOFIT; // in case we return early |
282 | ndof=hits.size()-2; |
283 | |
284 | // Loop over hits to calculate alpha, beta, gamma, and delta |
285 | // if a magnetic field map was given, use it to find average Z B-field |
286 | DHFHit_t *a = NULL__null; |
287 | for(unsigned int i=0;i<hits.size();i++){ |
288 | a = hits[i]; |
289 | float x=a->x; |
290 | float y=a->y; |
291 | alpha += x*x; |
292 | beta += y*y; |
293 | gamma += x*y; |
294 | deltax += x*(x*x+y*y)/2.0; |
295 | deltay += y*(x*x+y*y)/2.0; |
296 | } |
297 | |
298 | // Calculate x0,y0 - the center of the circle |
299 | double denom = alpha*beta-gamma*gamma; |
300 | if(fabs(denom)<1.0E-20)return UNRECOVERABLE_ERROR; |
301 | x0 = (deltax*beta-deltay*gamma)/denom; |
302 | //y0 = (deltay-gamma*x0)/beta; |
303 | y0 = (deltay*alpha-deltax*gamma)/denom; |
304 | |
305 | // Calculate the phi angle traversed by the track from the |
306 | // origin to the last hit. NOTE: this can be off by a multiple |
307 | // of 2pi! |
308 | double delta_phi=0.0; |
309 | if(a){ // a should be pointer to last hit from above loop |
310 | delta_phi = atan2(a->y-y0, a->x-x0); |
311 | if(delta_phi<0.0)delta_phi += M_TWO_PI6.28318530717958647692; |
312 | } |
313 | |
314 | // Momentum depends on magnetic field. If bfield has been |
315 | // set, we should use it to determine an average value of Bz |
316 | // for this track. Otherwise, assume -2T. |
317 | // Also assume a singly charged track (i.e. q=+/-1) |
318 | // The sign of the charge will be determined below. |
319 | Bz_avg=-2.0; |
320 | q = +1.0; |
321 | r0 = sqrt(x0*x0 + y0*y0); |
322 | p_trans = q*Bz_avg*r0*qBr2p0.003; // qBr2p converts to GeV/c |
323 | phi = atan2(y0,x0) - M_PI_21.57079632679489661923; |
324 | if(p_trans<0.0){ |
325 | p_trans = -p_trans; |
326 | } |
327 | if(phi<0)phi+=M_TWO_PI6.28318530717958647692; |
328 | if(phi>=M_TWO_PI6.28318530717958647692)phi-=M_TWO_PI6.28318530717958647692; |
329 | |
330 | // Calculate the chisq |
331 | ChisqCircle(); |
332 | chisq_source = CIRCLE; |
333 | |
334 | return NOERROR; |
335 | } |
336 | |
337 | //----------------- |
338 | // ChisqCircle |
339 | //----------------- |
340 | double DHelicalFit::ChisqCircle(void) |
341 | { |
342 | /// Calculate the chisq for the hits and current circle |
343 | /// parameters. |
344 | /// NOTE: This does not return the chi2/dof, just the |
345 | /// raw chi2 with errs set to 1.0 |
346 | chisq = 0.0; |
347 | for(unsigned int i=0;i<hits.size();i++){ |
348 | DHFHit_t *a = hits[i]; |
349 | float x = a->x - x0; |
350 | float y = a->y - y0; |
351 | float c = sqrt(x*x + y*y) - r0; |
352 | c *= c; |
353 | a->chisq = c; |
354 | chisq+=c; |
355 | } |
356 | |
357 | // Do NOT divide by DOF |
358 | return chisq; |
359 | } |
360 | |
361 | //----------------- |
362 | // FitCircleRiemann |
363 | //----------------- |
364 | jerror_t DHelicalFit::FitCircleRiemann(float z_vertex,float BeamRMS) |
365 | { |
366 | chisq_source = NOFIT; // in case we return early |
367 | |
368 | // Fake point at origin |
369 | float beam_var=BeamRMS*BeamRMS; |
370 | AddHitXYZ(0.,0.,z_vertex,beam_var,beam_var,0.); |
371 | |
372 | jerror_t err = FitCircleRiemann(); |
373 | if(err!=NOERROR)return err; |
374 | |
375 | // Number of degrees of freedom |
376 | ndof=hits.size()-3; |
377 | |
378 | // Momentum depends on magnetic field. If bfield has been |
379 | // set, we should use it to determine an average value of Bz |
380 | // for this track. Otherwise, assume -2T. |
381 | // Also assume a singly charged track (i.e. q=+/-1) |
382 | // The sign of the charge will be determined below. |
383 | Bz_avg=-2.0; |
384 | q = +1.0; |
385 | p_trans = q*Bz_avg*r0*qBr2p0.003; // qBr2p converts to GeV/c |
386 | if(p_trans<0.0){ |
387 | p_trans = -p_trans; |
388 | } |
389 | phi=atan2(-x0,y0); |
390 | if(phi<0)phi+=M_TWO_PI6.28318530717958647692; |
391 | if(phi>=M_TWO_PI6.28318530717958647692)phi-=M_TWO_PI6.28318530717958647692; |
392 | |
393 | // Normal vector for plane intersecting Riemann surface |
394 | normal.SetXYZ(N[0],N[1],N[2]); |
395 | |
396 | return NOERROR; |
397 | } |
398 | |
399 | //----------------- |
400 | // FitCircleRiemann |
401 | //----------------- |
402 | jerror_t DHelicalFit::FitCircleRiemann(void){ |
403 | /// Riemann Circle fit: points on a circle in x,y project onto a plane cutting |
404 | /// the circular paraboloid surface described by (x,y,x^2+y^2). Therefore the |
405 | /// task of fitting points in (x,y) to a circle is transormed to the task of |
406 | /// fitting planes in (x,y, w=x^2+y^2) space |
407 | /// |
408 | DMatrix X(hits.size(),3); |
409 | DMatrix Xavg(1,3); |
410 | DMatrix A(3,3); |
411 | // vector of ones |
412 | DMatrix OnesT(1,hits.size()); |
413 | double W_sum=0.; |
414 | DMatrix W(hits.size(),hits.size()); |
415 | |
416 | // Make sure hit list is ordered in z |
417 | std::sort(hits.begin(),hits.end(),RiemannFit_hit_cmp); |
418 | |
419 | // Covariance matrix |
420 | DMatrix CRPhi(hits.size(),hits.size()); |
421 | if (CovRPhi_==NULL__null){ |
422 | CovRPhi_ = new DMatrix(hits.size(),hits.size()); |
423 | for (unsigned int i=0;i<hits.size();i++){ |
424 | double Phi=atan2(hits[i]->y,hits[i]->x); |
425 | double cosPhi=cos(Phi); |
426 | double sinPhi=sin(Phi); |
427 | double Phi_cosPhi_minus_sinPhi=Phi*cosPhi-sinPhi; |
428 | double Phi_sinPhi_plus_cosPhi=Phi*sinPhi+cosPhi; |
429 | CovRPhi_->operator()(i,i) |
430 | =Phi_cosPhi_minus_sinPhi*Phi_cosPhi_minus_sinPhi*hits[i]->covx |
431 | +Phi_sinPhi_plus_cosPhi*Phi_sinPhi_plus_cosPhi*hits[i]->covy |
432 | +2.*Phi_sinPhi_plus_cosPhi*Phi_cosPhi_minus_sinPhi*hits[i]->covxy; |
433 | } |
434 | } |
435 | for (unsigned int i=0;i<hits.size();i++) |
436 | for (unsigned int j=0;j<hits.size();j++) |
437 | CRPhi(i,j)=CovRPhi_->operator()(i, j); |
438 | |
439 | // The goal is to find the eigenvector corresponding to the smallest |
440 | // eigenvalue of the equation |
441 | // lambda=n^T (X^T W X - W_sum Xavg^T Xavg)n |
442 | // where n is the normal vector to the plane slicing the cylindrical |
443 | // paraboloid described by the parameterization (x,y,w=x^2+y^2), |
444 | // and W is the weight matrix, assumed for now to be diagonal. |
445 | // In the absence of multiple scattering, W_sum is the sum of all the |
446 | // diagonal elements in W. |
447 | |
448 | for (unsigned int i=0;i<hits.size();i++){ |
449 | X(i,0)=hits[i]->x; |
450 | X(i,1)=hits[i]->y; |
451 | X(i,2)=hits[i]->x*hits[i]->x+hits[i]->y*hits[i]->y; |
452 | OnesT(0,i)=1.; |
453 | W(i,i)=1./CRPhi(i,i); |
454 | W_sum+=W(i,i); |
455 | } |
456 | Xavg=(1./W_sum)*(OnesT*(W*X)); |
457 | |
458 | A=DMatrix(DMatrix::kTransposed,X)*(W*X) |
459 | -W_sum*(DMatrix(DMatrix::kTransposed,Xavg)*Xavg); |
460 | if(!A.IsValid())return UNRECOVERABLE_ERROR; |
461 | |
462 | // The characteristic equation is |
463 | // lambda^3+B2*lambda^2+lambda*B1+B0=0 |
464 | // |
465 | double B2=-(A(0,0)+A(1,1)+A(2,2)); |
466 | double B1=A(0,0)*A(1,1)-A(1,0)*A(0,1)+A(0,0)*A(2,2)-A(2,0)*A(0,2) |
467 | +A(1,1)*A(2,2)-A(2,1)*A(1,2); |
468 | double B0=-A.Determinant(); |
469 | if(B0==0 || !finite(B0))return UNRECOVERABLE_ERROR; |
470 | |
471 | // The roots of the cubic equation are given by |
472 | // lambda1= -B2/3 + S+T |
473 | // lambda2= -B2/3 - (S+T)/2 + i sqrt(3)/2. (S-T) |
474 | // lambda3= -B2/3 - (S+T)/2 - i sqrt(3)/2. (S-T) |
475 | // where we define some temporary variables: |
476 | // S= (R+sqrt(Q^3+R^2))^(1/3) |
477 | // T= (R-sqrt(Q^3+R^2))^(1/3) |
478 | // Q=(3*B1-B2^2)/9 |
479 | // R=(9*B2*B1-27*B0-2*B2^3)/54 |
480 | // sum=S+T; |
481 | // diff=i*(S-T) |
482 | // We divide Q and R by a safety factor to prevent multiplying together |
483 | // enormous numbers that cause unreliable results. |
484 | |
485 | long double Q=(3.*B1-B2*B2)/9.e4; |
486 | long double R=(9.*B2*B1-27.*B0-2.*B2*B2*B2)/54.e6; |
487 | long double Q1=Q*Q*Q+R*R; |
488 | if (Q1<0) Q1=sqrt(-Q1); |
489 | else{ |
490 | return VALUE_OUT_OF_RANGE; |
491 | } |
492 | |
493 | // DeMoivre's theorem for fractional powers of complex numbers: |
494 | // (r*(cos(theta)+i sin(theta)))^(p/q) |
495 | // = r^(p/q)*(cos(p*theta/q)+i sin(p*theta/q)) |
496 | // |
497 | //double temp=100.*pow(R*R+Q1*Q1,0.16666666666666666667); |
498 | long double temp=100.*sqrt(cbrt(R*R+Q1*Q1)); |
499 | long double theta1=ONE_THIRD0.33333333333333333*atan2(Q1,R); |
500 | long double sum_over_2=temp*cos(theta1); |
501 | long double diff_over_2=-temp*sin(theta1); |
502 | // Third root |
503 | long double lambda_min=-ONE_THIRD0.33333333333333333*B2-sum_over_2+SQRT31.73205080756887719*diff_over_2; |
504 | |
505 | // Calculate the (normal) eigenvector corresponding to the eigenvalue lambda |
506 | N[0]=1.; |
507 | N[1]=(A(1,0)*A(0,2)-(A(0,0)-lambda_min)*A(1,2)) |
508 | /(A(0,1)*A(2,1)-(A(1,1)-lambda_min)*A(0,2)); |
509 | N[2]=(A(2,0)*(A(1,1)-lambda_min)-A(1,0)*A(2,1)) |
510 | /(A(1,2)*A(2,1)-(A(2,2)-lambda_min)*(A(1,1)-lambda_min)); |
511 | |
512 | // Normalize: n1^2+n2^2+n3^2=1 |
513 | double denom=sqrt(N[0]*N[0]+N[1]*N[1]+N[2]*N[2]); |
514 | for (int i=0;i<3;i++){ |
515 | N[i]/=denom; |
516 | } |
517 | |
518 | // Distance to origin |
519 | c_origin=-(N[0]*Xavg(0,0)+N[1]*Xavg(0,1)+N[2]*Xavg(0,2)); |
520 | |
521 | // Center and radius of the circle |
522 | double one_over_2Nz=1./(2.*N[2]); |
523 | //x0=-N[0]/2./N[2]; |
524 | //y0=-N[1]/2./N[2]; |
525 | //r0=sqrt(1.-N[2]*N[2]-4.*c_origin*N[2])/2./fabs(N[2]); |
526 | x0=-N[0]*one_over_2Nz; |
527 | y0=-N[1]*one_over_2Nz; |
528 | r0=sqrt(1.-N[2]*N[2]-4.*c_origin*N[2])*fabs(one_over_2Nz); |
529 | |
530 | // Phi value at "vertex" |
531 | phi=atan2(-x0,y0); |
532 | if(phi<0)phi+=M_TWO_PI6.28318530717958647692; |
533 | if(phi>=M_TWO_PI6.28318530717958647692)phi-=M_TWO_PI6.28318530717958647692; |
534 | |
535 | // Calculate the chisq |
536 | ChisqCircle(); |
537 | chisq_source = CIRCLE; |
538 | |
539 | return NOERROR; |
540 | } |
541 | |
542 | |
543 | //----------------- |
544 | // FitCircleRiemannCorrected |
545 | //----------------- |
546 | // Riemann Circle fit with correction for non-normal track incidence |
547 | // |
548 | jerror_t DHelicalFit::FitCircleRiemannCorrected(float rc){ |
549 | // Covariance matrices |
550 | DMatrix CRPhi(hits.size(),hits.size()); |
551 | DMatrix CR(hits.size(),hits.size()); |
552 | // Auxiliary matrices for correcting for non-normal track incidence to FDC |
553 | // The correction is |
554 | // CRPhi'= C*CRPhi*C+S*CR*S, where S(i,i)=R_i*kappa/2 |
555 | // and C(i,i)=sqrt(1-S(i,i)^2) |
556 | DMatrix C(hits.size(),hits.size()); |
557 | DMatrix S(hits.size(),hits.size()); |
558 | for (unsigned int i=0;i<hits.size();i++){ |
559 | S(i,i)=0.; |
560 | C(i,i)=1.; |
561 | if (rc>0){ |
562 | double rtemp=sqrt(hits[i]->x*hits[i]->x+hits[i]->y*hits[i]->y); |
563 | double stemp=rtemp/4./rc; |
564 | double ctemp=1.-stemp*stemp; |
565 | if (ctemp>0){ |
566 | S(i,i)=stemp; |
567 | C(i,i)=sqrt(ctemp); |
568 | } |
569 | } |
570 | } |
571 | |
572 | // Covariance matrices |
573 | if (CovRPhi_==NULL__null){ |
574 | CovRPhi_ = new DMatrix(hits.size(),hits.size()); |
575 | for (unsigned int i=0;i<hits.size();i++){ |
576 | double Phi=atan2(hits[i]->y,hits[i]->x); |
577 | double cosPhi=cos(Phi); |
578 | double sinPhi=sin(Phi); |
579 | double Phi_sinPhi_plus_cosPhi=Phi*sinPhi+cosPhi; |
580 | double Phi_cosPhi_minus_sinPhi=Phi*cosPhi-sinPhi; |
581 | CovRPhi_->operator()(i,i) |
582 | =Phi_cosPhi_minus_sinPhi*Phi_cosPhi_minus_sinPhi*hits[i]->covx |
583 | +Phi_sinPhi_plus_cosPhi*Phi_sinPhi_plus_cosPhi*hits[i]->covy |
584 | +2.*Phi_sinPhi_plus_cosPhi*Phi_cosPhi_minus_sinPhi*hits[i]->covxy; |
585 | } |
586 | } |
587 | if (CovR_==NULL__null){ |
588 | CovR_= new DMatrix(hits.size(),hits.size()); |
589 | for (unsigned int m=0;m<hits.size();m++){ |
590 | double Phi=atan2(hits[m]->y,hits[m]->x); |
591 | double cosPhi=cos(Phi); |
592 | double sinPhi=sin(Phi); |
593 | CovR_->operator()(m,m)=cosPhi*cosPhi*hits[m]->covx |
594 | +sinPhi*sinPhi*hits[m]->covy |
595 | +2.*sinPhi*cosPhi*hits[m]->covxy; |
596 | } |
597 | } |
598 | for (unsigned int i=0;i<hits.size();i++){ |
599 | for (unsigned int j=0;j<hits.size();j++){ |
600 | CR(i,j)=CovR_->operator()(i, j); |
601 | CRPhi(i,j)=CovRPhi_->operator()(i, j); |
602 | } |
603 | } |
604 | |
605 | // Correction for non-normal incidence of track on FDC |
606 | CRPhi=C*CRPhi*C+S*CR*S; |
607 | for (unsigned int i=0;i<hits.size();i++) |
608 | for (unsigned int j=0;j<hits.size();j++) |
609 | CovRPhi_->operator()(i,j)=CRPhi(i,j); |
610 | return FitCircleRiemann(); |
611 | } |
612 | |
613 | |
614 | //----------------- |
615 | // GetChargeRiemann |
616 | //----------------- |
617 | // Charge-finding routine with corrected CRPhi (see above) |
618 | // |
619 | jerror_t DHelicalFit::GetChargeRiemann(float rc_input){ |
620 | // Covariance matrices |
621 | DMatrix CRPhi(hits.size(),hits.size()); |
622 | DMatrix CR(hits.size(),hits.size()); |
623 | // Auxiliary matrices for correcting for non-normal track incidence to FDC |
624 | // The correction is |
625 | // CRPhi'= C*CRPhi*C+S*CR*S, where S(i,i)=R_i*kappa/2 |
626 | // and C(i,i)=sqrt(1-S(i,i)^2) |
627 | DMatrix C(hits.size(),hits.size()); |
628 | DMatrix S(hits.size(),hits.size()); |
629 | for (unsigned int i=0;i<hits.size();i++){ |
630 | double rtemp=sqrt(hits[i]->x*hits[i]->x+hits[i]->y*hits[i]->y); |
631 | double stemp=rtemp/(4.*rc_input); |
632 | double ctemp=1.-stemp*stemp; |
633 | if (ctemp>0){ |
634 | S(i,i)=stemp; |
635 | C(i,i)=sqrt(ctemp); |
636 | } |
637 | else{ |
638 | S(i,i)=0.; |
639 | C(i,i)=1; |
640 | } |
641 | } |
642 | |
643 | // Covariance matrices |
644 | if (CovRPhi_==NULL__null){ |
645 | CovRPhi_ = new DMatrix(hits.size(),hits.size()); |
646 | for (unsigned int i=0;i<hits.size();i++){ |
647 | double Phi=atan2(hits[i]->y,hits[i]->x); |
648 | double cosPhi=cos(Phi); |
649 | double sinPhi=sin(Phi); |
650 | double Phi_sinPhi_plus_cosPhi=Phi*sinPhi+cosPhi; |
651 | double Phi_cosPhi_minus_sinPhi=Phi*cosPhi-sinPhi; |
652 | CovRPhi_->operator()(i,i) |
653 | =Phi_cosPhi_minus_sinPhi*Phi_cosPhi_minus_sinPhi*hits[i]->covx |
654 | +Phi_sinPhi_plus_cosPhi*Phi_sinPhi_plus_cosPhi*hits[i]->covy |
655 | +2.*Phi_sinPhi_plus_cosPhi*Phi_cosPhi_minus_sinPhi*hits[i]->covxy; |
656 | } |
657 | } |
658 | if (CovR_==NULL__null){ |
659 | CovR_= new DMatrix(hits.size(),hits.size()); |
660 | for (unsigned int m=0;m<hits.size();m++){ |
661 | double Phi=atan2(hits[m]->y,hits[m]->x); |
662 | double cosPhi=cos(Phi); |
663 | double sinPhi=sin(Phi); |
664 | CovR_->operator()(m,m)=cosPhi*cosPhi*hits[m]->covx |
665 | +sinPhi*sinPhi*hits[m]->covy |
666 | +2.*sinPhi*cosPhi*hits[m]->covxy; |
667 | } |
668 | } |
669 | for (unsigned int i=0;i<hits.size();i++){ |
670 | for (unsigned int j=0;j<hits.size();j++){ |
671 | CR(i,j)=CovR_->operator()(i, j); |
672 | CRPhi(i,j)=CovRPhi_->operator()(i, j); |
673 | } |
674 | } |
675 | // Correction for non-normal incidence of track on FDC |
676 | CRPhi=C*CRPhi*C+S*CR*S; |
677 | for (unsigned int i=0;i<hits.size();i++) |
678 | for (unsigned int j=0;j<hits.size();j++) |
679 | CovRPhi_->operator()(i,j)=CRPhi(i,j); |
680 | return GetChargeRiemann(); |
681 | } |
682 | |
683 | //----------------- |
684 | // GetChargeRiemann |
685 | //----------------- |
686 | // Linear regression to find charge |
687 | // |
688 | jerror_t DHelicalFit::GetChargeRiemann(){ |
689 | // Covariance matrices |
690 | DMatrix CRPhi(hits.size(),hits.size()); |
691 | DMatrix CR(hits.size(),hits.size()); |
692 | if (CovRPhi_==NULL__null){ |
693 | CovRPhi_ = new DMatrix(hits.size(),hits.size()); |
694 | for (unsigned int i=0;i<hits.size();i++){ |
695 | double Phi=atan2(hits[i]->y,hits[i]->x); |
696 | double cosPhi=cos(Phi); |
697 | double sinPhi=sin(Phi); |
698 | double Phi_cosPhi_minus_sinPhi=Phi*cosPhi-sinPhi; |
699 | double Phi_sinPhi_plus_cosPhi=Phi*sinPhi+cosPhi; |
700 | CovRPhi_->operator()(i,i) |
701 | =Phi_cosPhi_minus_sinPhi*Phi_cosPhi_minus_sinPhi*hits[i]->covx |
702 | +Phi_sinPhi_plus_cosPhi*Phi_sinPhi_plus_cosPhi*hits[i]->covy |
703 | +2.*Phi_sinPhi_plus_cosPhi*Phi_cosPhi_minus_sinPhi*hits[i]->covxy; |
704 | } |
705 | } |
706 | if (CovR_==NULL__null){ |
707 | CovR_= new DMatrix(hits.size(),hits.size()); |
708 | for (unsigned int m=0;m<hits.size();m++){ |
709 | double Phi=atan2(hits[m]->y,hits[m]->x); |
710 | double cosPhi=cos(Phi); |
711 | double sinPhi=sin(Phi); |
712 | CovR_->operator()(m,m)=cosPhi*cosPhi*hits[m]->covx |
713 | +sinPhi*sinPhi*hits[m]->covy |
714 | +2.*sinPhi*cosPhi*hits[m]->covxy; |
715 | } |
716 | } |
717 | for (unsigned int i=0;i<hits.size();i++){ |
718 | for (unsigned int j=0;j<hits.size();j++){ |
719 | CR(i,j)=CovR_->operator()(i, j); |
720 | CRPhi(i,j)=CovRPhi_->operator()(i, j); |
721 | } |
722 | } |
723 | |
724 | double phi_old=atan2(hits[0]->y,hits[0]->x); |
725 | double sumv=0,sumy=0,sumx=0,sumxx=0,sumxy=0; |
726 | for (unsigned int k=0;k<hits.size();k++){ |
727 | DHFHit_t *hit=hits[k]; |
728 | double phi_z=atan2(hit->y,hit->x); |
729 | |
730 | // Check for problem regions near +pi and -pi |
731 | if (fabs(phi_z-phi_old)>M_PI3.14159265358979323846){ |
732 | if (phi_old<0) phi_z-=2.*M_PI3.14159265358979323846; |
733 | else phi_z+=2.*M_PI3.14159265358979323846; |
734 | } |
735 | double inv_var=(hit->x*hit->x+hit->y*hit->y) |
736 | /(CRPhi(k,k)+phi_z*phi_z*CR(k,k)); |
737 | sumv+=inv_var; |
738 | sumy+=phi_z*inv_var; |
739 | sumx+=hit->z*inv_var; |
740 | sumxx+=hit->z*hit->z*inv_var; |
741 | sumxy+=phi_z*hit->z*inv_var; |
742 | phi_old=phi_z; |
743 | } |
744 | double slope=(sumv*sumxy-sumy*sumx)/(sumv*sumxx-sumx*sumx); |
745 | |
746 | // Guess particle charge (+/-1); |
747 | q=+1.; |
748 | if (slope<0.) q= -1.; |
749 | |
750 | return NOERROR; |
751 | } |
752 | |
753 | //----------------- |
754 | // FitLineRiemann |
755 | //----------------- |
756 | jerror_t DHelicalFit::FitLineRiemann(){ |
757 | /// Riemann Line fit: linear regression of s on z to determine the tangent of |
758 | /// the dip angle and the z position of the closest approach to the beam line. |
759 | /// Computes intersection points along the helical path. |
760 | /// |
761 | // Get covariance matrix |
762 | DMatrix CR(hits.size(),hits.size()); |
763 | if (CovR_==NULL__null){ |
764 | CovR_= new DMatrix(hits.size(),hits.size()); |
765 | for (unsigned int m=0;m<hits.size();m++){ |
766 | double Phi=atan2(hits[m]->y,hits[m]->x); |
767 | double cosPhi=cos(Phi); |
768 | double sinPhi=sin(Phi); |
769 | CovR_->operator()(m,m)=cosPhi*cosPhi*hits[m]->covx |
770 | +sinPhi*sinPhi*hits[m]->covy |
771 | +2.*sinPhi*cosPhi*hits[m]->covxy; |
772 | } |
773 | } |
774 | for (unsigned int i=0;i<hits.size();i++) |
775 | for (unsigned int j=0;j<hits.size();j++) |
776 | CR(i,j)=CovR_->operator()(i, j); |
777 | |
778 | // Fill vector of intersection points |
779 | double denom= N[0]*N[0]+N[1]*N[1]; |
780 | vector<int>bad(hits.size()); |
781 | int numbad=0; |
782 | // projection vector |
783 | vector<DVector2>projections; |
784 | for (unsigned int m=0;m<hits.size();m++){ |
785 | double r2=hits[m]->x*hits[m]->x+hits[m]->y*hits[m]->y; |
786 | double numer=c_origin+r2*N[2]; |
787 | |
788 | if (r2==0){ |
789 | projections.push_back(DVector2(0.,0.)); |
790 | } |
791 | else{ |
792 | double ratio=numer/denom; |
793 | double x_int0=-N[0]*ratio; |
794 | double y_int0=-N[1]*ratio; |
795 | double temp=denom*r2-numer*numer; |
796 | if (temp<0){ // Skip point if the intersection gives nonsense |
797 | bad[m]=1; |
798 | numbad++; |
799 | projections.push_back(DVector2(x_int0,y_int0)); |
800 | continue; |
801 | } |
802 | temp=sqrt(temp)/denom; |
803 | |
804 | // Choose sign of square root based on proximity to actual measurements |
805 | double deltax=N[1]*temp; |
806 | double deltay=-N[0]*temp; |
807 | double x1=x_int0+deltax; |
808 | double y1=y_int0+deltay; |
809 | double x2=x_int0-deltax; |
810 | double y2=y_int0-deltay; |
811 | double diffx1=x1-hits[m]->x; |
812 | double diffy1=y1-hits[m]->y; |
813 | double diffx2=x2-hits[m]->x; |
814 | double diffy2=y2-hits[m]->y; |
815 | if (diffx1*diffx1+diffy1*diffy1 > diffx2*diffx2+diffy2*diffy2){ |
816 | //temphit->x=x2; |
817 | //temphit->y=y2; |
818 | projections.push_back(DVector2(x2,y2)); |
819 | } |
820 | else{ |
821 | //temphit->x=x1; |
822 | //temphit->y=y1; |
823 | projections.push_back(DVector2(x1,y1)); |
824 | } |
825 | } |
826 | //projections.push_back(temphit); |
827 | } |
828 | |
829 | // All arc lengths are measured relative to some reference plane with a hit. |
830 | // Don't use a "bad" hit for the reference... |
831 | unsigned int start=0; |
832 | for (unsigned int i=0;i<bad.size();i++){ |
833 | if (!bad[i]){ |
834 | start=i; |
835 | break; |
836 | } |
837 | } |
838 | |
839 | // Linear regression to find z0, tanl |
840 | unsigned int n=projections.size(); |
841 | double sumv=0.,sumx=0.,sumy=0.,sumxx=0.,sumxy=0.; |
842 | double sperp=0.,sperp_old=0., ratio=0, Delta; |
843 | double z_last=0.,z=0.; |
844 | DVector2 old_proj=projections[start]; |
845 | double two_r0=2.*r0; |
846 | for (unsigned int k=start;k<n;k++){ |
847 | if (!bad[k]){ |
848 | sperp_old=sperp; |
849 | z_last=z; |
850 | ratio=(projections[k]-old_proj).Mod()/(two_r0); |
851 | // Make sure the argument for the arcsin does not go out of range... |
852 | sperp=sperp_old+(ratio>1? two_r0*M_PI_21.57079632679489661923 : two_r0*asin(ratio)); |
853 | z=hits[k]->z; |
854 | |
855 | // Assume errors in s dominated by errors in R |
856 | double weight=1./CR(k,k); |
857 | sumv+=weight; |
858 | sumy+=sperp*weight; |
859 | sumx+=z*weight; |
860 | sumxx+=z*z*weight; |
861 | sumxy+=sperp*z*weight; |
862 | |
863 | // Store the current x and y projection values |
864 | old_proj=projections[k]; |
865 | } |
866 | } |
867 | Delta=sumv*sumxx-sumx*sumx; |
868 | double tanl_denom=sumv*sumxy-sumy*sumx; |
869 | if (fabs(Delta)<EPS1e-8 || fabs(tanl_denom)<EPS1e-8) return VALUE_OUT_OF_RANGE; |
870 | |
871 | // Track parameters tan(lambda) and z-vertex |
872 | tanl=-Delta/tanl_denom; |
873 | //z_vertex=(sumxx*sumy-sumx*sumxy)/Delta; |
874 | sperp-=sperp_old; |
875 | z_vertex=z_last-sperp*tanl; |
876 | |
877 | /* |
878 | if (z_vertex<Z_MIN){ |
879 | z_vertex=Z_MIN; |
880 | tanl=(z_last-Z_MIN)/sperp; |
881 | } |
882 | else if (z_vertex>Z_MAX){ |
883 | z_vertex=Z_MAX; |
884 | tanl=(z_last-Z_MAX)/sperp; |
885 | } |
886 | */ |
887 | theta=M_PI_21.57079632679489661923-atan(tanl); |
888 | |
889 | return NOERROR; |
890 | } |
891 | |
892 | //----------------- |
893 | // FitCircleStraightTrack |
894 | //----------------- |
895 | jerror_t DHelicalFit::FitCircleStraightTrack(void) |
896 | { |
897 | /// This is a last resort! |
898 | /// The circle fit can fail for high momentum tracks that are nearly |
899 | /// straight tracks. In these cases (when pt~2GeV or more) there |
900 | /// is not enough position resolution with wire positions only |
901 | /// to measure the curvature of the track. |
902 | /// We can though, fit the X/Y points with a straight line in order |
903 | /// to get phi and project out the stereo angles. |
904 | /// |
905 | /// For the radius of the circle (i.e. p_trans) we loop over momenta |
906 | /// from 1 to 9GeV and just use the one with the smallest chisq. |
907 | |
908 | // Average the x's and y's of the individual hits. We should be |
909 | // able to do this via linear regression, but I can't get it to |
910 | // work right now and I'm under the gun to get ready for a review |
911 | // so I take the easy way. Note that we don't average phi's since |
912 | // that will cause problems at the 0/2pi boundary. |
913 | double X=0.0, Y=0.0; |
914 | DHFHit_t *a = NULL__null; |
915 | for(unsigned int i=0;i<hits.size();i++){ |
916 | a = hits[i]; |
917 | double r = sqrt(pow((double)a->x,2.0) + pow((double)a->y, 2.0)); |
918 | // weight by r to give outer points more influence. Note that |
919 | // we really are really weighting by r^2 since x and y already |
920 | // have a magnitude component. |
921 | X += a->x*r; |
922 | Y += a->y*r; |
923 | } |
924 | phi = atan2(Y,X); |
925 | if(phi<0)phi+=M_TWO_PI6.28318530717958647692; |
926 | if(phi>=M_TWO_PI6.28318530717958647692)phi-=M_TWO_PI6.28318530717958647692; |
927 | |
928 | // Search the chi2 space for values for p_trans, x0, ... |
929 | SearchPtrans(9.0, 0.5); |
930 | |
931 | #if 0 |
932 | // We do a simple linear regression here to find phi. This is |
933 | // simplified by the intercept being zero (i.e. the track |
934 | // passes through the beamline). |
935 | float Sxx=0.0, Syy=0.0, Sxy=0.0; |
936 | chisq_source = NOFIT; // in case we return early |
937 | |
938 | // Loop over hits to calculate Sxx, Syy, and Sxy |
939 | DHFHit_t *a = NULL__null; |
940 | for(unsigned int i=0;i<hits.size();i++){ |
941 | a = hits[i]; |
942 | float x=a->x; |
943 | float y=a->y; |
944 | Sxx += x*x; |
945 | Syy += y*y; |
946 | Sxy += x*y; |
947 | } |
948 | double A = 2.0*Sxy; |
949 | double B = Sxx - Syy; |
950 | phi = B>A ? atan2(A,B)/2.0 : 1.0/atan2(B,A)/2.0; |
951 | if(phi<0)phi+=M_TWO_PI6.28318530717958647692; |
952 | if(phi>=M_TWO_PI6.28318530717958647692)phi-=M_TWO_PI6.28318530717958647692; |
953 | #endif |
954 | |
955 | return NOERROR; |
956 | } |
957 | |
958 | //----------------- |
959 | // SearchPtrans |
960 | //----------------- |
961 | void DHelicalFit::SearchPtrans(double ptrans_max, double ptrans_step) |
962 | { |
963 | /// Search the chi2 space as a function of the transverse momentum |
964 | /// for the minimal chi2. The values corresponding to the minmal |
965 | /// chi2 are left in chisq, x0, y0, r0, q, and p_trans. |
966 | /// |
967 | /// This does NOT optimize the value of p_trans. It simply |
968 | /// does a straight forward chisq calculation on a grid |
969 | /// and keeps the best one. It is intended for use in finding |
970 | /// a reasonable value for p_trans for straight tracks that |
971 | /// are contained to less than p_trans_max which is presumably |
972 | /// chosen based on a known θ. |
973 | |
974 | // Loop over several values for p_trans and calculate the |
975 | // chisq for each. |
976 | double min_chisq=1.0E6; |
977 | double min_x0=0.0, min_y0=0.0, min_r0=0.0; |
978 | for(double my_p_trans=ptrans_step; my_p_trans<=ptrans_max; my_p_trans+=ptrans_step){ |
979 | |
980 | r0 = my_p_trans/(0.003*2.0); |
981 | double alpha, my_chisq; |
982 | |
983 | // q = +1 |
984 | alpha = phi + M_PI_21.57079632679489661923; |
985 | x0 = r0*cos(alpha); |
986 | y0 = r0*sin(alpha); |
987 | my_chisq = ChisqCircle(); |
988 | if(my_chisq<min_chisq){ |
989 | min_chisq=my_chisq; |
990 | min_x0 = x0; |
991 | min_y0 = y0; |
992 | min_r0 = r0; |
993 | q = +1.0; |
994 | p_trans = my_p_trans; |
995 | } |
996 | |
997 | // q = -1 |
998 | alpha = phi - M_PI_21.57079632679489661923; |
999 | x0 = r0*cos(alpha); |
1000 | y0 = r0*sin(alpha); |
1001 | my_chisq = ChisqCircle(); |
1002 | if(my_chisq<min_chisq){ |
1003 | min_chisq=my_chisq; |
1004 | min_x0 = x0; |
1005 | min_y0 = y0; |
1006 | min_r0 = r0; |
1007 | q = -1.0; |
1008 | p_trans = my_p_trans; |
1009 | } |
1010 | } |
1011 | |
1012 | // Copy params from minimum chisq |
1013 | x0 = min_x0; |
1014 | y0 = min_y0; |
1015 | r0 = min_r0; |
1016 | } |
1017 | |
1018 | //----------------- |
1019 | // QuickPtrans |
1020 | //----------------- |
1021 | void DHelicalFit::QuickPtrans(void) |
1022 | { |
1023 | /// Quickly calculate a value of p_trans by looking |
1024 | /// for the hit furthest out and the hit closest |
1025 | /// to half that distance. Those 2 hits along with the |
1026 | /// origin are used to define a circle from which |
1027 | /// p_trans is calculated. |
1028 | |
1029 | // Find hit with largest R |
1030 | double R2max = 0.0; |
1031 | DHFHit_t *hit_max = NULL__null; |
1032 | for(unsigned int i=0; i<hits.size(); i++){ |
1033 | // use cross product to decide if hits is to left or right |
1034 | double x = hits[i]->x; |
1035 | double y = hits[i]->y; |
1036 | double R2 = x*x + y*y; |
1037 | if(R2>R2max){ |
1038 | R2max = R2; |
1039 | hit_max = hits[i]; |
1040 | } |
1041 | } |
1042 | |
1043 | // Bullet proof |
1044 | if(!hit_max)return; |
1045 | |
1046 | // Find hit closest to half-way out |
1047 | double Rmid = 0.0; |
1048 | double Rmax = sqrt(R2max); |
1049 | DHFHit_t *hit_mid = NULL__null; |
1050 | for(unsigned int i=0; i<hits.size(); i++){ |
1051 | // use cross product to decide if hits is to left or right |
1052 | double x = hits[i]->x; |
1053 | double y = hits[i]->y; |
1054 | double R = sqrt(x*x + y*y); |
1055 | if(fabs(R-Rmax/2.0) < fabs(Rmid-Rmax/2.0)){ |
1056 | Rmid = R; |
1057 | hit_mid = hits[i]; |
1058 | } |
1059 | } |
1060 | |
1061 | // Bullet proof |
1062 | if(!hit_mid)return; |
1063 | |
1064 | // Calculate p_trans from 3 points |
1065 | double x1 = hit_mid->x; |
1066 | double y1 = hit_mid->y; |
1067 | double x2 = hit_max->x; |
1068 | double y2 = hit_max->y; |
1069 | double r2 = sqrt(x2*x2 + y2*y2); |
1070 | double cos_phi = x2/r2; |
1071 | double sin_phi = y2/r2; |
1072 | double u1 = x1*cos_phi + y1*sin_phi; |
1073 | double v1 = -x1*sin_phi + y1*cos_phi; |
1074 | double u2 = x2*cos_phi + y2*sin_phi; |
1075 | double u0 = u2/2.0; |
1076 | double v0 = (u1*u1 + v1*v1 - u1*u2)/(2.0*v1); |
1077 | |
1078 | x0 = u0*cos_phi - v0*sin_phi; |
1079 | y0 = u0*sin_phi + v0*cos_phi; |
1080 | r0 = sqrt(x0*x0 + y0*y0); |
1081 | |
1082 | double B=-2.0; |
1083 | p_trans = fabs(0.003*B*r0); |
1084 | q = v1>0.0 ? -1.0:+1.0; |
1085 | } |
1086 | |
1087 | //----------------- |
1088 | // GuessChargeFromCircleFit |
1089 | //----------------- |
1090 | jerror_t DHelicalFit::GuessChargeFromCircleFit(void) |
1091 | { |
1092 | /// Adjust the sign of the charge (magnitude will stay 1) based on |
1093 | /// whether the hits tend to be to the right or to the left of the |
1094 | /// line going from the origin through the center of the circle. |
1095 | /// If the sign is flipped, the phi angle will also be shifted by |
1096 | /// +/- pi since the majority of hits are assumed to come from |
1097 | /// outgoing tracks. |
1098 | /// |
1099 | /// This is just a guess since tracks can bend all the way |
1100 | /// around and have hits that look exactly like tracks from an |
1101 | /// outgoing particle of opposite charge. The final charge should |
1102 | /// come from the sign of the dphi/dz slope. |
1103 | |
1104 | // Simply count the number of hits whose phi angle relative to the |
1105 | // phi of the center of the circle are greater than pi. |
1106 | unsigned int N=0; |
1107 | for(unsigned int i=0; i<hits.size(); i++){ |
1108 | // use cross product to decide if hits is to left or right |
1109 | double x = hits[i]->x; |
1110 | double y = hits[i]->y; |
1111 | if((x*y0 - y*x0) < 0.0)N++; |
1112 | } |
1113 | |
1114 | // Check if more hits are negative and make sign negative if so. |
1115 | if(N>hits.size()/2.0){ |
1116 | q = -1.0; |
1117 | phi += M_PI3.14159265358979323846; |
1118 | if(phi>M_TWO_PI6.28318530717958647692)phi-=M_TWO_PI6.28318530717958647692; |
1119 | } |
1120 | |
1121 | return NOERROR; |
1122 | } |
1123 | |
1124 | //----------------- |
1125 | // FitTrack |
1126 | //----------------- |
1127 | jerror_t DHelicalFit::FitTrack(void) |
1128 | { |
1129 | /// Find theta, sign of electric charge, total momentum and |
1130 | /// vertex z position. |
1131 | |
1132 | // Points must be in order of increasing Z |
1133 | sort(hits.begin(), hits.end(), DHFHitLessThanZ_C); |
1134 | |
1135 | // Fit to circle to get circle's center |
1136 | FitCircle(); |
1137 | |
1138 | // The thing that is really needed is dphi/dz (where phi is the angle |
1139 | // of the point as measured from the center of the circle, not the beam |
1140 | // line). The relation between phi and z is linear so we use linear |
1141 | // regression to find the slope (dphi/dz). The one complication is |
1142 | // that phi is periodic so the value obtained via the x and y of a |
1143 | // specific point can be off by an integral number of 2pis. Handle this |
1144 | // by assuming the first point is measured before a full rotation |
1145 | // has occurred. Subsequent points should always be within 2pi of |
1146 | // the previous point so we just need to calculate the relative phi |
1147 | // between succesive points and keep a running sum. We do this in |
1148 | // the first loop were we find the mean z and phi values. The regression |
1149 | // formulae are calculated in the second loop. |
1150 | |
1151 | // Calculate phi about circle center for each hit |
1152 | Fill_phi_circle(hits, x0, y0); |
1153 | |
1154 | // Linear regression for z/phi relation |
1155 | float Sxx=0.0, Syy=0.0, Sxy=0.0; |
1156 | for(unsigned int i=0;i<hits.size();i++){ |
1157 | DHFHit_t *a = hits[i]; |
1158 | float deltaZ = a->z - z_mean; |
1159 | float deltaPhi = a->phi_circle - phi_mean; |
1160 | Syy += deltaZ*deltaZ; |
1161 | Sxx += deltaPhi*deltaPhi; |
1162 | Sxy += deltaZ*deltaPhi; |
1163 | //cout<<__FILE__<<":"<<__LINE__<<" deltaZ="<<deltaZ<<" deltaPhi="<<deltaPhi<<" Sxy(i)="<<deltaZ*deltaPhi<<endl; |
1164 | } |
1165 | float dzdphi = Syy/Sxy; |
Value stored to 'dzdphi' during its initialization is never read | |
1166 | dzdphi = Syy/Sxy; |
1167 | z_vertex = z_mean - phi_mean*dzdphi; |
1168 | //cout<<__FILE__<<":"<<__LINE__<<" z_mean="<<z_mean<<" phi_mean="<<phi_mean<<" dphidz="<<dphidz<<" Sxy="<<Sxy<<" Syy="<<Syy<<" z_vertex="<<z_vertex<<endl; |
1169 | |
1170 | // Fill in the rest of the parameters |
1171 | return FillTrackParams(); |
1172 | } |
1173 | |
1174 | //----------------- |
1175 | // FitTrackRiemann |
1176 | //----------------- |
1177 | jerror_t DHelicalFit::FitTrackRiemann(float rc_input){ |
1178 | jerror_t error=NOERROR; |
1179 | |
1180 | if (CovR_!=NULL__null) delete CovR_; |
1181 | if (CovRPhi_!=NULL__null) delete CovRPhi_; |
1182 | CovR_=NULL__null; |
1183 | CovRPhi_=NULL__null; |
1184 | |
1185 | error=FitCircleRiemannCorrected(rc_input); |
1186 | if (error!=NOERROR) return error; |
1187 | error=FitLineRiemann(); |
1188 | GetChargeRiemann(); |
1189 | |
1190 | // Shift phi by pi if the charge is negative |
1191 | if (q<0) phi+=M_PI3.14159265358979323846; |
1192 | |
1193 | return error; |
1194 | } |
1195 | |
1196 | |
1197 | jerror_t DHelicalFit::FitCircleAndLineRiemann(float rc_input){ |
1198 | jerror_t error=NOERROR; |
1199 | |
1200 | if (CovR_!=NULL__null) delete CovR_; |
1201 | if (CovRPhi_!=NULL__null) delete CovRPhi_; |
1202 | CovR_=NULL__null; |
1203 | CovRPhi_=NULL__null; |
1204 | |
1205 | error=FitCircleRiemannCorrected(rc_input); |
1206 | if (error!=NOERROR) return error; |
1207 | error=FitLineRiemann(); |
1208 | |
1209 | return error; |
1210 | } |
1211 | |
1212 | |
1213 | |
1214 | |
1215 | //----------------- |
1216 | // FitTrack_FixedZvertex |
1217 | //----------------- |
1218 | jerror_t DHelicalFit::FitTrack_FixedZvertex(float z_vertex) |
1219 | { |
1220 | /// Fit the points, but hold the z_vertex fixed at the specified value. |
1221 | /// |
1222 | /// This just calls FitCircle and FitLine_FixedZvertex to find all parameters |
1223 | /// of the track |
1224 | |
1225 | // Fit to circle to get circle's center |
1226 | FitCircle(); |
1227 | |
1228 | // Fit to line in phi-z plane and return error |
1229 | return FitLine_FixedZvertex(z_vertex); |
1230 | } |
1231 | |
1232 | //----------------- |
1233 | // FitLine_FixedZvertex |
1234 | //----------------- |
1235 | jerror_t DHelicalFit::FitLine_FixedZvertex(float z_vertex) |
1236 | { |
1237 | /// Fit the points, but hold the z_vertex fixed at the specified value. |
1238 | /// |
1239 | /// This assumes FitCircle has already been called and the values |
1240 | /// in x0 and y0 are valid. |
1241 | /// |
1242 | /// This just fits the phi-z angle by minimizing the chi-squared |
1243 | /// using a linear regression technique. As it turns out, the |
1244 | /// chi-squared weights points by their distances squared which |
1245 | /// leads to a quadratic equation for cot(theta) (where theta is |
1246 | /// the angle in the phi-z plane). To pick the right solution of |
1247 | /// the quadratic equation, we pick the one closest to the linearly |
1248 | /// weighted angle. One could argue that we should just use the |
1249 | /// linear weighting here rather than the square weighting. The |
1250 | /// choice depends though on how likely the "out-lier" points |
1251 | /// are to really be from this track. If they are likely, to |
1252 | /// be, then we would want to leverage their "longer" lever arms |
1253 | /// with the squared weighting. If they are more likely to be bad |
1254 | /// points, then we would want to minimize their influence with |
1255 | /// a linear (or maybe even root) weighting. It is expected than |
1256 | /// for our use, the points are all likely valid so we use the |
1257 | /// square weighting. |
1258 | |
1259 | // Points must be in order of increasing Z |
1260 | sort(hits.begin(), hits.end(), DHFHitLessThanZ_C); |
1261 | |
1262 | // Fit is being done for a fixed Z-vertex |
1263 | this->z_vertex = z_vertex; |
1264 | |
1265 | // Calculate phi about circle center for each hit |
1266 | Fill_phi_circle(hits, x0, y0); |
1267 | |
1268 | // Do linear regression on phi-Z |
1269 | float Sx=0, Sy=0; |
1270 | float Sxx=0, Syy=0, Sxy = 0; |
1271 | float r0 = sqrt(x0*x0 + y0*y0); |
1272 | for(unsigned int i=0; i<hits.size(); i++){ |
1273 | DHFHit_t *a = hits[i]; |
1274 | |
1275 | float dz = a->z - z_vertex; |
1276 | float dphi = a->phi_circle; |
1277 | Sx += dz; |
1278 | Sy += r0*dphi; |
1279 | Syy += r0*dphi*r0*dphi; |
1280 | Sxx += dz*dz; |
1281 | Sxy += r0*dphi*dz; |
1282 | } |
1283 | |
1284 | float k = (Syy-Sxx)/(2.0*Sxy); |
1285 | float s = sqrt(k*k + 1.0); |
1286 | float cot_theta1 = -k+s; |
1287 | float cot_theta2 = -k-s; |
1288 | float cot_theta_lin = Sx/Sy; |
1289 | float cot_theta; |
1290 | if(fabs(cot_theta_lin-cot_theta1) < fabs(cot_theta_lin-cot_theta2)){ |
1291 | cot_theta = cot_theta1; |
1292 | }else{ |
1293 | cot_theta = cot_theta2; |
1294 | } |
1295 | |
1296 | dzdphi = r0*cot_theta; |
1297 | |
1298 | // Fill in the rest of the paramters |
1299 | return FillTrackParams(); |
1300 | } |
1301 | |
1302 | //------------------------------------------------------------------ |
1303 | // Fill_phi_circle |
1304 | //------------------------------------------------------------------ |
1305 | jerror_t DHelicalFit::Fill_phi_circle(vector<DHFHit_t*> hits, float x0, float y0) |
1306 | { |
1307 | float x_last = -x0; |
1308 | float y_last = -y0; |
1309 | float phi_last = 0.0; |
1310 | z_mean = phi_mean = 0.0; |
1311 | for(unsigned int i=0; i<hits.size(); i++){ |
1312 | DHFHit_t *a = hits[i]; |
1313 | |
1314 | float dx = a->x - x0; |
1315 | float dy = a->y - y0; |
1316 | float dphi = atan2f(dx*y_last - dy*x_last, dx*x_last + dy*y_last)atan2((double)dx*y_last - dy*x_last,(double)dx*x_last + dy*y_last ); |
1317 | float my_phi = phi_last +dphi; |
1318 | a->phi_circle = my_phi; |
1319 | |
1320 | z_mean += a->z; |
1321 | phi_mean += my_phi; |
1322 | |
1323 | x_last = dx; |
1324 | y_last = dy; |
1325 | phi_last = my_phi; |
1326 | } |
1327 | z_mean /= (float)hits.size(); |
1328 | phi_mean /= (float)hits.size(); |
1329 | |
1330 | return NOERROR; |
1331 | } |
1332 | |
1333 | //------------------------------------------------------------------ |
1334 | // FillTrackParams |
1335 | //------------------------------------------------------------------ |
1336 | jerror_t DHelicalFit::FillTrackParams(void) |
1337 | { |
1338 | /// Fill in and tweak some parameters like q, phi, theta using |
1339 | /// other values already set in the class. This is used by |
1340 | /// both FitTrack() and FitTrack_FixedZvertex(). |
1341 | |
1342 | float r0 = sqrt(x0*x0 + y0*y0); |
1343 | theta = atan(r0/fabs(dzdphi)); |
1344 | |
1345 | // The sign of the electric charge will be opposite that |
1346 | // of dphi/dz. Also, the value of phi will be PI out of phase |
1347 | if(dzdphi<0.0){ |
1348 | phi += M_PI3.14159265358979323846; |
1349 | if(phi<0)phi+=M_TWO_PI6.28318530717958647692; |
1350 | if(phi>=M_TWO_PI6.28318530717958647692)phi-=M_TWO_PI6.28318530717958647692; |
1351 | }else{ |
1352 | q = -q; |
1353 | } |
1354 | |
1355 | // Re-calculate chi-sq (needs to be done) |
1356 | chisq_source = TRACK; |
1357 | |
1358 | // Up to now, the fit has assumed a forward going particle. In |
1359 | // other words, if the particle is going backwards, the helix does |
1360 | // actually still go through the points, but only if extended |
1361 | // backward in time. We use the average z of the hits compared |
1362 | // to the reconstructed vertex to determine if it was back-scattered. |
1363 | if(z_mean<z_vertex){ |
1364 | // back scattered particle |
1365 | theta = M_PI3.14159265358979323846 - theta; |
1366 | phi += M_PI3.14159265358979323846; |
1367 | if(phi<0)phi+=M_TWO_PI6.28318530717958647692; |
1368 | if(phi>=M_TWO_PI6.28318530717958647692)phi-=M_TWO_PI6.28318530717958647692; |
1369 | q = -q; |
1370 | } |
1371 | |
1372 | // There is a problem with theta for data generated with a uniform |
1373 | // field. The following factor is emprical until I can figure out |
1374 | // what the source of the descrepancy is. |
1375 | //theta += 0.1*pow((double)sin(theta),2.0); |
1376 | |
1377 | p = fabs(p_trans/sin(theta)); |
1378 | |
1379 | return NOERROR; |
1380 | } |
1381 | |
1382 | //------------------------------------------------------------------ |
1383 | |
1384 | //------------------------------------------------------------------ |
1385 | jerror_t DHelicalFit::Print(void) const |
1386 | { |
1387 | cout<<"-- DHelicalFit Params ---------------"<<endl; |
1388 | cout<<" x0 = "<<x0<<endl; |
1389 | cout<<" y0 = "<<y0<<endl; |
1390 | cout<<" q = "<<q<<endl; |
1391 | cout<<" p = "<<p<<endl; |
1392 | cout<<" p_trans = "<<p_trans<<endl; |
1393 | cout<<" phi = "<<phi<<endl; |
1394 | cout<<" theta = "<<theta<<endl; |
1395 | cout<<" z_vertex = "<<z_vertex<<endl; |
1396 | cout<<" chisq = "<<chisq<<endl; |
1397 | cout<<"chisq_source = "; |
1398 | switch(chisq_source){ |
1399 | case NOFIT: cout<<"NOFIT"; break; |
1400 | case CIRCLE: cout<<"CIRCLE"; break; |
1401 | case TRACK: cout<<"TRACK"; break; |
1402 | } |
1403 | cout<<endl; |
1404 | cout<<" Bz(avg) = "<<Bz_avg<<endl; |
1405 | |
1406 | return NOERROR; |
1407 | } |
1408 | |
1409 | |
1410 | //------------------------------------------------------------------ |
1411 | // Dump |
1412 | //------------------------------------------------------------------ |
1413 | jerror_t DHelicalFit::Dump(void) const |
1414 | { |
1415 | Print(); |
1416 | |
1417 | for(unsigned int i=0;i<hits.size();i++){ |
1418 | DHFHit_t *v = hits[i]; |
1419 | cout<<" x="<<v->x<<" y="<<v->y<<" z="<<v->z; |
1420 | cout<<" phi_circle="<<v->phi_circle<<" chisq="<<v->chisq<<endl; |
1421 | } |
1422 | |
1423 | return NOERROR; |
1424 | } |