StRoot: Stl3Util/ftf/FtfBaseTrack.cxx Source File

00001 //:>----------------------------------------------------------------- 00002 //: FILE: FtfBaseTrack.cxx 00003 //: HISTORY: 00004 //: 28oct1996 version 1.00 00005 //: 11aug1999 ppy primary flag in FtfPara replace with vertexConstrainedFit 00006 //: 11aug1999 ppy primary flag in track filled now 00007 //: 3sep1999 ppy fitLine, dpsi cannot be greater than 1. Check introduced 00008 //: 5oct1999 ppy fitLine, bug corrected 00009 //: 6oct1999 ppy Root switch added 00010 //: 9mar2000 ppy extrapolation methods added 00011 //: 9mar2000 ppy lots of changes to use void pointers 00012 //: 28mar2000 ppy closestApproach split in two methods 00013 //: 19jul2000 ppy fitHelix, for secondary tracks r0,phi0,z0 00014 //: is given at a point on the helix close to 00015 //: the inner most point. Previously it was given 00016 //: exactly at the inner most point. Residuals 00017 //: are slightly better with this change 00018 //:<------------------------------------------------------------------ 00019 //:>------------------------------------------------------------------ 00020 //: CLASS: FtfBaseTrack 00021 //: DESCRIPTION: Basic Description of a track P 00022 //: AUTHOR: ppy - Pablo Yepes, yepes@rice.edu 00023 //:>------------------------------------------------------------------ 00024 #include "Stl3Util/ftf/FtfBaseTrack.h" 00025 00026 00027 00028 void ftfMatrixDiagonal ( double *h, double &h11, double &h22, double &h33 ) ; 00029 00030 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00031 // Track Initialization 00032 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00033 FtfBaseTrack::FtfBaseTrack ( ){ 00034 firstHit = 0 ; 00035 lastHit = 0 ; 00036 } 00037 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00038 // Fit a circle 00039 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00040 int FtfBaseTrack::fitHelix ( ) 00041 { 00042 if ( fitCircle ( ) ){ 00043 ftfLog ( " Problem in Fit_Circle " ) ; 00044 return 1 ; 00045 } 00046 // 00047 // Fit line in s-z plane now 00048 // 00049 if ( fitLine ( )) { 00050 ftfLog ( " Problem fitting a line " ) ; 00051 return 1 ; 00052 } 00053 return 0 ; 00054 } 00055 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00056 // End of Fit Helix 00057 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00058 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00059 // 00060 // Fits circle parameters using algorithm 00061 // described by ChErnov and Oskov in Computer Physics 00062 // Communications. 00063 // 00064 // Written in FORTRAN by Jawluen Tang, Physics department , UT-Austin 00065 // Moved to C by Pablo Yepes 00066 //--------------------------------------------------------------- 00067 int FtfBaseTrack::fitCircle ( ) 00068 { 00069 double wsum = 0.0 ; 00070 double xav = 0.0 ; 00071 double yav = 0.0 ; 00072 // 00073 // Loop over hits calculating average 00074 // 00075 for ( startLoop() ; done() ; nextHit() ) { 00076 00077 FtfBaseHit* cHit = (FtfBaseHit *)currentHit ; 00078 cHit->wxy = 1.0F/ (double)(cHit->dx*cHit->dx + cHit->dy*cHit->dy) ; 00079 wsum += cHit->wxy ; 00080 xav += cHit->wxy * cHit->x ; 00081 yav += cHit->wxy * cHit->y ; 00082 } 00083 if ( getPara()->vertexConstrainedFit ) { 00084 wsum += getPara()->xyWeightVertex ; 00085 xav += getPara()->xVertex ; 00086 yav += getPara()->yVertex ; 00087 } 00088 xav = xav / wsum ; 00089 yav = yav / wsum ; 00090 // 00091 // CALCULATE <X**2>, <XY>, AND <Y**2> WITH <X> = 0, & <Y> = 0 00092 // 00093 double xxav = 0.0 ; 00094 double xyav = 0.0 ; 00095 double yyav = 0.0 ; 00096 double xi, yi ; 00097 00098 for ( startLoop() ; done() ; nextHit() ) { 00099 FtfBaseHit* cHit = (FtfBaseHit *)currentHit ; 00100 xi = cHit->x - xav ; 00101 yi = cHit->y - yav ; 00102 xxav += xi * xi * cHit->wxy ; 00103 xyav += xi * yi * cHit->wxy ; 00104 yyav += yi * yi * cHit->wxy ; 00105 } 00106 if ( getPara()->vertexConstrainedFit ) { 00107 xi = getPara()->xVertex - xav ; 00108 yi = getPara()->yVertex - yav ; 00109 xxav += xi * xi * getPara()->xyWeightVertex ; 00110 xyav += xi * yi * getPara()->xyWeightVertex ; 00111 yyav += yi * yi * getPara()->xyWeightVertex ; 00112 } 00113 xxav = xxav / wsum ; 00114 xyav = xyav / wsum ; 00115 yyav = yyav / wsum ; 00116 // 00117 //--> ROTATE COORDINATES SO THAT <XY> = 0 00118 // 00119 //--> SIGN(C**2 - S**2) = SIGN(XXAV - YYAV) > 00120 //--> & > ==> NEW : (XXAV-YYAV) > 0 00121 //--> SIGN(S) = SIGN(XYAV) > 00122 00123 double a = fabs( xxav - yyav ) ; 00124 double b = 4.0 * xyav * xyav ; 00125 00126 double asqpb = a * a + b ; 00127 double rasqpb = sqrt ( asqpb) ; 00128 00129 double splus = 1.0 + a / rasqpb ; 00130 double sminus = b / (asqpb * splus) ; 00131 00132 splus = sqrt (0.5 * splus ) ; 00133 sminus = sqrt (0.5 * sminus) ; 00134 // 00135 //-> FIRST REQUIRE : SIGN(C**2 - S**2) = SIGN(XXAV - YYAV) 00136 // 00137 double sinrot, cosrot ; 00138 if ( xxav <= yyav ) { 00139 cosrot = sminus ; 00140 sinrot = splus ; 00141 } 00142 else { 00143 cosrot = splus ; 00144 sinrot = sminus ; 00145 } 00146 // 00147 //-> REQUIRE : SIGN(S) = SIGN(XYAV) * SIGN(C) (ASSUMING SIGN(C) > 0) 00148 // 00149 if ( xyav < 0.0 ) sinrot = - sinrot ; 00150 // 00151 //--> WE NOW HAVE THE SMALLEST ANGLE THAT GUARANTEES <X**2> > <Y**2> 00152 //--> TO GET THE SIGN OF THE CHARGE RIGHT, THE NEW X-AXIS MUST POINT 00153 //--> OUTWARD FROM THE ORGIN. WE ARE FREE TO CHANGE SIGNS OF BOTH 00154 //--> COSROT AND SINROT SIMULTANEOUSLY TO ACCOMPLISH THIS. 00155 // 00156 //--> CHOOSE SIGN OF C WISELY TO BE ABLE TO GET THE SIGN OF THE CHARGE 00157 // 00158 if ( cosrot*xav+sinrot*yav < 0.0 ) { 00159 cosrot = -cosrot ; 00160 sinrot = -sinrot ; 00161 } 00162 // 00163 //-> NOW GET <R**2> AND RSCALE= SQRT(<R**2>) 00164 // 00165 double rrav = xxav + yyav ; 00166 double rscale = ::sqrt(rrav) ; 00167 00168 xxav = 0.0 ; 00169 yyav = 0.0 ; 00170 xyav = 0.0 ; 00171 double xrrav = 0.0 ; 00172 double yrrav = 0.0 ; 00173 double rrrrav = 0.0 ; 00174 00175 double xixi, yiyi, riri, wiriri, xold, yold ; 00176 for ( startLoop() ; done() ; nextHit() ) { 00177 FtfBaseHit* cHit = (FtfBaseHit *)currentHit ; 00178 xold = cHit->x - xav ; 00179 yold = cHit->y - yav ; 00180 // 00181 //--> ROTATE SO THAT <XY> = 0 & DIVIDE BY RSCALE SO THAT <R**2> = 1 00182 // 00183 xi = ( cosrot * xold + sinrot * yold ) / rscale ; 00184 yi = ( -sinrot * xold + cosrot * yold ) / rscale ; 00185 00186 xixi = xi * xi ; 00187 yiyi = yi * yi ; 00188 riri = xixi + yiyi ; 00189 wiriri = cHit->wxy * riri ; 00190 00191 xyav += cHit->wxy * xi * yi ; 00192 xxav += cHit->wxy * xixi ; 00193 yyav += cHit->wxy * yiyi ; 00194 00195 xrrav += wiriri * xi ; 00196 yrrav += wiriri * yi ; 00197 rrrrav += wiriri * riri ; 00198 } 00199 // 00200 // Include vertex if required 00201 // 00202 if ( getPara()->vertexConstrainedFit ) { 00203 xold = getPara()->xVertex - xav ; 00204 yold = getPara()->yVertex - yav ; 00205 // 00206 //--> ROTATE SO THAT <XY> = 0 & DIVIDE BY RSCALE SO THAT <R**2> = 1 00207 // 00208 xi = ( cosrot * xold + sinrot * yold ) / rscale ; 00209 yi = ( -sinrot * xold + cosrot * yold ) / rscale ; 00210 00211 xixi = xi * xi ; 00212 yiyi = yi * yi ; 00213 riri = xixi + yiyi ; 00214 wiriri = getPara()->xyWeightVertex * riri ; 00215 00216 xyav += getPara()->xyWeightVertex * xi * yi ; 00217 xxav += getPara()->xyWeightVertex * xixi ; 00218 yyav += getPara()->xyWeightVertex * yiyi ; 00219 00220 xrrav += wiriri * xi ; 00221 yrrav += wiriri * yi ; 00222 rrrrav += wiriri * riri ; 00223 } 00224 // 00225 // 00226 // 00227 //--> DIVIDE BY WSUM TO MAKE AVERAGES 00228 // 00229 xxav = xxav / wsum ; 00230 yyav = yyav / wsum ; 00231 xrrav = xrrav / wsum ; 00232 yrrav = yrrav / wsum ; 00233 rrrrav = rrrrav / wsum ; 00234 xyav = xyav / wsum ; 00235 00236 int const ntry = 5 ; 00237 // 00238 //--> USE THESE TO GET THE COEFFICIENTS OF THE 4-TH ORDER POLYNIMIAL 00239 //--> DON'T PANIC - THE THIRD ORDER TERM IS ZERO ! 00240 // 00241 double xrrxrr = xrrav * xrrav ; 00242 double yrryrr = yrrav * yrrav ; 00243 double rrrrm1 = rrrrav - 1.0 ; 00244 double xxyy = xxav * yyav ; 00245 00246 double c0 = rrrrm1*xxyy - xrrxrr*yyav - yrryrr*xxav ; 00247 double c1 = - rrrrm1 + xrrxrr + yrryrr - 4.0*xxyy ; 00248 double c2 = 4.0 + rrrrm1 - 4.0*xxyy ; 00249 double c4 = - 4.0 ; 00250 // 00251 //--> COEFFICIENTS OF THE DERIVATIVE - USED IN NEWTON-RAPHSON ITERATIONS 00252 // 00253 double c2d = 2.0 * c2 ; 00254 double c4d = 4.0 * c4 ; 00255 // 00256 //--> 0'TH VALUE OF LAMDA - LINEAR INTERPOLATION BETWEEN P(0) & P(YYAV) 00257 // 00258 // LAMDA = YYAV * C0 / (C0 + YRRSQ * (XXAV-YYAV)) 00259 double lamda = 0.0 ; 00260 double dlamda = 0.0 ; 00261 // 00262 double chiscl = wsum * rscale * rscale ; 00263 double dlamax = 0.001 / chiscl ; 00264 00265 double p, pd ; 00266 for ( int itry = 1 ; itry <= ntry ; itry++ ) { 00267 p = c0 + lamda * (c1 + lamda * (c2 + lamda * lamda * c4 )) ; 00268 pd = (c1 + lamda * (c2d + lamda * lamda * c4d)) ; 00269 dlamda = -p / pd ; 00270 lamda = lamda + dlamda ; 00271 if (fabs(dlamda)< dlamax) break ; 00272 } 00273 00274 chi2[0] = (double)(chiscl * lamda) ; 00275 // double dchisq = chiscl * dlamda ; 00276 // 00277 //--> NOW CALCULATE THE MATRIX ELEMENTS FOR ALPHA, BETA & KAPPA 00278 // 00279 double h11 = xxav - lamda ; 00280 double h14 = xrrav ; 00281 double h22 = yyav - lamda ; 00282 double h24 = yrrav ; 00283 double h34 = 1.0 + 2.0*lamda ; 00284 if ( h11 == 0.0 || h22 == 0.0 ){ 00285 ftfLog ( " Problems fitting a circle " ) ; 00286 return 1 ; 00287 } 00288 double rootsq = (h14*h14)/(h11*h11) + 4.0*h34 ; 00289 00290 double ratio, kappa, beta ; 00291 if ( fabs(h22) > fabs(h24) ) { 00292 ratio = h24 / h22 ; 00293 rootsq = ratio * ratio + rootsq ; 00294 kappa = 1.0 / ::sqrt(rootsq) ; 00295 beta = - ratio * kappa ; 00296 } 00297 else { 00298 ratio = h22 / h24 ; 00299 rootsq = 1.0 + ratio * ratio * rootsq ; 00300 beta = 1.0 / ::sqrt(rootsq) ; 00301 if ( h24 > 0 ) beta = - beta ; 00302 kappa = -ratio * beta ; 00303 } 00304 double alpha = - (h14/h11) * kappa ; 00305 // 00306 //--> transform these into the lab coordinate system 00307 //--> first get kappa and back to real dimensions 00308 // 00309 double kappa1 = kappa / rscale ; 00310 double dbro = 0.5 / kappa1 ; 00311 // 00312 //--> next rotate alpha and beta and scale 00313 // 00314 double alphar = (cosrot * alpha - sinrot * beta)* dbro ; 00315 double betar = (sinrot * alpha + cosrot * beta)* dbro ; 00316 // 00317 //--> then translate by (xav,yav) 00318 // 00319 double acent = (double)(xav - alphar) ; 00320 double bcent = (double)(yav - betar ) ; 00321 double radius = (double)dbro ; 00322 // 00323 // Get charge 00324 // 00325 q = ( ( yrrav < 0 ) ? 1 : -1 ) ; 00326 pt = (double)(2.9979e-3 * getPara()->bField * radius ) ; 00327 // 00328 // Get other track parameters 00329 // 00330 double x0, y0 ; 00331 if ( getPara()->vertexConstrainedFit ) { 00332 flag = 1 ; // primary track flag 00333 x0 = getPara()->xVertex ; 00334 y0 = getPara()->yVertex ; 00335 phi0 = getPara()->phiVertex ; 00336 r0 = getPara()->rVertex ; 00337 psi = (double)atan2(bcent-y0,acent-x0) ; 00338 psi = psi + q * 0.5F * pi ; 00339 if ( psi < 0 ) psi = psi + twoPi ; 00340 } 00341 else { 00342 FtfBaseHit* lHit = (FtfBaseHit *)lastHit ; 00343 flag = 0 ; // primary track flag 00344 psi = (double)atan2(bcent-(lHit->y),acent-(lHit->x)) ; 00345 x0 = acent - radius * cos(psi); 00346 y0 = bcent - radius * sin(psi); 00347 psi = psi + q * 0.5F * pi ; 00348 phi0 = atan2(y0,x0); 00349 if ( phi0 < 0 ) phi0 += twoPi ; 00350 r0 = sqrt ( x0 * x0 + y0 * y0 ) ; 00351 } 00352 // 00353 // 00354 // Get errors from fast fit 00355 // 00356 if ( getPara()->getErrors ) getErrorsCircleFit ( acent, bcent, radius ) ; 00357 // 00358 return 0 ; 00359 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00360 // End Fit Circle 00361 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00362 } 00363 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00364 // Fit Line in s-z plane 00365 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00366 int FtfBaseTrack::fitLine ( ) 00367 { 00368 // 00369 // initialization 00370 // 00371 double sum = 0.F ; 00372 double ss = 0.F ; 00373 double sz = 0.F ; 00374 double sss = 0.F ; 00375 double ssz = 0.F ; 00376 // 00377 // find sum , sums ,sumz, sumss 00378 // 00379 double dx, dy ; 00380 double radius = (double)(pt / ( 2.9979e-3 * getPara()->bField ) ) ; 00381 if ( getPara()->vertexConstrainedFit ) { 00382 dx = ((FtfBaseHit *)firstHit)->x - getPara()->xVertex ; 00383 dy = ((FtfBaseHit *)firstHit)->y - getPara()->yVertex ; 00384 } 00385 else { 00386 dx = ((FtfBaseHit *)firstHit)->x - ((FtfBaseHit *)lastHit)->x ; 00387 dy = ((FtfBaseHit *)firstHit)->y - ((FtfBaseHit *)lastHit)->y ; 00388 } 00389 double localPsi = 0.5F * sqrt ( dx*dx + dy*dy ) / radius ; 00390 double total_s ; 00391 if ( fabs(localPsi) < 1. ) { 00392 total_s = 2.0F * radius * asin ( localPsi ) ; 00393 } 00394 else { 00395 // ftfLog ( "FtfBaseTrack: fitLine: problem calculating s \n" ) ; 00396 total_s = 2.0F * radius * M_PI ; 00397 } 00398 00399 // 00400 FtfBaseHit *previousHit = 0 ; 00401 00402 for ( startLoop() ; done() ; nextHit() ) { 00403 FtfBaseHit* cHit = (FtfBaseHit *)currentHit ; 00404 if ( currentHit != firstHit ) { 00405 dx = cHit->x - previousHit->x ; 00406 dy = cHit->y - previousHit->y ; 00407 dpsi = 0.5F * (double)sqrt ( dx*dx + dy*dy ) / radius ; 00408 if ( dpsi > 1.) { 00409 //ftfLog ( "FtfBaseHit::fitLine(): dpsi=%f\n", dpsi); 00410 dpsi = 1.; 00411 } 00412 cHit->s = previousHit->s - 2.0F * radius * (double)asin ( dpsi ) ; 00413 } 00414 else 00415 cHit->s = total_s ; 00416 00417 sum += cHit->wz ; 00418 ss += cHit->wz * cHit->s ; 00419 sz += cHit->wz * cHit->z ; 00420 sss += cHit->wz * cHit->s * cHit->s ; 00421 ssz += cHit->wz * cHit->s * cHit->z ; 00422 previousHit = cHit ; 00423 } 00424 00425 double det = sum * sss - ss * ss; 00426 if ( fabs(det) < 1e-20){ 00427 chi2[1] = 99999.F ; 00428 return 0 ; 00429 } 00430 // 00431 // compute the best fitted parameters A,B 00432 // 00433 tanl = (double)((sum * ssz - ss * sz ) / det ); 00434 z0 = (double)((sz * sss - ssz * ss ) / det ); 00435 // 00436 // calculate chi-square 00437 // 00438 chi2[1] = 0.F ; 00439 double r1 ; 00440 for ( startLoop() ; done() ; nextHit() ) { 00441 FtfBaseHit* cHit = (FtfBaseHit *)currentHit ; 00442 r1 = cHit->z - tanl * cHit->s - z0 ; 00443 chi2[1] += (double) ( (double)cHit->wz * (r1 * r1) ); 00444 } 00445 // 00446 // calculate estimated variance 00447 // varsq=chi/(double(n)-2.) 00448 // calculate covariance matrix 00449 // siga=::sqrt(varsq*sxx/det) 00450 // sigb=::sqrt(varsq*sum/det) 00451 // 00452 dtanl = (double) ( sum / det ); 00453 dz0 = (double) ( sss / det ); 00454 00455 return 0 ; 00456 } 00457 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00458 // End Fit Line 00459 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00460 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00461 // CIRCOV - a covariance matrix calculation program for circle fitting 00462 // DESCRIPTION: 00463 // Compute the covariance matrix of an effective circle fitting algorithm 00464 // The circle equation is (X(I)-A)**2 + (Y(I)-B)**2 = R**2. 00465 // The functional minimum is W(I)*[(X(I)-A)**2+(Y(I)-B)**2-R*R]**2/(R*R) 00466 // For details about the algorithm, see 00467 // N.I. CHERNOV, G.A. OSOSKOV, COMPUT. PHYS. COMMUN. 33(1984) 329-333 00468 // INPUT ARGUMENTS: */ 00469 // A - Best fitted circle center in X axis, REAL 00470 // B - Best fitted circle center in Y axis, REAL 00471 // R - Best fitted radius REAL 00472 // 00473 // From a routine written in Fortran by AUTHOR: 00474 // Jawluen Tang, Physics department , UT-Austin 00475 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00476 int FtfBaseTrack::getErrorsCircleFit ( double a, double b, double r ) { 00477 00478 double h[9] = { 0. }; 00479 double dx, dy ; 00480 double h11, h22, h33 ; 00481 static int j ; 00482 static double ratio, c1, s1; 00483 static double hyp; 00484 00485 00486 for (j = 0; j < 9; j++ ) { 00487 h[j] = 0.; 00488 } 00489 // 00490 // If circle fit was not used the 00491 // errors in the real space need to 00492 // be calculated 00493 // 00494 if ( pt < getPara()->ptMinHelixFit ) { 00495 for ( startLoop() ; done() ; nextHit() ) { 00496 00497 FtfBaseHit* cHit = (FtfBaseHit *)currentHit ; 00498 cHit->wxy = 1.0F/ (double)(cHit->dx*cHit->dx + cHit->dy*cHit->dy) ; 00499 } 00500 } 00501 // 00502 // Loop over points in fit 00503 // 00504 for ( startLoop() ; done() ; nextHit() ) { 00505 FtfBaseHit* cHit = (FtfBaseHit *)currentHit ; 00506 dx = cHit->x - a; 00507 dy = cHit->y - b; 00508 hyp = (double)::sqrt( dx * dx + dy * dy ); 00509 s1 = dx / hyp; 00510 c1 = dy / hyp; 00511 ratio = r / hyp; 00512 h[0] += cHit->wxy * (ratio * (s1 * s1 - 1) + 1); 00513 h[1] += cHit->wxy * ratio * s1 * c1; 00514 h[2] += cHit->wxy * s1; 00515 h[4] += cHit->wxy * (ratio * (c1 * c1 - 1) + 1); 00516 h[5] += cHit->wxy * c1; 00517 h[8] += cHit->wxy ; 00518 } 00519 h[3] = h[1]; 00520 h[6] = h[2]; 00521 h[7] = h[5]; 00522 00523 ftfMatrixDiagonal ( h, h11, h22, h33 ) ; 00524 // 00525 // Calculate pt error now 00526 // 00527 dpt = (double)(2.9979e-3 * getPara()->bField * h33 ); 00528 // 00529 // Get error in psi now 00530 // 00531 if ( getPara()->vertexConstrainedFit ) { 00532 dx = a ; 00533 dy = b ; 00534 } 00535 else { 00536 dx = ((FtfBaseHit *)lastHit)->x + a - ((FtfBaseHit *)firstHit)->x ; 00537 dy = ((FtfBaseHit *)lastHit)->y + b + ((FtfBaseHit *)firstHit)->y ; 00538 } 00539 double w = dy / dx ; 00540 dpsi = (double)(( 1. / ( 1. + w*w ) ) * ( h22 / dx - dy * h11 / ( dx*dx ) )) ; 00541 00542 return 0 ; 00543 } 00544 00545 //************************************************************************* 00546 // Prints one track information 00547 //************************************************************************* 00548 void FtfBaseTrack::Print ( int level ) 00549 { 00550 double pmom, pz; 00551 /* 00552 -----> Print info 00553 */ 00554 if ( level > 9 ) { 00555 pz = pt * tanl ; 00556 pmom = (double)sqrt ( pz * pz + pt * pt ) ; 00557 ftfLog ( " \n =======> Track %d <======== ", id ) ; 00558 ftfLog ( " \n p, pt, q %7.2f %7.2f %2d ", pmom, pt, q ) ; 00559 } 00560 if ( level > 19 ) { 00561 ftfLog ( " \n r0, z0, nHits %7.2f %7.2f %d ", r0, z0, nHits ) ; 00562 ftfLog ( " \n phi0, psi, tanl %7.2f %7.2f %7.2f ", phi0, psi, tanl ) ; 00563 } 00564 else ftfLog ( "\n " ) ; 00565 00566 if ( level > 29 ) { 00567 ftfLog ( " \n chi2 (s,z) %6.2e %6.2e ", chi2[0], 00568 chi2[1] ) ; 00569 } 00570 else ftfLog ( "\n " ) ; 00571 00572 00573 if ( fmod((double)level,10.) > 0 ) { 00574 ftfLog ( " \n *** Clusters in this track *** " ) ; 00575 ((FtfBaseHit *)firstHit)->print ( 10 ) ; 00576 for ( startLoop() ; done() ; nextHit() ) { 00577 ((FtfBaseHit *)currentHit)->print ( 1 ) ; 00578 } 00579 } 00580 ftfLog ( "\n " ) ; 00581 } 00582 /*:>-------------------------------------------------------------------- 00583 **: METHOD: Finds point of closest approach 00584 **: 00585 **: AUTHOR: ppy - P.P. Yepes, yepes@rice.edu 00586 **: ARGUMENTS: 00587 **: IN: xBeam, yBeam: beam position 00588 **: 00589 **: RETURNS: 00590 **: tHit - Point closest approach to center 00591 *:>------------------------------------------------------------------*/ 00592 Ftf3DHit FtfBaseTrack::closestApproach ( double xBeam, double yBeam ) { 00593 double rc, xc, yc ; 00594 return getClosest ( xBeam, yBeam, rc, xc, yc ) ; 00595 } 00596 /*:>-------------------------------------------------------------------- 00597 **: METHOD: Finds point of closest approach 00598 **: 00599 **: AUTHOR: ppy - P.P. Yepes, yepes@rice.edu 00600 **: ARGUMENTS: 00601 **: IN: xBeam, yBeam: beam position 00602 **: OUT: 00603 **: rc, xc, yc track circle radius and center 00604 **: 00605 **: RETURNS: 00606 **: tHit - Point closest approach to center 00607 *:>------------------------------------------------------------------*/ 00608 Ftf3DHit FtfBaseTrack::getClosest ( double xBeam, double yBeam, 00609 double &rc, double &xc, double &yc ) { 00610 double xp, yp, zp ; 00611 xp = yp = 0. ; 00612 //-------------------------------------------------------- 00613 // Get track parameters 00614 //-------------------------------------------------------- 00615 double tPhi0 = psi + double(q) * 0.5 * M_PI / fabs((double)q) ; 00616 00617 double x0 = r0 * cos(phi0) ; 00618 double y0 = r0 * sin(phi0) ; 00619 rc = pt / ( bFactor * getPara()->bField ) ; 00620 xc = x0 - rc * cos(tPhi0) ; 00621 yc = y0 - rc * sin(tPhi0) ; 00622 00623 getClosest ( xBeam, yBeam, rc, xc, yc, xp, yp ) ; 00624 00625 //----------------------------------------------------------------- 00626 // Get the z coordinate 00627 //----------------------------------------------------------------- 00628 double angle = atan2 ( (yp-yc), (xp-xc) ) ; 00629 if ( angle < 0. ) angle = angle + 2.0 * M_PI ; 00630 00631 double dangle = angle - tPhi0 ; 00632 00633 if ( fabs(dangle) < 1.e-4 ) dangle = 0 ; // Problems with -0.000 values 00634 dangle = fmod ( dangle, 2.0 * M_PI ) ; 00635 if ( (float(q) * dangle) < 0 ) 00636 dangle = dangle + float(q) * 2. * M_PI ; 00637 00638 double stot = fabs(dangle) * rc ; 00639 zp = z0 - stot * tanl ; 00640 00641 Ftf3DHit tHit(xp,yp,zp) ; 00642 return tHit ; 00643 } 00644 /*:>-------------------------------------------------------------------- 00645 **: METHOD: Finds point of closest approach 00646 **: 00647 **: AUTHOR: ppy - P.P. Yepes, yepes@rice.edu 00648 **: ARGUMENTS: 00649 **: IN: xBeam, yBeam: beam position 00650 **: rc, xc, yc : track circle radius and center 00651 **: OUT: 00652 **: double xClosest, yClosest 00653 **: 00654 **: RETURNS: 0 if Ok 00655 *:>------------------------------------------------------------------*/ 00656 int FtfBaseTrack::getClosest ( double xBeam, double yBeam, 00657 double rc, double xc, double yc, 00658 double& xClosest, double& yClosest ) { 00659 //---------------------------------------------------------- 00660 // Shift center to respect beam axis 00661 //---------------------------------------------------------- 00662 double dx = xc - xBeam ; 00663 double dy = yc - yBeam ; 00664 //---------------------------------------------------------- 00665 // Solve the equations 00666 //---------------------------------------------------------- 00667 double fact = rc / sqrt ( dx * dx + dy * dy ) ; 00668 double f1 = 1. + fact ; 00669 double f2 = 1. - fact ; 00670 00671 double dx1 = dx * f1 ; 00672 double dy1 = dy * f1 ; 00673 double d1 = sqrt ( dx1 * dx1 + dy1 * dy1 ) ; 00674 00675 double dx2 = dx * f2 ; 00676 double dy2 = dy * f2 ; 00677 double d2 = sqrt ( dx2 * dx2 + dy2 * dy2 ) ; 00678 //--------------------------------------------------------------- 00679 // Choose the closest 00680 //--------------------------------------------------------------- 00681 if ( d1 < d2 ) { 00682 xClosest = dx1 + xBeam ; 00683 yClosest = dy1 + yBeam ; 00684 } 00685 else { 00686 xClosest = dx2 + xBeam ; 00687 yClosest = dy2 + yBeam ; 00688 } 00689 return 0 ; 00690 } 00691 /*:>-------------------------------------------------------------------- 00692 **: METHOD: Extrapolates track to cylinder with radius r 00693 **: 00694 **: 00695 **: AUTHOR: ppy - P.P. Yepes, yepes@rice.edu 00696 **: ARGUMENTS: 00697 **: IN: 00698 **: track - Global track pointer 00699 **: r - Cylinder radius 00700 **: OUT: 00701 **: x,y,z - Extrapolated point 00702 **: xc,yc,rr - Center and radius track circle in x-y plane 00703 **: 00704 **: RETURNS: 0=ok, <>0 error 00705 **:>------------------------------------------------------------------*/ 00706 Ftf3DHit FtfBaseTrack::extraRadius ( double r ) { 00707 double phi ; 00708 // 00709 // Default values 00710 // 00711 double x, y, z, rc, xc, yc ; 00712 x = y = z = 0.F ; 00713 // 00714 // If error return with error 00715 // 00716 Ftf3DHit tHit(0,0,0) ; 00717 if ( extraRCyl ( r, phi, z, rc, xc, yc ) ) return tHit ; 00718 // 00719 // Otherwise get point in cartesian coordinates 00720 // 00721 x = r * cos(phi) ; 00722 y = r * sin(phi) ; 00723 tHit.x = x ; 00724 tHit.y = y ; 00725 tHit.z = z ; 00726 00727 return tHit ; 00728 } 00729 /*:>-------------------------------------------------------------------- 00730 **: METHOD: Extrapolates track to cylinder with radius r 00731 **: 00732 **: 00733 **: AUTHOR: ppy - P.P. Yepes, yepes@rice.edu 00734 **: ARGUMENTS: 00735 **: IN: 00736 **: r - Cylinder radius 00737 **: OUT: 00738 **: phi,z - Extrapolated point 00739 **: xc,yc,rc - Center and radius track circle in x-y plane 00740 **: 00741 **: RETURNS: 0=ok, <>0 error 00742 **:>------------------------------------------------------------------*/ 00743 int FtfBaseTrack::extraRCyl ( double &r, double &phi, double &z, 00744 double &rc, double &xc, double &yc ) { 00745 00746 double td ; 00747 double fac1,sfac, fac2 ; 00748 //-------------------------------------------------------- 00749 // Get track parameters 00750 //-------------------------------------------------------- 00751 double tPhi0 = psi + double(q) * 0.5 * M_PI / fabs((double)q) ; 00752 double x0 = r0 * cos(phi0) ; 00753 double y0 = r0 * sin(phi0) ; 00754 rc = fabs(pt) / ( bFactor * getPara()->bField ) ; 00755 xc = x0 - rc * cos(tPhi0) ; 00756 yc = y0 - rc * sin(tPhi0) ; 00757 // 00758 // Check helix and cylinder intersect 00759 // 00760 fac1 = xc*xc + yc*yc ; 00761 sfac = ::sqrt( fac1 ) ; 00762 // 00763 // If they don't intersect return 00764 // Trick to solve equation of intersection of two circles 00765 // rotate coordinates to have both circles with centers on x axis 00766 // pretty simple system of equations, then rotate back 00767 // 00768 if ( fabs(sfac-rc) > r || fabs(sfac+rc) < r ) { 00769 // ftfLog ( "particle does not intersect \n" ) ; 00770 return 1 ; 00771 } 00772 // 00773 // Find intersection 00774 // 00775 fac2 = ( r*r + fac1 - rc*rc) / (2.00 * r * sfac ) ; 00776 phi = atan2(yc,xc) + float(q)*acos(fac2) ; 00777 td = atan2(r*sin(phi) - yc,r*cos(phi) - xc) ; 00778 // Intersection in z 00779 00780 if ( td < 0 ) td = td + 2. * M_PI ; 00781 double dangle = tPhi0 - td ; 00782 dangle = fmod ( dangle, 2.0 * M_PI ) ; 00783 if ( r < r0 ) dangle *= -1 ; 00784 // ftfLog ( "dangle %f q %d \n", dangle, q ) ; 00785 if ( (float(q) * dangle) < 0 ) dangle = dangle + float(q) * 2. * M_PI ; 00786 00787 double stot = fabs(dangle) * rc ; 00788 // ftfLog ( "dangle %f z0 %f stot %f \n", dangle, z0, stot ) ; 00789 if ( r > r0 ) z = z0 + stot * tanl ; 00790 else z = z0 - stot * tanl ; 00791 // 00792 // That's it 00793 // 00794 return 0 ; 00795 } 00796 /*:>-------------------------------------------------------------------- 00797 **: METHOD: Calculates intersection of track with plane define by line 00798 **: y = a x + b and the z 00799 **: 00800 **: AUTHOR: ppy - P.P. Yepes, yepes@rice.edu 00801 **: ARGUMENTS: 00802 **: IN: 00803 **: a, b - Line parameters 00804 **: OUT: 00805 **: crossPoint - intersection point 00806 **: 00807 **: RETURNS: 0=ok, <>0 track does not cross the plane 00808 **:>------------------------------------------------------------------*/ 00809 00810 //-------------------------------------------------------------------- 00811 // tom: split into two methods: 00812 // one gets both intersections within one turn of the helix, 00813 // the other selects the intersection closest to the origin. 00814 // (Provided for compatibility with older programs.) 00815 //-------------------------------------------------------------------- 00816 00817 int FtfBaseTrack::intersectorZLine ( double a, double b, Ftf3DHit& cross ) { 00818 Ftf3DHit cross1, cross2; 00819 00820 if (0 != intersectorZLine(a, b, cross1, cross2)) 00821 return 1; 00822 00823 double r1sq = cross1.x*cross1.x + cross1.y*cross1.y; 00824 double r2sq = cross2.x*cross2.x + cross2.y*cross2.y; 00825 00826 if (r1sq < r2sq) { // cross 1 is closer to the beamline 00827 cross = cross1; 00828 } else { 00829 cross = cross2; 00830 } 00831 00832 return 0; 00833 } 00834 00835 int FtfBaseTrack::intersectorZLine ( double a, double b, 00836 Ftf3DHit& cross1, Ftf3DHit& cross2 ) { 00837 // 00838 // Calculate circle center and radius 00839 // 00840 double x0 = r0 * cos(phi0) ; // x first point on track 00841 double y0 = r0 * sin(phi0) ; // y first point on track 00842 double trackPhi0 = psi + q * 0.5 * M_PI / fabs((double)q) ;// 00843 double rc = pt / ( bFactor * getPara()->bField ) ; 00844 double xc = x0 - rc * cos(trackPhi0) ; 00845 double yc = y0 - rc * sin(trackPhi0) ; 00846 00847 double ycPrime = yc - b ; 00848 double aa = ( 1. + a * a ) ; 00849 double bb = -2. * ( xc + a * ycPrime ) ; 00850 double cc = ( xc * xc + ycPrime * ycPrime - rc * rc ) ; 00851 00852 double racine = bb * bb - 4. * aa * cc ; 00853 if ( racine < 0 ) { 00854 cross1.set(0,0,0); 00855 cross2.set(0,0,0); 00856 return 1 ; 00857 } 00858 double rootRacine = ::sqrt(racine) ; 00859 00860 double oneOverA = 1./aa; 00861 // 00862 // First solution 00863 // 00864 double x1 = 0.5 * oneOverA * ( -1. * bb + rootRacine ) ; 00865 double y1 = a * x1 + b ; 00866 // 00867 // Second solution 00868 // 00869 double x2 = 0.5 * oneOverA * ( -1. * bb - rootRacine ) ; 00870 double y2 = a * x2 + b ; 00871 00872 00873 //------------------------------------------------------------------- 00874 // Get the first z coordinate 00875 //------------------------------------------------------------------- 00876 double angle, dangle, stot; 00877 00878 angle = atan2 ( (y1-yc), (x1-xc) ) ; 00879 if ( angle < 0. ) angle = angle + 2.0 * M_PI ; 00880 00881 dangle = angle - trackPhi0 ; 00882 dangle = fmod ( dangle, 2.0 * M_PI ) ; 00883 00884 if ( (q * dangle) > 0 ) dangle = dangle - q * 2. * M_PI ; 00885 00886 stot = fabs(dangle) * rc ; 00887 double z1 = z0 + stot * tanl ; 00888 00889 cross1.set ( x1, y1, z1 ) ; 00890 00891 00892 //------------------------------------------------------------------- 00893 // Get the second z coordinate 00894 //------------------------------------------------------------------- 00895 angle = atan2 ( (y2-yc), (x2-xc) ) ; 00896 if ( angle < 0. ) angle = angle + 2.0 * M_PI ; 00897 00898 dangle = angle - trackPhi0 ; 00899 dangle = fmod ( dangle, 2.0 * M_PI ) ; 00900 00901 if ( (q * dangle) > 0 ) dangle = dangle - q * 2. * M_PI ; 00902 00903 stot = fabs(dangle) * rc ; 00904 double z2 = z0 + stot * tanl ; 00905 00906 cross2.set ( x2, y2, z2 ) ; 00907 00908 return 0 ; 00909 } 00910 00911 00912 00913 /*:>-------------------------------------------------------------------- 00914 **: METHOD: Calculates intersection of track with plane define by line 00915 **: x = a and the z 00916 **: 00917 **: AUTHOR: ppy - P.P. Yepes, yepes@rice.edu 00918 **: ARGUMENTS: 00919 **: IN: 00920 **: a - Line parameter 00921 **: OUT: 00922 **: crossPoint - intersection point 00923 **: 00924 **: RETURNS: 0=ok, <>0 track does not cross the plane 00925 **:>------------------------------------------------------------------*/ 00926 int FtfBaseTrack::intersectorYCteLine ( double a, Ftf3DHit& cross ) { 00927 // 00928 //calculate circle center and radius 00929 // 00930 double x0 = r0*cos(phi0); 00931 double y0 = r0*sin(phi0); 00932 double trackPhi0= psi + q*0.5*M_PI/fabs(q); 00933 double rcoc = pt / ( bFactor * getPara()->bField ) ; 00934 double xcoc = x0 - (rcoc*cos(trackPhi0)); 00935 double ycoc = y0 - (rcoc*sin(trackPhi0)); 00936 // 00937 // Calculate circle center and radius 00938 // 00939 double xHit = a ; 00940 00941 double f1 = (xHit-xcoc)*(xHit-xcoc); 00942 double r2 = rcoc*rcoc; 00943 if ( f1 > r2 ) { 00944 return 1 ; 00945 } 00946 00947 double sf2 = ::sqrt(r2-f1); 00948 double y1 = ycoc + sf2; 00949 double y2 = ycoc - sf2; 00950 double yHit = y1; 00951 if ( fabs(y2) < fabs(y1) ) yHit=y2; 00952 00953 //Get z coordinate: 00954 double angle = atan2 ( (yHit-ycoc), (xHit-xcoc) ) ; 00955 if ( angle < 0. ) angle = angle + 2.0 * M_PI ; 00956 00957 double dangle = angle - trackPhi0 ; 00958 dangle = fmod ( dangle, 2.0 * M_PI ) ; 00959 if ( (q * dangle) > 0 ) dangle = dangle - q * 2. * M_PI ; 00960 00961 double stot = fabs(dangle) * rcoc ; 00962 double zHit = z0 + stot * tanl; 00963 00964 cross.set(xHit,yHit,zHit); 00965 // 00966 return 0 ; 00967 } 00968 // 00969 /*:>-------------------------------------------------------------------- 00970 **: METHOD: Calculates trajectory length between two points on a track 00971 **: 00972 **: AUTHOR: ppy - P.P. Yepes, yepes@rice.edu 00973 **: ARGUMENTS: 00974 **: IN: 00975 **: track - Track pointer 00976 **: x1, y1 - Point 1 00977 **: x2, y2 - Point 2 00978 **: OUT: 00979 **: 00980 **: RETURNS: 0=ok, <>0 error 00981 **:>------------------------------------------------------------------*/ 00982 double FtfBaseTrack::arcLength ( double x1, double y1, double x2, double y2 ) 00983 { 00984 double x0, y0, xc, yc, rc ; 00985 double angle_1, angle_2, d_angle, sleng_xy, sleng ; 00986 /*---------------------------------------------------------- 00987 Get track parameters 00988 ----------------------------------------------------------*/ 00989 00990 x0 = r0 * cos(phi0) ; 00991 y0 = r0 * sin(phi0) ; 00992 rc = pt / ( bFactor * getPara()->bField ) ; 00993 double tPhi0 = psi + double(q) * 0.5 * M_PI / fabs((double)q) ; 00994 xc = x0 - rc * cos(tPhi0) ; 00995 yc = y0 - rc * sin(tPhi0) ; 00996 /* 00997 Get angle difference 00998 */ 00999 angle_1 = atan2 ( (y1-yc), (x1-xc) ) ; 01000 if ( angle_1 < 0 ) angle_1 = angle_1 + 2. * M_PI ; 01001 angle_2 = atan2 ( (y2-yc), (x2-xc) ) ; 01002 if ( angle_2 < 0 ) angle_2 = angle_2 + 2. * M_PI ; 01003 d_angle = double(q) * ( angle_1 - angle_2 ) ; 01004 d_angle = fmod ( d_angle, 2. * M_PI ) ; 01005 if ( d_angle < 0 ) d_angle = d_angle + 2. * M_PI ; 01006 /*---------------------------------------------------------- 01007 Get total angle and total trajectory 01008 ----------------------------------------------------------*/ 01009 sleng_xy = fabs ( rc ) * d_angle ; 01010 sleng = sleng_xy * sqrt ( 1.0 + tanl * tanl ) ; 01011 return sleng ; 01012 01013 } 01014 /*:>-------------------------------------------------------------------- 01015 **: METHOD: Phi rotates the track 01016 **: 01017 **: AUTHOR: ppy - P.P. Yepes, yepes@rice.edu 01018 **: ARGUMENTS: 01019 **: IN: 01020 **: deltaPhi - Angle to rotate in phi 01021 **: 01022 **: RETURNS: 0=ok, <>0 error 01023 **:>------------------------------------------------------------------*/ 01024 int FtfBaseTrack::phiRotate ( double deltaPhi ) { 01025 phi0 += deltaPhi ; 01026 if ( phi0 > 2. * M_PI ) phi0 -= 2. * M_PI ; 01027 if ( phi0 < 0 ) phi0 += 2. * M_PI ; 01028 psi += deltaPhi ; 01029 if ( psi > 2. * M_PI ) psi -= 2. * M_PI ; 01030 if ( psi < 0 ) psi += 2. * M_PI ; 01031 01032 return 0 ; 01033 } 01034 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 01035 // Fit a circle 01036 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 01037 int FtfBaseTrack::refitHelix ( int mod, int modEqual, int rowMin, int rowMax ) { 01038 typedef FtfBaseHit* hitPointer; 01039 01040 if ( nHits < 1 || nHits > 500 ) { 01041 ftfLog ( "FtfBaseTrack:refitHelix: nHits %d out of range \n", nHits ) ; 01042 return 1 ; 01043 } 01044 // 01045 hitPointer* hitP = new hitPointer[nHits]; 01046 int nHitsOrig = nHits ; 01047 01048 int counter = 0 ; 01049 for ( startLoop() ; done() ; nextHit() ) { 01050 hitP[counter] = (FtfBaseHit *)currentHit ; 01051 // hitP[counter]->print(1); 01052 counter++ ; 01053 } 01054 01055 int row ; 01056 firstHit = 0 ; 01057 counter = 0 ; 01058 for ( int i = 0 ; i < nHitsOrig ; i++ ) { 01059 row = hitP[i]->row ; 01060 hitP[i]->nextTrackHit = 0 ; 01061 if ( row%mod != modEqual ) continue ; 01062 if ( row < rowMin || row > rowMax ) continue ; 01063 01064 if ( firstHit == 0 ) firstHit = hitP[i] ; 01065 else 01066 ((FtfBaseHit *)lastHit)->nextTrackHit = (void *)hitP[i] ; 01067 lastHit = (void *)hitP[i] ; 01068 counter++ ; 01069 } 01070 nHits = counter ; 01071 01072 int problem = 0 ; 01073 if ( nHits > 5 ) fitHelix ( ) ; 01074 else problem = 1 ; 01075 // 01076 // Put hits back 01077 // 01078 firstHit = 0 ; 01079 for ( int i = 0 ; i < nHitsOrig ; i++ ) { 01080 row = hitP[i]->row ; 01081 if ( firstHit == 0 ) firstHit = hitP[i] ; 01082 else 01083 ((FtfBaseHit *)lastHit)->nextTrackHit = (void *)hitP[i] ; 01084 lastHit = (void *)hitP[i] ; 01085 } 01086 nHits = nHitsOrig ; 01087 01088 delete []hitP ; 01089 01090 return problem ; 01091 } 01092 /*:>-------------------------------------------------------------------- 01093 **: METHOD: Updates track parameters to point of intersection with 01094 **: cylinder of radius r 01095 **: 01096 **: 01097 **: AUTHOR: ppy - P.P. Yepes, yepes@rice.edu 01098 **: ARGUMENTS: 01099 **: IN: 01100 **: radius - Cylinder radius to extrapolate track 01101 **: OUT: 01102 **: 01103 **:>------------------------------------------------------------------*/ 01104 void FtfBaseTrack::updateToRadius ( double radius ) { 01105 01106 double phiExtra, zExtra, rCircle, xCircleCenter, yCircleCenter ; 01107 01108 int ok = extraRCyl ( radius, phiExtra, zExtra, rCircle, xCircleCenter, yCircleCenter ) ; 01109 if ( ok ) { 01110 // ftfLog ( "FtfBaseTrack::updateToRadius: track %d does not intersect radius %f\n", 01111 // id, radius ) ; 01112 return ; 01113 } 01114 // 01115 // Check extrapolation falls inside volume 01116 // 01117 // if ( fabs(zExtra) > getPara()->zMax ) { 01118 // ftfLog ( "problem extrapolating track \n" ) ; 01119 // return ; 01120 // } 01121 double xExtra = radius * cos(phiExtra) ; 01122 double yExtra = radius * sin(phiExtra) ; 01123 01124 double tPhi = atan2(yExtra-yCircleCenter,xExtra-xCircleCenter); 01125 01126 // if ( tPhi < 0 ) tPhi += 2. * M_PI ; 01127 01128 double tPsi = tPhi - double(q) * 0.5 * M_PI / fabs((double)q) ; 01129 if ( tPsi > 2. * M_PI ) tPsi -= 2. * M_PI ; 01130 if ( tPsi < 0. ) tPsi += 2. * M_PI ; 01131 // 01132 // Update track parameters 01133 // 01134 r0 = radius ; 01135 phi0 = phiExtra ; 01136 z0 = zExtra ; 01137 psi = tPsi ; 01138 } 01139 /*:>-------------------------------------------------------------------- 01140 **: METHOD: Updates track parameters to point of closest approach 01141 **: 01142 **: 01143 **: AUTHOR: ppy - P.P. Yepes, yepes@rice.edu 01144 **: ARGUMENTS: 01145 **: IN: 01146 **: xBeam - x Beam axis 01147 **: yBeam - y Beam axis 01148 **: rMax - radius of point of closest approach beyond which no update 01149 **: 01150 **:>------------------------------------------------------------------*/ 01151 void FtfBaseTrack::updateToClosestApproach ( double xBeam, double yBeam, double rMax) { 01152 double rc, xc, yc ; 01153 Ftf3DHit closest = getClosest ( xBeam, yBeam, rc, xc, yc ) ; 01154 // ftfLog ( "FtfBaseTrack::updateClosestApproach: closest x y z %f %f %f \n", 01155 // closest.x, closest.y, closest.z ) ; 01156 01157 double radius = ::sqrt(closest.x*closest.x+closest.y*closest.y); 01158 if ( radius > rMax ) return ; 01159 // 01160 double tPhi = atan2(closest.y-yc,closest.x-xc); 01161 01162 // if ( tPhi < 0 ) tPhi += 2. * M_PI ; 01163 01164 double tPsi = tPhi - double(q) * 0.5 * M_PI / fabs((double)q) ; 01165 if ( tPsi > 2. * M_PI ) tPsi -= 2. * M_PI ; 01166 if ( tPsi < 0. ) tPsi += 2. * M_PI ; 01167 // 01168 // Update track parameters 01169 // 01170 r0 = radius ; 01171 phi0 = atan2(closest.y,closest.x) ; 01172 if ( phi0 < 0 ) phi0 += 2. * M_PI ; 01173 z0 = closest.z ; 01174 psi = tPsi ; 01175 } 01176 01177 01178 /*:>-------------------------------------------------------------------- 01179 **: METHOD: Extrapolation of a track to the given pathlenght 01180 **: 01181 **: 01182 **: 01183 **: AUTHOR: JB, Jens Berger, jberger@bnl.gov 01184 **: ARGUMENTS: 01185 **: IN: 01186 **: doubel pathlength: pathlength of the track, staring at x0,y0,z0 01187 **: OUT: 01188 **: one of pablos fancy Ftf3DHit 01189 **: 01190 **:>------------------------------------------------------------------*/ 01191 Ftf3DHit FtfBaseTrack::extrapolate2PathLength(double pathlength) 01192 { 01193 // some helix para 01194 01195 // BFilef scaled with 0.01 01196 double Bfield=getPara()->bField * 0.01; 01197 double c=0.3; 01198 double lambda=atan(tanl); 01199 double kapa=(c*q*Bfield)/pt; 01200 double heli=-((q*Bfield)/fabs(q*Bfield)); 01201 double phi=psi-heli*M_PI/2; 01202 // staring point 01203 double x0=r0*cos(phi0); 01204 double y0=r0*sin(phi0); 01205 01206 Ftf3DHit CoordOfExtrapol; 01207 CoordOfExtrapol.x = x0 + (1/kapa)*(cos(phi+heli*pathlength*kapa*cos(lambda))-cos(phi)); 01208 CoordOfExtrapol.y = y0 + (1/kapa)*(sin(phi+heli*pathlength*kapa*cos(lambda))-sin(phi)); 01209 CoordOfExtrapol.z = z0 + pathlength*sin(lambda); 01210 01211 // debug 01212 //cout<<"x0:"<<x0<<" y0:"<<y0<<" z0:"<<" lambda:"<<lambda<<z0<<" phi:"<<phi<<" h:"<<heli<<" kapa:"<<kapa<<endl; 01213 //cout<<" xs:"<<CoordOfExtrapol.x<<" ys:"<<CoordOfExtrapol.y<<" zs:"<<CoordOfExtrapol.z<<endl; 01214 01215 01216 // return coord of extrapolation 01217 return CoordOfExtrapol; 01218 } 01219 01220 /*:>-------------------------------------------------------------------- 01221 **: METHOD: Get the Radius of the track, helix at X,Y center 01222 **: 01223 **: 01224 **: 01225 **: AUTHOR: JB, Jens Berger, jberger@bnl.gov 01226 **: ARGUMENTS: 01227 **: IN: 01228 **: nada 01229 **: OUT: 01230 **: radius 01231 **: 01232 **:>------------------------------------------------------------------*/ 01233 double FtfBaseTrack::getRadius() 01234 { 01235 double radius=pt / ( bFactor * getPara()->bField ); 01236 01237 return radius; 01238 } 01239 01240 /*:>-------------------------------------------------------------------- 01241 **: METHOD: Get y center of the corresponding radius of the track, helix 01242 **: 01243 **: 01244 **: 01245 **: AUTHOR: JB, Jens Berger, jberger@bnl.gov 01246 **: ARGUMENTS: 01247 **: IN: 01248 **: nada 01249 **: OUT: 01250 **: x center 01251 **: 01252 **:>------------------------------------------------------------------*/ 01253 double FtfBaseTrack::getXCenter() 01254 { 01255 double tPhi0 = psi + double(q) * 0.5 * M_PI / fabs((double)q) ; 01256 // staring point 01257 double x0=r0*cos(phi0); 01258 01259 return ( x0- getRadius() * cos(tPhi0) ); 01260 } 01261 /*:>-------------------------------------------------------------------- 01262 **: METHOD: Get y center of the corresponding radius of the track, helix 01263 **: 01264 **: 01265 **: 01266 **: AUTHOR: JB, Jens Berger, jberger@bnl.gov 01267 **: ARGUMENTS: 01268 **: IN: 01269 **: nada 01270 **: OUT: 01271 **: y center 01272 **: 01273 **:>------------------------------------------------------------------*/ 01274 double FtfBaseTrack::getYCenter() 01275 { 01276 double tPhi0 = psi + double(q) * 0.5 * M_PI / fabs((double)q) ; 01277 // staring point 01278 double y0=r0*sin(phi0); 01279 01280 return ( y0 - getRadius() * sin(tPhi0) ); 01281 } 01282