| d54804bf |
1 | // $Id$ |
| 2 | //*************************************************************************** |
| 3 | // This file is property of and copyright by the ALICE HLT Project * |
| 4 | // ALICE Experiment at CERN, All rights reserved. * |
| 5 | // * |
| 6 | // Primary Authors: Sergey Gorbunov <sergey.gorbunov@kip.uni-heidelberg.de> * |
| 7 | // Ivan Kisel <kisel@kip.uni-heidelberg.de> * |
| 8 | // for The ALICE HLT Project. * |
| 9 | // * |
| 10 | // Permission to use, copy, modify and distribute this software and its * |
| 11 | // documentation strictly for non-commercial purposes is hereby granted * |
| 12 | // without fee, provided that the above copyright notice appears in all * |
| 13 | // copies and that both the copyright notice and this permission notice * |
| 14 | // appear in the supporting documentation. The authors make no claims * |
| 15 | // about the suitability of this software for any purpose. It is * |
| 16 | // provided "as is" without express or implied warranty. * |
| 17 | //*************************************************************************** |
| 18 | |
| 19 | #include "AliHLTTPCCATrackParam.h" |
| 20 | #include "TMath.h" |
| eb30eb49 |
21 | #include "Riostream.h" |
| d54804bf |
22 | |
| 23 | //ClassImp(AliHLTTPCCATrackParam) |
| 24 | |
| 25 | // |
| 26 | // Circle in XY: |
| 27 | // |
| 28 | // R = 1/TMath::Abs(Kappa); |
| 29 | // Xc = X - sin(Phi)/Kappa; |
| 30 | // Yc = Y + cos(Phi)/Kappa; |
| 31 | // |
| 32 | |
| 33 | |
| 34 | |
| 35 | void AliHLTTPCCATrackParam::ConstructXY3( const Float_t x[3], const Float_t y[3], |
| 36 | const Float_t sigmaY2[3], Float_t CosPhi0 ) |
| 37 | { |
| 38 | //* Construct the track in XY plane by 3 points |
| 39 | |
| 40 | Float_t x0 = x[0]; |
| 41 | Float_t y0 = y[0]; |
| 42 | Float_t x1 = x[1] - x0; |
| 43 | Float_t y1 = y[1] - y0; |
| 44 | Float_t x2 = x[2] - x0; |
| 45 | Float_t y2 = y[2] - y0; |
| 46 | |
| 47 | Float_t a1 = x1*x1 + y1*y1; |
| 48 | Float_t a2 = x2*x2 + y2*y2; |
| 49 | Float_t a = 2*(x1*y2 - y1*x2); |
| 50 | Float_t lx = a1*y2 - a2*y1; |
| 51 | Float_t ly = -a1*x2 + a2*x1; |
| 52 | Float_t l = TMath::Sqrt(lx*lx + ly*ly); |
| 53 | |
| 54 | Float_t li = 1./l; |
| 55 | Float_t li2 = li*li; |
| 56 | Float_t li3 = li2*li; |
| 57 | Float_t cosPhi = ly*li; |
| 58 | |
| 59 | Float_t sinPhi = -lx*li; |
| 60 | Float_t kappa = a/l; |
| 61 | |
| 62 | Float_t dlx = a2 - a1; // D lx / D y0 |
| 63 | Float_t dly = -a; // D ly / D y0 |
| 64 | Float_t dA = 2*(x2 - x1); // D a / D y0 |
| 65 | Float_t dl = (lx*dlx + ly*dly)*li; |
| 66 | |
| 67 | // D sinPhi,kappa / D y0 |
| 68 | |
| 69 | Float_t d0[2] = { -(dlx*ly-lx*dly)*ly*li3, (dA*l-a*dl)*li2 }; |
| 70 | |
| 71 | // D sinPhi,kappa / D y1 |
| 72 | |
| 73 | dlx = -a2 + 2*y1*y2; |
| 74 | dly = -2*x2*y1; |
| 75 | dA = -2*x2; |
| 76 | dl = (lx*dlx + ly*dly)*li; |
| 77 | |
| 78 | Float_t d1[2] = { -(dlx*ly-lx*dly)*ly*li3, (dA*l-a*dl)*li2 }; |
| 79 | |
| 80 | // D sinPhi,kappa / D y2 |
| 81 | |
| 82 | dlx = a1 - 2*y1*y2; |
| 83 | dly = -2*x1*y2; |
| 84 | dA = 2*x1; |
| 85 | dl = (lx*dlx + ly*dly)*li; |
| 86 | |
| 87 | Float_t d2[2] = { -(dlx*ly-lx*dly)*ly*li3, (dA*l-a*dl)*li2 }; |
| 88 | |
| 89 | if( CosPhi0*cosPhi <0 ){ |
| 90 | cosPhi = -cosPhi; |
| 91 | sinPhi = -sinPhi; |
| 92 | kappa = -kappa; |
| 93 | d0[0] = -d0[0]; |
| 94 | d0[1] = -d0[1]; |
| 95 | d1[0] = -d1[0]; |
| 96 | d1[1] = -d1[1]; |
| 97 | d2[0] = -d2[0]; |
| 98 | d2[1] = -d2[1]; |
| 99 | } |
| 100 | |
| 101 | X() = x0; |
| 102 | Y() = y0; |
| 103 | SinPhi() = sinPhi; |
| 104 | Kappa() = kappa; |
| 105 | CosPhi() = cosPhi; |
| 106 | |
| 107 | Float_t s0 = sigmaY2[0]; |
| 108 | Float_t s1 = sigmaY2[1]; |
| 109 | Float_t s2 = sigmaY2[2]; |
| eb30eb49 |
110 | |
| d54804bf |
111 | fC[0] = s0; |
| 112 | fC[1] = 0; |
| eb30eb49 |
113 | fC[2] = 100.; |
| d54804bf |
114 | |
| 115 | fC[3] = d0[0]*s0; |
| 116 | fC[4] = 0; |
| 117 | fC[5] = d0[0]*d0[0]*s0 + d1[0]*d1[0]*s1 + d2[0]*d2[0]*s2; |
| 118 | |
| 119 | fC[6] = 0; |
| 120 | fC[7] = 0; |
| 121 | fC[8] = 0; |
| eb30eb49 |
122 | fC[9] = 100.; |
| d54804bf |
123 | |
| 124 | fC[10] = d0[1]*s0; |
| 125 | fC[11] = 0; |
| 126 | fC[12] = d0[0]*d0[1]*s0 + d1[0]*d1[1]*s1 + d2[0]*d2[1]*s2; |
| 127 | fC[13] = 0; |
| 128 | fC[14] = d0[1]*d0[1]*s0 + d1[1]*d1[1]*s1 + d2[1]*d2[1]*s2; |
| 129 | } |
| 130 | |
| 131 | |
| 132 | Float_t AliHLTTPCCATrackParam::GetS( Float_t x, Float_t y ) const |
| 133 | { |
| 134 | //* Get XY path length to the given point |
| 135 | |
| 136 | Float_t k = GetKappa(); |
| 137 | Float_t ex = GetCosPhi(); |
| 138 | Float_t ey = GetSinPhi(); |
| 139 | x-= GetX(); |
| 140 | y-= GetY(); |
| 141 | Float_t dS = x*ex + y*ey; |
| 142 | if( TMath::Abs(k)>1.e-4 ) dS = TMath::ATan2( k*dS, 1+k*(x*ey-y*ex) )/k; |
| 143 | return dS; |
| 144 | } |
| 145 | |
| 146 | void AliHLTTPCCATrackParam::GetDCAPoint( Float_t x, Float_t y, Float_t z, |
| 147 | Float_t &xp, Float_t &yp, Float_t &zp ) const |
| 148 | { |
| 149 | //* Get the track point closest to the (x,y,z) |
| 150 | |
| 151 | Float_t x0 = GetX(); |
| 152 | Float_t y0 = GetY(); |
| 153 | Float_t k = GetKappa(); |
| 154 | Float_t ex = GetCosPhi(); |
| 155 | Float_t ey = GetSinPhi(); |
| 156 | Float_t dx = x - x0; |
| 157 | Float_t dy = y - y0; |
| 158 | Float_t ax = dx*k+ey; |
| 159 | Float_t ay = dy*k-ex; |
| 160 | Float_t a = sqrt( ax*ax+ay*ay ); |
| 161 | xp = x0 + (dx - ey*( (dx*dx+dy*dy)*k - 2*(-dx*ey+dy*ex) )/(a+1) )/a; |
| 162 | yp = y0 + (dy + ex*( (dx*dx+dy*dy)*k - 2*(-dx*ey+dy*ex) )/(a+1) )/a; |
| 163 | Float_t s = GetS(x,y); |
| 164 | zp = GetZ() + GetDzDs()*s; |
| 165 | if( TMath::Abs(k)>1.e-2 ){ |
| 166 | Float_t dZ = TMath::Abs( GetDzDs()*TMath::TwoPi()/k ); |
| 167 | if( dZ>.1 ){ |
| 168 | zp+= TMath::Nint((z-zp)/dZ)*dZ; |
| 169 | } |
| 170 | } |
| 171 | } |
| 172 | |
| 173 | void AliHLTTPCCATrackParam::ConstructXYZ3( const Float_t p0[5], const Float_t p1[5], |
| 174 | const Float_t p2[5], |
| 175 | Float_t CosPhi0, Float_t t0[] ) |
| 176 | { |
| 177 | //* Construct the track in XYZ by 3 points |
| 178 | |
| 179 | Float_t px[3] = { p0[0], p1[0], p2[0] }; |
| 180 | Float_t py[3] = { p0[1], p1[1], p2[1] }; |
| 181 | Float_t pz[3] = { p0[2], p1[2], p2[2] }; |
| 182 | Float_t ps2y[3] = { p0[3]*p0[3], p1[3]*p1[3], p2[3]*p2[3] }; |
| 183 | Float_t ps2z[3] = { p0[4]*p0[4], p1[4]*p1[4], p2[4]*p2[4] }; |
| 184 | |
| 185 | Float_t kold = t0 ?t0[4] :0; |
| 186 | ConstructXY3( px, py, ps2y, CosPhi0 ); |
| 187 | |
| 188 | Float_t pS[3] = { GetS(px[0],py[0]), GetS(px[1],py[1]), GetS(px[2],py[2]) }; |
| 189 | Float_t k = Kappa(); |
| 190 | if( TMath::Abs(k)>1.e-2 ){ |
| 191 | Float_t dS = TMath::Abs( TMath::TwoPi()/k ); |
| 192 | pS[1]+= TMath::Nint( (pS[0]-pS[1])/dS )*dS; // not more than half turn |
| 193 | pS[2]+= TMath::Nint( (pS[1]-pS[2])/dS )*dS; |
| 194 | if( t0 ){ |
| 195 | Float_t dZ = TMath::Abs(t0[3]*dS); |
| 196 | if( TMath::Abs(dZ)>1. ){ |
| 197 | Float_t dsDz = 1./t0[3]; |
| 198 | if( kold*k<0 ) dsDz = -dsDz; |
| 199 | Float_t s0 = (pz[0]-t0[1])*dsDz; |
| 200 | Float_t s1 = (pz[1]-t0[1])*dsDz; |
| 201 | Float_t s2 = (pz[2]-t0[1])*dsDz; |
| 202 | pS[0]+= TMath::Nint( (s0-pS[0])/dS )*dS ; |
| 203 | pS[1]+= TMath::Nint( (s1-pS[1])/dS )*dS ; |
| 204 | pS[2]+= TMath::Nint( (s2-pS[2])/dS )*dS ; |
| 205 | } |
| 206 | } |
| 207 | } |
| 208 | |
| 209 | Float_t s = pS[0] + pS[1] + pS[2]; |
| 210 | Float_t z = pz[0] + pz[1] + pz[2]; |
| 211 | Float_t sz = pS[0]*pz[0] + pS[1]*pz[1] + pS[2]*pz[2]; |
| 212 | Float_t ss = pS[0]*pS[0] + pS[1]*pS[1] + pS[2]*pS[2]; |
| 213 | |
| 214 | Float_t a = 3*ss-s*s; |
| 215 | Z() = (z*ss-sz*s)/a; // z0 |
| 216 | DzDs() = (3*sz-z*s)/a; // t = dz/ds |
| 217 | |
| 218 | Float_t dz0[3] = {ss - pS[0]*s,ss - pS[1]*s,ss - pS[2]*s }; |
| 219 | Float_t dt [3] = {3*pS[0] - s, 3*pS[1] - s, 3*pS[2] - s }; |
| 220 | |
| 221 | fC[2] = (dz0[0]*dz0[0]*ps2z[0] + dz0[1]*dz0[1]*ps2z[1] + dz0[2]*dz0[2]*ps2z[2])/a/a; |
| 222 | fC[7]= (dz0[0]*dt [0]*ps2z[0] + dz0[1]*dt [1]*ps2z[1] + dz0[2]*dt [2]*ps2z[2])/a/a; |
| 223 | fC[9]= (dt [0]*dt [0]*ps2z[0] + dt [1]*dt [1]*ps2z[1] + dt [2]*dt [2]*ps2z[2])/a/a; |
| 224 | } |
| 225 | |
| 226 | |
| 227 | Bool_t AliHLTTPCCATrackParam::TransportToX( Float_t x ) |
| 228 | { |
| 229 | //* Transport the track parameters to X=x |
| 230 | |
| 231 | Bool_t ret = 1; |
| 232 | |
| 233 | Float_t x0 = X(); |
| 234 | //Float_t y0 = Y(); |
| 235 | Float_t k = Kappa(); |
| 236 | Float_t ex = CosPhi(); |
| 237 | Float_t ey = SinPhi(); |
| 238 | Float_t dx = x - x0; |
| 239 | |
| 240 | Float_t ey1 = k*dx + ey; |
| 241 | Float_t ex1; |
| 242 | if( TMath::Abs(ey1)>1 ){ // no intersection -> check the border |
| 243 | ey1 = ( ey1>0 ) ?1 :-1; |
| 244 | ex1 = 0; |
| 245 | dx = ( TMath::Abs(k)>1.e-4) ? ( (ey1-ey)/k ) :0; |
| 246 | |
| 247 | Float_t ddx = TMath::Abs(x0+dx - x)*k*k; |
| 248 | Float_t hx[] = {0, -k, 1+ey }; |
| 249 | Float_t sx2 = hx[1]*hx[1]*fC[ 3] + hx[2]*hx[2]*fC[ 5]; |
| 250 | if( ddx*ddx>3.5*3.5*sx2 ) ret = 0; // x not withing the error |
| 251 | ret = 0; // any case |
| 252 | return ret; |
| 253 | }else{ |
| 254 | ex1 = TMath::Sqrt(1 - ey1*ey1); |
| 255 | if( ex<0 ) ex1 = -ex1; |
| 256 | } |
| 257 | |
| 258 | Float_t dx2 = dx*dx; |
| 259 | CosPhi() = ex1; |
| 260 | Float_t ss = ey+ey1; |
| 261 | Float_t cc = ex+ex1; |
| 262 | Float_t tg = 0; |
| 263 | if( TMath::Abs(cc)>1.e-4 ) tg = ss/cc; // tan((phi1+phi)/2) |
| 264 | else ret = 0; |
| 265 | Float_t dy = dx*tg; |
| 266 | Float_t dl = dx*TMath::Sqrt(1+tg*tg); |
| 267 | |
| 268 | if( cc<0 ) dl = -dl; |
| 269 | Float_t dSin = dl*k/2; |
| 270 | if( dSin > 1 ) dSin = 1; |
| 271 | if( dSin <-1 ) dSin = -1; |
| eb30eb49 |
272 | Float_t dS = ( TMath::Abs(k)>1.e-4) ? (2*TMath::ASin(dSin)/k) :dl; |
| d54804bf |
273 | Float_t dz = dS*DzDs(); |
| 274 | |
| 275 | Float_t cci = 0, exi = 0, ex1i = 0; |
| 276 | if( TMath::Abs(cc)>1.e-4 ) cci = 1./cc; |
| 277 | else ret = 0; |
| 278 | if( TMath::Abs(ex)>1.e-4 ) exi = 1./ex; |
| 279 | else ret = 0; |
| 280 | if( TMath::Abs(ex1)>1.e-4 ) ex1i = 1./ex1; |
| 281 | else ret = 0; |
| 282 | |
| 283 | if( !ret ) return ret; |
| eb30eb49 |
284 | |
| d54804bf |
285 | X() += dx; |
| 286 | fP[0]+= dy; |
| 287 | fP[1]+= dz; |
| 288 | fP[2] = ey1; |
| 289 | fP[3] = fP[3]; |
| 290 | fP[4] = fP[4]; |
| 291 | |
| 292 | Float_t h2 = dx*(1+ ex*ex1 + ey*ey1 )*cci*exi*ex1i; |
| 293 | Float_t h4 = dx2*(cc + ss*ey1*ex1i )*cci*cci; |
| 294 | |
| 295 | Float_t c00 = fC[0]; |
| 296 | Float_t c10 = fC[1]; |
| 297 | Float_t c11 = fC[2]; |
| 298 | Float_t c20 = fC[3]; |
| 299 | Float_t c21 = fC[4]; |
| 300 | Float_t c22 = fC[5]; |
| 301 | Float_t c30 = fC[6]; |
| 302 | Float_t c31 = fC[7]; |
| 303 | Float_t c32 = fC[8]; |
| 304 | Float_t c33 = fC[9]; |
| 305 | Float_t c40 = fC[10]; |
| 306 | Float_t c41 = fC[11]; |
| 307 | Float_t c42 = fC[12]; |
| 308 | Float_t c43 = fC[13]; |
| 309 | Float_t c44 = fC[14]; |
| 310 | |
| 311 | //Float_t H0[5] = { 1,0, h2, 0, h4 }; |
| 312 | //Float_t H1[5] = { 0, 1, 0, dS, 0 }; |
| 313 | //Float_t H2[5] = { 0, 0, 1, 0, dx }; |
| 314 | //Float_t H3[5] = { 0, 0, 0, 1, 0 }; |
| 315 | //Float_t H4[5] = { 0, 0, 0, 0, 1 }; |
| 316 | |
| 317 | |
| 318 | fC[0]=( c00 + h2*h2*c22 + h4*h4*c44 |
| 319 | + 2*( h2*c20 + h4*c40 + h2*h4*c42 ) ); |
| 320 | |
| 321 | fC[1]= c10 + h2*c21 + h4*c41 + dS*(c30 + h2*c32 + h4*c43); |
| 322 | fC[2]= c11 + 2*dS*c31 + dS*dS*c33; |
| 323 | |
| 324 | fC[3]= c20 + h2*c22 + h4*c42 + dx*( c40 + h2*c42 + h4*c44); |
| 325 | fC[4]= c21 + dS*c32 + dx*(c41 + dS*c43); |
| 326 | fC[5]= c22 +2*dx*c42 + dx2*c44; |
| 327 | |
| 328 | fC[6]= c30 + h2*c32 + h4*c43; |
| 329 | fC[7]= c31 + dS*c33; |
| 330 | fC[8]= c32 + dx*c43; |
| 331 | fC[9]= c33; |
| 332 | |
| 333 | fC[10]= c40 + h2*c42 + h4*c44; |
| 334 | fC[11]= c41 + dS*c43; |
| 335 | fC[12]= c42 + dx*c44; |
| 336 | fC[13]= c43; |
| 337 | fC[14]= c44; |
| 338 | |
| 339 | return ret; |
| 340 | } |
| 341 | |
| 0fad9753 |
342 | Bool_t AliHLTTPCCATrackParam::TransportToXWithMaterial( Float_t x, Float_t Bz ) |
| eb30eb49 |
343 | { |
| 344 | AliHLTTPCCATrackFitParam par; |
| 345 | CalculateFitParameters( par, Bz ); |
| 0fad9753 |
346 | return TransportToXWithMaterial(x, par ); |
| eb30eb49 |
347 | } |
| d54804bf |
348 | |
| 349 | |
| eb30eb49 |
350 | Bool_t AliHLTTPCCATrackParam::TransportToXWithMaterial( Float_t x, AliHLTTPCCATrackFitParam &par ) |
| d54804bf |
351 | { |
| eb30eb49 |
352 | //* Transport the track parameters to X=x |
| d54804bf |
353 | |
| 354 | Bool_t ret = 1; |
| 355 | |
| eb30eb49 |
356 | Float_t oldX=GetX(); |
| 357 | |
| 358 | Float_t x0 = X(); |
| 359 | //Float_t y0 = Y(); |
| 360 | Float_t k = Kappa(); |
| 361 | Float_t ex = CosPhi(); |
| 362 | Float_t ey = SinPhi(); |
| 363 | Float_t dx = x - x0; |
| 364 | |
| 365 | Float_t ey1 = k*dx + ey; |
| 366 | Float_t ex1; |
| 367 | if( TMath::Abs(ey1)>.99 ){ // no intersection -> check the border |
| 368 | ey1 = ( ey1>0 ) ?1 :-1; |
| 369 | ex1 = 0; |
| 370 | dx = ( TMath::Abs(k)>1.e-4) ? ( (ey1-ey)/k ) :0; |
| 371 | |
| 372 | Float_t ddx = TMath::Abs(x0+dx - x)*k*k; |
| 373 | Float_t hx[] = {0, -k, 1+ey }; |
| 374 | Float_t sx2 = hx[1]*hx[1]*fC[ 3] + hx[2]*hx[2]*fC[ 5]; |
| 375 | if( ddx*ddx>3.5*3.5*sx2 ) ret = 0; // x not withing the error |
| 376 | ret = 0; // any case |
| 377 | return ret; |
| 378 | }else{ |
| 379 | ex1 = TMath::Sqrt(1 - ey1*ey1); |
| 380 | if( ex<0 ) ex1 = -ex1; |
| 381 | } |
| 382 | |
| 383 | Float_t dx2 = dx*dx; |
| 384 | CosPhi() = ex1; |
| 385 | Float_t ss = ey+ey1; |
| 386 | Float_t cc = ex+ex1; |
| 387 | Float_t tg = 0; |
| 388 | if( TMath::Abs(cc)>1.e-4 ) tg = ss/cc; // tan((phi1+phi)/2) |
| 389 | else ret = 0; |
| 390 | Float_t dy = dx*tg; |
| 391 | Float_t dl = dx*TMath::Sqrt(1+tg*tg); |
| 392 | |
| 393 | if( cc<0 ) dl = -dl; |
| 394 | Float_t dSin = dl*k/2; |
| 395 | if( dSin > 1 ) dSin = 1; |
| 396 | if( dSin <-1 ) dSin = -1; |
| 397 | Float_t dS = ( TMath::Abs(k)>1.e-4) ? (2*TMath::ASin(dSin)/k) :dl; |
| 398 | Float_t dz = dS*DzDs(); |
| 399 | |
| 400 | Float_t cci = 0, exi = 0, ex1i = 0; |
| 401 | if( TMath::Abs(cc)>1.e-4 ) cci = 1./cc; |
| 402 | else ret = 0; |
| 403 | if( TMath::Abs(ex)>1.e-4 ) exi = 1./ex; |
| 404 | else ret = 0; |
| 405 | if( TMath::Abs(ex1)>1.e-4 ) ex1i = 1./ex1; |
| 406 | else ret = 0; |
| 407 | |
| 408 | if( !ret ) return ret; |
| 409 | |
| 410 | X() += dx; |
| 411 | fP[0]+= dy; |
| 412 | fP[1]+= dz; |
| 413 | fP[2] = ey1; |
| 414 | fP[3] = fP[3]; |
| 415 | fP[4] = fP[4]; |
| 416 | |
| 417 | Float_t h2 = dx*(1+ ex*ex1 + ey*ey1 )*cci*exi*ex1i; |
| 418 | Float_t h4 = dx2*(cc + ss*ey1*ex1i )*cci*cci; |
| 419 | |
| 420 | Float_t c00 = fC[0]; |
| 421 | Float_t c10 = fC[1]; |
| 422 | Float_t c11 = fC[2]; |
| 423 | Float_t c20 = fC[3]; |
| 424 | Float_t c21 = fC[4]; |
| 425 | Float_t c22 = fC[5]; |
| 426 | Float_t c30 = fC[6]; |
| 427 | Float_t c31 = fC[7]; |
| 428 | Float_t c32 = fC[8]; |
| 429 | Float_t c33 = fC[9]; |
| 430 | Float_t c40 = fC[10]; |
| 431 | Float_t c41 = fC[11]; |
| 432 | Float_t c42 = fC[12]; |
| 433 | Float_t c43 = fC[13]; |
| 434 | Float_t c44 = fC[14]; |
| 435 | |
| 436 | //Float_t H0[5] = { 1,0, h2, 0, h4 }; |
| 437 | //Float_t H1[5] = { 0, 1, 0, dS, 0 }; |
| 438 | //Float_t H2[5] = { 0, 0, 1, 0, dx }; |
| 439 | //Float_t H3[5] = { 0, 0, 0, 1, 0 }; |
| 440 | //Float_t H4[5] = { 0, 0, 0, 0, 1 }; |
| 441 | |
| 442 | |
| 443 | fC[0]=( c00 + h2*h2*c22 + h4*h4*c44 |
| 444 | + 2*( h2*c20 + h4*c40 + h2*h4*c42 ) ); |
| 445 | |
| 446 | fC[1]= c10 + h2*c21 + h4*c41 + dS*(c30 + h2*c32 + h4*c43); |
| 447 | fC[2]= c11 + 2*dS*c31 + dS*dS*c33; |
| 448 | |
| 449 | fC[3]= c20 + h2*c22 + h4*c42 + dx*( c40 + h2*c42 + h4*c44); |
| 450 | fC[4]= c21 + dS*c32 + dx*(c41 + dS*c43); |
| 451 | fC[5]= c22 +2*dx*c42 + dx2*c44; |
| 452 | |
| 453 | fC[6]= c30 + h2*c32 + h4*c43; |
| 454 | fC[7]= c31 + dS*c33; |
| 455 | fC[8]= c32 + dx*c43; |
| 456 | fC[9]= c33; |
| 457 | |
| 458 | fC[10]= c40 + h2*c42 + h4*c44; |
| 459 | fC[11]= c41 + dS*c43; |
| 460 | fC[12]= c42 + dx*c44; |
| 461 | fC[13]= c43; |
| 462 | fC[14]= c44; |
| 463 | |
| 464 | Float_t d = TMath::Sqrt(dS*dS + dz*dz ); |
| 465 | |
| 466 | if (oldX > GetX() ) d = -d; |
| 467 | { |
| 468 | Float_t rho=0.9e-3; |
| 469 | Float_t radLen=28.94; |
| 470 | CorrectForMeanMaterial(d*rho/radLen,d*rho,par); |
| 471 | } |
| 472 | |
| 473 | return ret; |
| 474 | } |
| 475 | |
| 476 | |
| 477 | |
| 478 | Float_t AliHLTTPCCATrackParam::ApproximateBetheBloch( Float_t beta2 ) |
| 479 | { |
| 480 | //------------------------------------------------------------------ |
| 481 | // This is an approximation of the Bethe-Bloch formula with |
| 482 | // the density effect taken into account at beta*gamma > 3.5 |
| 483 | // (the approximation is reasonable only for solid materials) |
| 484 | //------------------------------------------------------------------ |
| 485 | if (beta2 >= 1) return 0; |
| 486 | |
| 487 | if (beta2/(1-beta2)>3.5*3.5) |
| 488 | return 0.153e-3/beta2*( log(3.5*5940)+0.5*log(beta2/(1-beta2)) - beta2); |
| 489 | return 0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2); |
| 490 | } |
| 491 | |
| 492 | |
| 493 | void AliHLTTPCCATrackParam::CalculateFitParameters( AliHLTTPCCATrackFitParam &par, Float_t Bz, Float_t mass ) |
| 494 | { |
| 495 | //*! |
| 496 | |
| 497 | const Float_t kCLight = 0.000299792458; |
| 498 | Float_t c = Bz*kCLight; |
| 499 | Float_t p2 = (1.+ fP[3]*fP[3])*c*c; |
| 500 | Float_t k2 = fP[4]*fP[4]; |
| 501 | Float_t beta2= p2 / (p2 + mass*mass*k2); |
| 502 | Float_t Bethe = ApproximateBetheBloch(beta2); |
| 503 | |
| 504 | Float_t P2 = (k2>1.e-8) ?p2/k2 :10000; // impuls 2 |
| 505 | par.fBethe = Bethe; |
| 506 | par.fE = TMath::Sqrt( P2 + mass*mass); |
| 507 | par.fTheta2 = 14.1*14.1/(beta2*p2*1e6)*k2; |
| 508 | par.fEP2 = par.fE/p2*k2; |
| 509 | |
| 510 | // Approximate energy loss fluctuation (M.Ivanov) |
| 511 | |
| 512 | const Float_t knst=0.07; // To be tuned. |
| 513 | par.fSigmadE2 = knst*par.fEP2*fP[4]; |
| 514 | par.fSigmadE2 = par.fSigmadE2 * par.fSigmadE2; |
| 515 | |
| 516 | par.fK22 = (1. + fP[3]*fP[3]); |
| 517 | par.fK33 = par.fK22*par.fK22; |
| 518 | par.fK43 = fP[3]*fP[4]*par.fK22; |
| 519 | par.fK44 = fP[3]*fP[3]*fP[4]*fP[4]; |
| 520 | } |
| 521 | |
| 522 | |
| 523 | Bool_t AliHLTTPCCATrackParam::CorrectForMeanMaterial( Float_t xOverX0, Float_t xTimesRho, AliHLTTPCCATrackFitParam &par ) |
| 524 | { |
| 525 | //------------------------------------------------------------------ |
| 526 | // This function corrects the track parameters for the crossed material. |
| 527 | // "xOverX0" - X/X0, the thickness in units of the radiation length. |
| 528 | // "xTimesRho" - is the product length*density (g/cm^2). |
| 529 | //------------------------------------------------------------------ |
| 530 | |
| 531 | Float_t &fC22=fC[5]; |
| 532 | Float_t &fC33=fC[9]; |
| 533 | Float_t &fC40=fC[10]; |
| 534 | Float_t &fC41=fC[11]; |
| 535 | Float_t &fC42=fC[12]; |
| 536 | Float_t &fC43=fC[13]; |
| 537 | Float_t &fC44=fC[14]; |
| 538 | |
| 539 | //Energy losses************************ |
| 540 | |
| 541 | Float_t dE = par.fBethe*xTimesRho; |
| 542 | if ( TMath::Abs(dE) > 0.3*par.fE ) return 0; //30% energy loss is too much! |
| 543 | Float_t corr = (1.- par.fEP2*dE); |
| 544 | if( corr<0.3 ) return 0; |
| 545 | fP[4]*= corr; |
| 546 | fC40*= corr; |
| 547 | fC41*= corr; |
| 548 | fC42*= corr; |
| 549 | fC43*= corr; |
| 550 | fC44*= corr*corr; |
| 551 | fC44+= par.fSigmadE2*TMath::Abs(dE); |
| 552 | |
| 553 | |
| 554 | //Multiple scattering****************** |
| 555 | |
| 556 | Float_t theta2 = par.fTheta2*TMath::Abs(xOverX0); |
| 557 | fC22 += theta2*par.fK22*(1.- fP[2]*fP[2]); |
| 558 | fC33 += theta2*par.fK33; |
| 559 | fC43 += theta2*par.fK43; |
| 560 | fC44 += theta2*par.fK44; |
| 561 | |
| 562 | return 1; |
| 563 | } |
| 564 | |
| 565 | |
| 566 | |
| 567 | #include "Riostream.h" |
| 568 | |
| 569 | Bool_t AliHLTTPCCATrackParam::Rotate( Float_t alpha ) |
| 570 | { |
| 571 | //* Rotate the coordinate system in XY on the angle alpha |
| 572 | |
| d54804bf |
573 | Float_t cA = TMath::Cos( alpha ); |
| 574 | Float_t sA = TMath::Sin( alpha ); |
| 575 | Float_t x = X(), y= Y(), sP= SinPhi(), cP= CosPhi(); |
| eb30eb49 |
576 | Float_t cosPhi = cP*cA + sP*sA; |
| 577 | Float_t sinPhi =-cP*sA + sP*cA; |
| 578 | |
| 579 | if( TMath::Abs(sinPhi)>.99 || TMath::Abs(cosPhi)<1.e-2 || TMath::Abs(cP)<1.e-2 ) return 0; |
| 580 | |
| 581 | Float_t j0 = cP/cosPhi; |
| 582 | Float_t j2 = cosPhi/cP; |
| 583 | |
| d54804bf |
584 | X() = x*cA + y*sA; |
| 585 | Y() = -x*sA + y*cA; |
| eb30eb49 |
586 | CosPhi() = cosPhi; |
| 587 | SinPhi() = sinPhi; |
| d54804bf |
588 | |
| d54804bf |
589 | |
| 590 | //Float_t J[5][5] = { { j0, 0, 0, 0, 0 }, // Y |
| 591 | // { 0, 1, 0, 0, 0 }, // Z |
| 592 | // { 0, 0, j2, 0, 0 }, // SinPhi |
| 593 | // { 0, 0, 0, 1, 0 }, // DzDs |
| 594 | // { 0, 0, 0, 0, 1 } }; // Kappa |
| eb30eb49 |
595 | //cout<<"alpha="<<alpha<<" "<<x<<" "<<y<<" "<<sP<<" "<<cP<<" "<<j0<<" "<<j2<<endl; |
| 596 | //cout<<" "<<fC[0]<<" "<<fC[1]<<" "<<fC[6]<<" "<<fC[10]<<" "<<fC[4]<<" "<<fC[5]<<" "<<fC[8]<<" "<<fC[12]<<endl; |
| d54804bf |
597 | fC[0]*= j0*j0; |
| 598 | fC[1]*= j0; |
| 599 | //fC[3]*= j0; |
| 600 | fC[6]*= j0; |
| 601 | fC[10]*= j0; |
| 602 | |
| 603 | //fC[3]*= j2; |
| 604 | fC[4]*= j2; |
| 605 | fC[5]*= j2*j2; |
| 606 | fC[8]*= j2; |
| 607 | fC[12]*= j2; |
| eb30eb49 |
608 | //cout<<" "<<fC[0]<<" "<<fC[1]<<" "<<fC[6]<<" "<<fC[10]<<" "<<fC[4]<<" "<<fC[5]<<" "<<fC[8]<<" "<<fC[12]<<endl; |
| 609 | return 1; |
| d54804bf |
610 | } |
| 611 | |
| 612 | |
| eb30eb49 |
613 | Bool_t AliHLTTPCCATrackParam::Filter2( Float_t y, Float_t z, Float_t err2Y, Float_t err2Z ) |
| d54804bf |
614 | { |
| 615 | //* Add the y,z measurement with the Kalman filter |
| 616 | |
| 617 | Float_t |
| 618 | c00 = fC[ 0], |
| 619 | c10 = fC[ 1], c11 = fC[ 2], |
| 620 | c20 = fC[ 3], c21 = fC[ 4], |
| 621 | c30 = fC[ 6], c31 = fC[ 7], |
| 622 | c40 = fC[10], c41 = fC[11]; |
| 623 | |
| 624 | Float_t |
| 625 | z0 = y-fP[0], |
| 626 | z1 = z-fP[1]; |
| 627 | |
| eb30eb49 |
628 | Float_t v[3] = {err2Y, 0, err2Z}; |
| d54804bf |
629 | |
| 630 | Float_t mS[3] = { c00+v[0], c10+v[1], c11+v[2] }; |
| 631 | |
| 632 | Float_t mSi[3]; |
| 633 | Float_t det = (mS[0]*mS[2] - mS[1]*mS[1]); |
| 634 | |
| eb30eb49 |
635 | if( det < 1.e-8 ) return 0; |
| d54804bf |
636 | det = 1./det; |
| 637 | mSi[0] = mS[2]*det; |
| 638 | mSi[1] = -mS[1]*det; |
| 639 | mSi[2] = mS[0]*det; |
| 640 | |
| d54804bf |
641 | // K = CHtS |
| 642 | |
| 643 | Float_t k00, k01 , k10, k11, k20, k21, k30, k31, k40, k41; |
| 644 | |
| 645 | k00 = c00*mSi[0] + c10*mSi[1]; k01 = c00*mSi[1] + c10*mSi[2]; |
| 646 | k10 = c10*mSi[0] + c11*mSi[1]; k11 = c10*mSi[1] + c11*mSi[2]; |
| 647 | k20 = c20*mSi[0] + c21*mSi[1]; k21 = c20*mSi[1] + c21*mSi[2]; |
| 648 | k30 = c30*mSi[0] + c31*mSi[1]; k31 = c30*mSi[1] + c31*mSi[2] ; |
| 649 | k40 = c40*mSi[0] + c41*mSi[1]; k41 = c40*mSi[1] + c41*mSi[2] ; |
| 650 | |
| 651 | Float_t sinPhi = fP[2] + k20*z0 + k21*z1 ; |
| eb30eb49 |
652 | if( TMath::Abs(sinPhi)>=0.99 ) return 0; |
| 653 | |
| 654 | fNDF += 2; |
| 655 | fChi2 += ( +(mSi[0]*z0 + mSi[1]*z1 )*z0 |
| 656 | +(mSi[1]*z0 + mSi[2]*z1 )*z1 ); |
| d54804bf |
657 | |
| 658 | fP[ 0]+= k00*z0 + k01*z1 ; |
| 659 | fP[ 1]+= k10*z0 + k11*z1 ; |
| eb30eb49 |
660 | fP[ 2] = sinPhi; |
| d54804bf |
661 | fP[ 3]+= k30*z0 + k31*z1 ; |
| 662 | fP[ 4]+= k40*z0 + k41*z1 ; |
| 663 | |
| 664 | |
| 665 | fC[ 0]-= k00*c00 + k01*c10 ; |
| 666 | |
| 667 | fC[ 1]-= k10*c00 + k11*c10 ; |
| 668 | fC[ 2]-= k10*c10 + k11*c11 ; |
| 669 | |
| 670 | fC[ 3]-= k20*c00 + k21*c10 ; |
| 671 | fC[ 4]-= k20*c10 + k21*c11 ; |
| 672 | fC[ 5]-= k20*c20 + k21*c21 ; |
| 673 | |
| 674 | fC[ 6]-= k30*c00 + k31*c10 ; |
| 675 | fC[ 7]-= k30*c10 + k31*c11 ; |
| 676 | fC[ 8]-= k30*c20 + k31*c21 ; |
| 677 | fC[ 9]-= k30*c30 + k31*c31 ; |
| 678 | |
| 679 | fC[10]-= k40*c00 + k41*c10 ; |
| 680 | fC[11]-= k40*c10 + k41*c11 ; |
| 681 | fC[12]-= k40*c20 + k41*c21 ; |
| 682 | fC[13]-= k40*c30 + k41*c31 ; |
| 683 | fC[14]-= k40*c40 + k41*c41 ; |
| 684 | |
| 685 | if( CosPhi()>=0 ){ |
| 686 | CosPhi() = TMath::Sqrt(1-SinPhi()*SinPhi()); |
| 687 | }else{ |
| 688 | CosPhi() = -TMath::Sqrt(1-SinPhi()*SinPhi()); |
| 689 | } |
| eb30eb49 |
690 | return 1; |
| 691 | } |
| 692 | |
| 693 | void AliHLTTPCCATrackParam::FilterY( Float_t y, Float_t erry ) |
| 694 | { |
| 695 | //* Add the y measurement with the Kalman filter |
| 696 | |
| 697 | Float_t |
| 698 | c00 = fC[ 0], |
| 699 | c10 = fC[ 1], |
| 700 | c20 = fC[ 3], |
| 701 | c30 = fC[ 6], |
| 702 | c40 = fC[10]; |
| 703 | |
| 704 | Float_t |
| 705 | z0 = y-fP[0]; |
| 706 | |
| 707 | Float_t s = { c00+erry*erry }; |
| 708 | if( TMath::Abs(s)<1.e-4 ) return; |
| 709 | |
| 710 | Float_t si = 1/s; |
| 711 | |
| 712 | fNDF += 1; |
| 713 | fChi2 += si*z0*z0; |
| 714 | |
| 715 | // K = CHtS |
| 716 | |
| 717 | Float_t k0, k1 , k2, k3, k4; |
| d54804bf |
718 | |
| eb30eb49 |
719 | k0 = c00*si; |
| 720 | k1 = c10*si; |
| 721 | k2 = c20*si; |
| 722 | k3 = c30*si; |
| 723 | k4 = c40*si; |
| 724 | |
| 725 | Float_t sinPhi = fP[2] + k2*z0 ; |
| 726 | if( TMath::Abs(sinPhi)>=0.99 ) return; |
| 727 | |
| 728 | fP[ 0]+= k0*z0 ; |
| 729 | fP[ 1]+= k1*z0 ; |
| 730 | fP[ 2] = sinPhi; |
| 731 | fP[ 3]+= k3*z0 ; |
| 732 | fP[ 4]+= k4*z0 ; |
| 733 | |
| 734 | fC[ 0]-= k0*c00; |
| 735 | |
| 736 | fC[ 1]-= k1*c00; |
| 737 | fC[ 2]-= k1*c10; |
| 738 | |
| 739 | fC[ 3]-= k2*c00; |
| 740 | fC[ 4]-= k2*c10; |
| 741 | fC[ 5]-= k2*c20; |
| 742 | |
| 743 | fC[ 6]-= k3*c00; |
| 744 | fC[ 7]-= k3*c10; |
| 745 | fC[ 8]-= k3*c20; |
| 746 | fC[ 9]-= k3*c30; |
| 747 | |
| 748 | fC[10]-= k4*c00; |
| 749 | fC[11]-= k4*c10; |
| 750 | fC[12]-= k4*c20; |
| 751 | fC[13]-= k4*c30; |
| 752 | fC[14]-= k4*c40; |
| 753 | |
| 754 | if( CosPhi()>=0 ){ |
| 755 | CosPhi() = TMath::Sqrt(1-SinPhi()*SinPhi()); |
| 756 | }else{ |
| 757 | CosPhi() = -TMath::Sqrt(1-SinPhi()*SinPhi()); |
| 758 | } |
| 759 | |
| 760 | } |
| 761 | |
| 762 | void AliHLTTPCCATrackParam::FilterZ( Float_t z, Float_t errz ) |
| 763 | { |
| 764 | //* Add the z measurement with the Kalman filter |
| 765 | |
| 766 | Float_t |
| 767 | c01 = fC[ 1], |
| 768 | c11 = fC[ 2], |
| 769 | c21 = fC[ 4], |
| 770 | c31 = fC[ 7], |
| 771 | c41 = fC[11]; |
| 772 | |
| 773 | Float_t |
| 774 | z1 = z-fP[1]; |
| 775 | |
| 776 | Float_t s = c11 + errz*errz; |
| 777 | if( TMath::Abs(s)<1.e-4 ) return; |
| 778 | |
| 779 | Float_t si = 1./s; |
| 780 | |
| 781 | fNDF += 1; |
| 782 | fChi2 += si*z1*z1; |
| 783 | |
| 784 | // K = CHtS |
| 785 | |
| 786 | Float_t k0, k1 , k2, k3, k4; |
| 787 | |
| 788 | k0 = 0;//c01*si; |
| 789 | k1 = c11*si; |
| 790 | k2 = 0;//c21*si; |
| 791 | k3 = c31*si; |
| 792 | k4 = 0;//c41*si; |
| 793 | |
| 794 | Float_t sinPhi = fP[2] + k2*z1 ; |
| 795 | if( TMath::Abs(sinPhi)>=0.99 ) return; |
| 796 | |
| 797 | fP[ 0]+= k0*z1 ; |
| 798 | fP[ 1]+= k1*z1 ; |
| 799 | fP[ 2] = sinPhi; |
| 800 | fP[ 3]+= k3*z1 ; |
| 801 | fP[ 4]+= k4*z1 ; |
| 802 | |
| 803 | |
| 804 | fC[ 0]-= k0*c01 ; |
| 805 | |
| 806 | fC[ 1]-= k1*c01 ; |
| 807 | fC[ 2]-= k1*c11 ; |
| 808 | |
| 809 | fC[ 3]-= k2*c01 ; |
| 810 | fC[ 4]-= k2*c11 ; |
| 811 | fC[ 5]-= k2*c21 ; |
| 812 | |
| 813 | fC[ 6]-= k3*c01 ; |
| 814 | fC[ 7]-= k3*c11 ; |
| 815 | fC[ 8]-= k3*c21 ; |
| 816 | fC[ 9]-= k3*c31 ; |
| 817 | |
| 818 | fC[10]-= k4*c01 ; |
| 819 | fC[11]-= k4*c11 ; |
| 820 | fC[12]-= k4*c21 ; |
| 821 | fC[13]-= k4*c31 ; |
| 822 | fC[14]-= k4*c41 ; |
| 823 | |
| 824 | if( CosPhi()>=0 ){ |
| 825 | CosPhi() = TMath::Sqrt(1-SinPhi()*SinPhi()); |
| 826 | }else{ |
| 827 | CosPhi() = -TMath::Sqrt(1-SinPhi()*SinPhi()); |
| 828 | } |
| 829 | |
| 830 | } |
| 831 | |
| 832 | void AliHLTTPCCATrackParam::Print() const |
| 833 | { |
| 834 | cout<<"track: "<<GetX()<<" "<<GetY()<<" "<<GetZ()<<" "<<GetSinPhi()<<" "<<GetDzDs()<<" "<<GetKappa()<<endl; |
| 835 | cout<<"errs2: "<<GetErr2Y()<<" "<<GetErr2Z()<<" "<<GetErr2SinPhi()<<" "<<GetErr2DzDs()<<" "<<GetErr2Kappa()<<endl; |
| d54804bf |
836 | } |
git://git.uio.no / u / mrichter / AliRoot.git / blame