Fit mathematics improved, obsollete GBTracker cleaned up
[u/mrichter/AliRoot.git] / HLT / TPCLib / tracking-ca / AliHLTTPCCATrackParam.cxx
CommitLineData
d54804bf 1// $Id$
ce565086 2// **************************************************************************
d54804bf 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. *
ce565086 17// *
d54804bf 18//***************************************************************************
19
ce565086 20
d54804bf 21#include "AliHLTTPCCATrackParam.h"
00d07bcd 22#include "AliHLTTPCCAMath.h"
15d2e9cf 23#include "AliHLTTPCCATrackLinearisation.h"
693d2443 24#include <iostream>
d54804bf 25
26//
27// Circle in XY:
15d2e9cf 28//
29// kCLight = 0.000299792458;
30// Kappa = Bz*kCLight*QPt;
d54804bf 31// R = 1/TMath::Abs(Kappa);
32// Xc = X - sin(Phi)/Kappa;
33// Yc = Y + cos(Phi)/Kappa;
34//
35
693d2443 36GPUd() Float_t AliHLTTPCCATrackParam::GetDist2( const AliHLTTPCCATrackParam &t ) const
37{
38 // get squared distance between tracks
39
40 Float_t dx = GetX() - t.GetX();
41 Float_t dy = GetY() - t.GetY();
42 Float_t dz = GetZ() - t.GetZ();
43 return dx*dx + dy*dy + dz*dz;
44}
45
46GPUd() Float_t AliHLTTPCCATrackParam::GetDistXZ2( const AliHLTTPCCATrackParam &t ) const
47{
48 // get squared distance between tracks in X&Z
49
50 Float_t dx = GetX() - t.GetX();
51 Float_t dz = GetZ() - t.GetZ();
52 return dx*dx + dz*dz;
53}
d54804bf 54
55
15d2e9cf 56GPUd() Float_t AliHLTTPCCATrackParam::GetS( Float_t x, Float_t y, Float_t Bz ) const
d54804bf 57{
58 //* Get XY path length to the given point
59
15d2e9cf 60 Float_t k = GetKappa( Bz );
d54804bf 61 Float_t ex = GetCosPhi();
62 Float_t ey = GetSinPhi();
63 x-= GetX();
64 y-= GetY();
65 Float_t dS = x*ex + y*ey;
00d07bcd 66 if( CAMath::Abs(k)>1.e-4 ) dS = CAMath::ATan2( k*dS, 1+k*(x*ey-y*ex) )/k;
d54804bf 67 return dS;
68}
69
00d07bcd 70GPUd() void AliHLTTPCCATrackParam::GetDCAPoint( Float_t x, Float_t y, Float_t z,
15d2e9cf 71 Float_t &xp, Float_t &yp, Float_t &zp,
72 Float_t Bz ) const
d54804bf 73{
74 //* Get the track point closest to the (x,y,z)
75
76 Float_t x0 = GetX();
77 Float_t y0 = GetY();
15d2e9cf 78 Float_t k = GetKappa( Bz );
d54804bf 79 Float_t ex = GetCosPhi();
80 Float_t ey = GetSinPhi();
81 Float_t dx = x - x0;
82 Float_t dy = y - y0;
83 Float_t ax = dx*k+ey;
84 Float_t ay = dy*k-ex;
85 Float_t a = sqrt( ax*ax+ay*ay );
86 xp = x0 + (dx - ey*( (dx*dx+dy*dy)*k - 2*(-dx*ey+dy*ex) )/(a+1) )/a;
87 yp = y0 + (dy + ex*( (dx*dx+dy*dy)*k - 2*(-dx*ey+dy*ex) )/(a+1) )/a;
15d2e9cf 88 Float_t s = GetS(x,y, Bz);
d54804bf 89 zp = GetZ() + GetDzDs()*s;
00d07bcd 90 if( CAMath::Abs(k)>1.e-2 ){
91 Float_t dZ = CAMath::Abs( GetDzDs()*CAMath::TwoPi()/k );
d54804bf 92 if( dZ>.1 ){
00d07bcd 93 zp+= CAMath::Nint((z-zp)/dZ)*dZ;
d54804bf 94 }
95 }
96}
97
d54804bf 98
15d2e9cf 99//*
100//* Transport routines
101//*
d54804bf 102
103
15d2e9cf 104GPUd() Bool_t AliHLTTPCCATrackParam::TransportToX( Float_t x, AliHLTTPCCATrackLinearisation &t0, Float_t Bz, Float_t maxSinPhi, Float_t *DL )
d54804bf 105{
15d2e9cf 106 //* Transport the track parameters to X=x, using linearization at t0, and the field value Bz
107 //* maxSinPhi is the max. allowed value for |t0.SinPhi()|
108 //* linearisation of trajectory t0 is also transported to X=x,
109 //* returns 1 if OK
110 //*
d54804bf 111
15d2e9cf 112 Float_t ex = t0.CosPhi();
113 Float_t ey = t0.SinPhi();
114 Float_t k = t0.QPt()*Bz;
115 Float_t dx = x - X();
d54804bf 116
117 Float_t ey1 = k*dx + ey;
118 Float_t ex1;
15d2e9cf 119
120 // check for intersection with X=x
121
122 if( CAMath::Abs(ey1)>maxSinPhi ) return 0;
123
124 ex1 = CAMath::Sqrt(1 - ey1*ey1);
125 if( ex<0 ) ex1 = -ex1;
d54804bf 126
127 Float_t dx2 = dx*dx;
d54804bf 128 Float_t ss = ey+ey1;
129 Float_t cc = ex+ex1;
00d07bcd 130
131 if( CAMath::Abs(cc)<1.e-4 || CAMath::Abs(ex)<1.e-4 || CAMath::Abs(ex1)<1.e-4 ) return 0;
132
133 Float_t tg = ss/cc; // tan((phi1+phi)/2)
134
d54804bf 135 Float_t dy = dx*tg;
00d07bcd 136 Float_t dl = dx*CAMath::Sqrt(1+tg*tg);
d54804bf 137
138 if( cc<0 ) dl = -dl;
139 Float_t dSin = dl*k/2;
140 if( dSin > 1 ) dSin = 1;
141 if( dSin <-1 ) dSin = -1;
00d07bcd 142 Float_t dS = ( CAMath::Abs(k)>1.e-4) ? (2*CAMath::ASin(dSin)/k) :dl;
15d2e9cf 143 Float_t dz = dS*t0.DzDs();
144
145 if( DL ) *DL = -dS*CAMath::Sqrt(1 + t0.DzDs()*t0.DzDs() );
d54804bf 146
00d07bcd 147 Float_t cci = 1./cc;
148 Float_t exi = 1./ex;
149 Float_t ex1i = 1./ex1;
150
15d2e9cf 151 Float_t d[5] = { 0,
152 0,
153 fP[2]-t0.SinPhi(),
154 fP[3]-t0.DzDs(),
155 fP[4]-t0.QPt() };
d54804bf 156
157 //Float_t H0[5] = { 1,0, h2, 0, h4 };
158 //Float_t H1[5] = { 0, 1, 0, dS, 0 };
15d2e9cf 159 //Float_t H2[5] = { 0, 0, 1, 0, dxBz };
d54804bf 160 //Float_t H3[5] = { 0, 0, 0, 1, 0 };
161 //Float_t H4[5] = { 0, 0, 0, 0, 1 };
693d2443 162
15d2e9cf 163 Float_t h2 = dx*(1+ey*ey1 + ex*ex1)*exi*ex1i*cci;
164 Float_t h4 = dx2*(cc + ss*ey1*ex1i )*cci*cci*Bz;
165 Float_t dxBz = dx*Bz;
693d2443 166
15d2e9cf 167 t0.SetCosPhi( ex1 );
168 t0.SetSinPhi( ey1 );
693d2443 169
15d2e9cf 170 fX = X() + dx;
171 fP[0] = Y() + dy + h2*d[2] + h4*d[4];
172 fP[1] = Z() + dz + dS*d[3];
173 fP[2] = t0.SinPhi() + d[2] + dxBz*d[4];
4687b8fc 174
693d2443 175 Float_t c00 = fC[0];
176 Float_t c10 = fC[1];
177 Float_t c11 = fC[2];
178 Float_t c20 = fC[3];
179 Float_t c21 = fC[4];
180 Float_t c22 = fC[5];
181 Float_t c30 = fC[6];
182 Float_t c31 = fC[7];
183 Float_t c32 = fC[8];
184 Float_t c33 = fC[9];
185 Float_t c40 = fC[10];
186 Float_t c41 = fC[11];
187 Float_t c42 = fC[12];
188 Float_t c43 = fC[13];
189 Float_t c44 = fC[14];
190
693d2443 191 fC[0]=( c00 + h2*h2*c22 + h4*h4*c44
192 + 2*( h2*c20 + h4*c40 + h2*h4*c42 ) );
193
194 fC[1]= c10 + h2*c21 + h4*c41 + dS*(c30 + h2*c32 + h4*c43);
195 fC[2]= c11 + 2*dS*c31 + dS*dS*c33;
196
15d2e9cf 197 fC[3]= c20 + h2*c22 + h4*c42 + dxBz*( c40 + h2*c42 + h4*c44);
198 fC[4]= c21 + dS*c32 + dxBz*(c41 + dS*c43);
199 fC[5]= c22 +2*dxBz*c42 + dxBz*dxBz*c44;
693d2443 200
201 fC[6]= c30 + h2*c32 + h4*c43;
202 fC[7]= c31 + dS*c33;
15d2e9cf 203 fC[8]= c32 + dxBz*c43;
693d2443 204 fC[9]= c33;
205
206 fC[10]= c40 + h2*c42 + h4*c44;
207 fC[11]= c41 + dS*c43;
15d2e9cf 208 fC[12]= c42 + dxBz*c44;
693d2443 209 fC[13]= c43;
210 fC[14]= c44;
211
212 return 1;
213}
214
215
15d2e9cf 216GPUd() Bool_t AliHLTTPCCATrackParam::TransportToX( Float_t x, Float_t sinPhi0, Float_t cosPhi0, Float_t Bz, Float_t maxSinPhi )
693d2443 217{
15d2e9cf 218 //* Transport the track parameters to X=x, using linearization at phi0 with 0 curvature,
219 //* and the field value Bz
220 //* maxSinPhi is the max. allowed value for |t0.SinPhi()|
221 //* linearisation of trajectory t0 is also transported to X=x,
222 //* returns 1 if OK
223 //*
224
225 Float_t ex = cosPhi0;
226 Float_t ey = sinPhi0;
227 Float_t dx = x - X();
693d2443 228
229 if( CAMath::Abs(ex)<1.e-4 ) return 0;
230 Float_t exi = 1./ex;
4687b8fc 231
15d2e9cf 232 Float_t dxBz = dx*Bz;
233 Float_t dS = dx*exi;
234 Float_t h2 = dS*exi*exi;
235 Float_t h4 = .5*h2*dxBz;
236
693d2443 237 //Float_t H0[5] = { 1,0, h2, 0, h4 };
238 //Float_t H1[5] = { 0, 1, 0, dS, 0 };
15d2e9cf 239 //Float_t H2[5] = { 0, 0, 1, 0, dxBz };
693d2443 240 //Float_t H3[5] = { 0, 0, 0, 1, 0 };
241 //Float_t H4[5] = { 0, 0, 0, 0, 1 };
242
15d2e9cf 243 Float_t sinPhi = SinPhi() + dxBz*QPt();
244 if( maxSinPhi>0 && CAMath::Abs(sinPhi)>maxSinPhi ) return 0;
693d2443 245
15d2e9cf 246 fX = X() + dx;
247 fP[0]+= dS*ey + h2*( SinPhi() - ey ) + h4*QPt();
248 fP[1]+= dS*DzDs();
249 fP[2] = sinPhi;
693d2443 250
693d2443 251
693d2443 252 Float_t c00 = fC[0];
253 Float_t c10 = fC[1];
254 Float_t c11 = fC[2];
255 Float_t c20 = fC[3];
256 Float_t c21 = fC[4];
257 Float_t c22 = fC[5];
258 Float_t c30 = fC[6];
259 Float_t c31 = fC[7];
260 Float_t c32 = fC[8];
261 Float_t c33 = fC[9];
262 Float_t c40 = fC[10];
263 Float_t c41 = fC[11];
264 Float_t c42 = fC[12];
265 Float_t c43 = fC[13];
266 Float_t c44 = fC[14];
267
693d2443 268
269 fC[0]=( c00 + h2*h2*c22 + h4*h4*c44
270 + 2*( h2*c20 + h4*c40 + h2*h4*c42 ) );
4687b8fc 271
693d2443 272 fC[1]= c10 + h2*c21 + h4*c41 + dS*(c30 + h2*c32 + h4*c43);
273 fC[2]= c11 + 2*dS*c31 + dS*dS*c33;
274
15d2e9cf 275 fC[3]= c20 + h2*c22 + h4*c42 + dxBz*( c40 + h2*c42 + h4*c44);
276 fC[4]= c21 + dS*c32 + dxBz*(c41 + dS*c43);
277 fC[5]= c22 +2*dxBz*c42 + dxBz*dxBz*c44;
693d2443 278
279 fC[6]= c30 + h2*c32 + h4*c43;
280 fC[7]= c31 + dS*c33;
15d2e9cf 281 fC[8]= c32 + dxBz*c43;
693d2443 282 fC[9]= c33;
283
284 fC[10]= c40 + h2*c42 + h4*c44;
285 fC[11]= c41 + dS*c43;
15d2e9cf 286 fC[12]= c42 + dxBz*c44;
693d2443 287 fC[13]= c43;
288 fC[14]= c44;
289
290 return 1;
291}
4687b8fc 292
293
d54804bf 294
295
d54804bf 296
d54804bf 297
15d2e9cf 298GPUd() Bool_t AliHLTTPCCATrackParam::TransportToX( Float_t x, Float_t Bz, Float_t maxSinPhi )
299{
300 //* Transport the track parameters to X=x
eb30eb49 301
15d2e9cf 302 AliHLTTPCCATrackLinearisation t0(*this);
eb30eb49 303
15d2e9cf 304 return TransportToX( x, t0, Bz, maxSinPhi );
305}
eb30eb49 306
eb30eb49 307
eb30eb49 308
15d2e9cf 309GPUd() Bool_t AliHLTTPCCATrackParam::TransportToXWithMaterial( Float_t x, AliHLTTPCCATrackLinearisation &t0, AliHLTTPCCATrackFitParam &par, Float_t Bz, Float_t maxSinPhi )
310{
311 //* Transport the track parameters to X=x taking into account material budget
eb30eb49 312
15d2e9cf 313 const Float_t kRho = 1.025e-3;//0.9e-3;
314 const Float_t kRadLen=29.532;//28.94;
315 const Float_t kRhoOverRadLen = kRho/kRadLen;
316 Float_t dl;
eb30eb49 317
15d2e9cf 318 if( !TransportToX( x, t0, Bz, maxSinPhi, &dl ) ) return 0;
eb30eb49 319
15d2e9cf 320 CorrectForMeanMaterial(dl*kRhoOverRadLen,dl*kRho,par);
321 return 1;
322}
eb30eb49 323
eb30eb49 324
15d2e9cf 325GPUd() Bool_t AliHLTTPCCATrackParam::TransportToXWithMaterial( Float_t x, AliHLTTPCCATrackFitParam &par, Float_t Bz, Float_t maxSinPhi )
326{
327 //* Transport the track parameters to X=x taking into account material budget
eb30eb49 328
15d2e9cf 329 AliHLTTPCCATrackLinearisation t0(*this);
330 return TransportToXWithMaterial( x, t0, par, Bz, maxSinPhi );
331}
eb30eb49 332
15d2e9cf 333GPUd() Bool_t AliHLTTPCCATrackParam::TransportToXWithMaterial( Float_t x, Float_t Bz, Float_t maxSinPhi )
334{
335 //* Transport the track parameters to X=x taking into account material budget
eb30eb49 336
15d2e9cf 337 AliHLTTPCCATrackFitParam par;
338 CalculateFitParameters( par );
339 return TransportToXWithMaterial(x, par, Bz, maxSinPhi );
340}
eb30eb49 341
eb30eb49 342
15d2e9cf 343//*
344//* Multiple scattering and energy losses
345//*
eb30eb49 346
eb30eb49 347
15d2e9cf 348Float_t AliHLTTPCCATrackParam::BetheBlochGeant(Float_t bg2,
349 Float_t kp0,
350 Float_t kp1,
351 Float_t kp2,
352 Float_t kp3,
353 Float_t kp4)
354{
355 //
356 // This is the parameterization of the Bethe-Bloch formula inspired by Geant.
357 //
358 // bg2 - (beta*gamma)^2
359 // kp0 - density [g/cm^3]
360 // kp1 - density effect first junction point
361 // kp2 - density effect second junction point
362 // kp3 - mean excitation energy [GeV]
363 // kp4 - mean Z/A
364 //
365 // The default values for the kp* parameters are for silicon.
366 // The returned value is in [GeV/(g/cm^2)].
367 //
368
369 const Float_t mK = 0.307075e-3; // [GeV*cm^2/g]
370 const Float_t me = 0.511e-3; // [GeV/c^2]
371 const Float_t rho = kp0;
372 const Float_t x0 = kp1*2.303;
373 const Float_t x1 = kp2*2.303;
374 const Float_t mI = kp3;
375 const Float_t mZA = kp4;
376 const Float_t maxT= 2*me*bg2; // neglecting the electron mass
377
378 //*** Density effect
379 Float_t d2=0.;
380 const Float_t x=0.5*TMath::Log(bg2);
381 const Float_t lhwI=TMath::Log(28.816*1e-9*TMath::Sqrt(rho*mZA)/mI);
382 if (x > x1) {
383 d2 = lhwI + x - 0.5;
384 } else if (x > x0) {
385 const Float_t r=(x1-x)/(x1-x0);
386 d2 = lhwI + x - 0.5 + (0.5 - lhwI - x0)*r*r*r;
eb30eb49 387 }
388
15d2e9cf 389 return mK*mZA*(1+bg2)/bg2*(0.5*TMath::Log(2*me*bg2*maxT/(mI*mI)) - bg2/(1+bg2) - d2);
eb30eb49 390}
391
15d2e9cf 392Float_t AliHLTTPCCATrackParam::BetheBlochSolid(Float_t bg)
693d2443 393{
15d2e9cf 394 //------------------------------------------------------------------
395 // This is an approximation of the Bethe-Bloch formula,
396 // reasonable for solid materials.
397 // All the parameters are, in fact, for Si.
398 // The returned value is in [GeV]
399 //------------------------------------------------------------------
693d2443 400
15d2e9cf 401 return BetheBlochGeant(bg);
402}
693d2443 403
15d2e9cf 404Float_t AliHLTTPCCATrackParam::BetheBlochGas(Float_t bg)
405{
406 //------------------------------------------------------------------
407 // This is an approximation of the Bethe-Bloch formula,
408 // reasonable for gas materials.
409 // All the parameters are, in fact, for Ne.
410 // The returned value is in [GeV]
411 //------------------------------------------------------------------
693d2443 412
15d2e9cf 413 const Float_t rho = 0.9e-3;
414 const Float_t x0 = 2.;
415 const Float_t x1 = 4.;
416 const Float_t mI = 140.e-9;
417 const Float_t mZA = 0.49555;
693d2443 418
15d2e9cf 419 return BetheBlochGeant(bg,rho,x0,x1,mI,mZA);
420}
693d2443 421
693d2443 422
693d2443 423
424
00d07bcd 425GPUd() Float_t AliHLTTPCCATrackParam::ApproximateBetheBloch( Float_t beta2 )
eb30eb49 426{
427 //------------------------------------------------------------------
428 // This is an approximation of the Bethe-Bloch formula with
429 // the density effect taken into account at beta*gamma > 3.5
430 // (the approximation is reasonable only for solid materials)
431 //------------------------------------------------------------------
432 if (beta2 >= 1) return 0;
433
434 if (beta2/(1-beta2)>3.5*3.5)
435 return 0.153e-3/beta2*( log(3.5*5940)+0.5*log(beta2/(1-beta2)) - beta2);
436 return 0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2);
437}
438
439
15d2e9cf 440GPUd() void AliHLTTPCCATrackParam::CalculateFitParameters( AliHLTTPCCATrackFitParam &par, Float_t mass )
eb30eb49 441{
442 //*!
443
15d2e9cf 444 Float_t p2 = (1.+ fP[3]*fP[3]);
eb30eb49 445 Float_t k2 = fP[4]*fP[4];
15d2e9cf 446 Float_t mass2 = mass*mass;
447 Float_t beta2= p2 / (p2 + mass2*k2);
eb30eb49 448
15d2e9cf 449 Float_t pp2 = (k2>1.e-8) ?p2/k2 :10000; // impuls 2
450
451 //par.fBethe = BetheBlochGas( pp2/mass2);
452 par.fBethe = ApproximateBetheBloch( pp2/mass2);
453 par.fE = CAMath::Sqrt( pp2 + mass2);
454 par.fTheta2 = 14.1*14.1/(beta2*pp2*1e6);
455 par.fEP2 = par.fE/pp2;
eb30eb49 456
457 // Approximate energy loss fluctuation (M.Ivanov)
458
459 const Float_t knst=0.07; // To be tuned.
15d2e9cf 460 par.fSigmadE2 = knst*par.fEP2*fP[4];
eb30eb49 461 par.fSigmadE2 = par.fSigmadE2 * par.fSigmadE2;
15d2e9cf 462
eb30eb49 463 par.fK22 = (1. + fP[3]*fP[3]);
464 par.fK33 = par.fK22*par.fK22;
465 par.fK43 = fP[3]*fP[4]*par.fK22;
466 par.fK44 = fP[3]*fP[3]*fP[4]*fP[4];
15d2e9cf 467
eb30eb49 468}
469
470
ce565086 471GPUd() Bool_t AliHLTTPCCATrackParam::CorrectForMeanMaterial( Float_t xOverX0, Float_t xTimesRho, const AliHLTTPCCATrackFitParam &par )
eb30eb49 472{
473 //------------------------------------------------------------------
474 // This function corrects the track parameters for the crossed material.
475 // "xOverX0" - X/X0, the thickness in units of the radiation length.
476 // "xTimesRho" - is the product length*density (g/cm^2).
477 //------------------------------------------------------------------
478
479 Float_t &fC22=fC[5];
480 Float_t &fC33=fC[9];
481 Float_t &fC40=fC[10];
482 Float_t &fC41=fC[11];
483 Float_t &fC42=fC[12];
484 Float_t &fC43=fC[13];
485 Float_t &fC44=fC[14];
486
487 //Energy losses************************
15d2e9cf 488
eb30eb49 489 Float_t dE = par.fBethe*xTimesRho;
00d07bcd 490 if ( CAMath::Abs(dE) > 0.3*par.fE ) return 0; //30% energy loss is too much!
eb30eb49 491 Float_t corr = (1.- par.fEP2*dE);
15d2e9cf 492 if( corr<0.3 || corr>1.3 ) return 0;
493
eb30eb49 494 fP[4]*= corr;
495 fC40*= corr;
496 fC41*= corr;
497 fC42*= corr;
498 fC43*= corr;
499 fC44*= corr*corr;
00d07bcd 500 fC44+= par.fSigmadE2*CAMath::Abs(dE);
eb30eb49 501
502 //Multiple scattering******************
503
00d07bcd 504 Float_t theta2 = par.fTheta2*CAMath::Abs(xOverX0);
eb30eb49 505 fC22 += theta2*par.fK22*(1.- fP[2]*fP[2]);
506 fC33 += theta2*par.fK33;
693d2443 507 fC43 += theta2*par.fK43;
eb30eb49 508 fC44 += theta2*par.fK44;
15d2e9cf 509
eb30eb49 510 return 1;
511}
512
513
15d2e9cf 514//*
515//* Rotation
516//*
517
eb30eb49 518
693d2443 519GPUd() Bool_t AliHLTTPCCATrackParam::Rotate( Float_t alpha, Float_t maxSinPhi )
eb30eb49 520{
521 //* Rotate the coordinate system in XY on the angle alpha
522
00d07bcd 523 Float_t cA = CAMath::Cos( alpha );
524 Float_t sA = CAMath::Sin( alpha );
15d2e9cf 525 Float_t x = X(), y= Y(), sP= SinPhi(), cP= GetCosPhi();
eb30eb49 526 Float_t cosPhi = cP*cA + sP*sA;
527 Float_t sinPhi =-cP*sA + sP*cA;
528
693d2443 529 if( CAMath::Abs(sinPhi)> maxSinPhi || CAMath::Abs(cosPhi)<1.e-2 || CAMath::Abs(cP)<1.e-2 ) return 0;
eb30eb49 530
531 Float_t j0 = cP/cosPhi;
532 Float_t j2 = cosPhi/cP;
533
693d2443 534 SetX( x*cA + y*sA );
535 SetY(-x*sA + y*cA );
15d2e9cf 536 SetSignCosPhi( cosPhi );
693d2443 537 SetSinPhi( sinPhi );
d54804bf 538
d54804bf 539
540 //Float_t J[5][5] = { { j0, 0, 0, 0, 0 }, // Y
541 // { 0, 1, 0, 0, 0 }, // Z
542 // { 0, 0, j2, 0, 0 }, // SinPhi
543 // { 0, 0, 0, 1, 0 }, // DzDs
544 // { 0, 0, 0, 0, 1 } }; // Kappa
eb30eb49 545 //cout<<"alpha="<<alpha<<" "<<x<<" "<<y<<" "<<sP<<" "<<cP<<" "<<j0<<" "<<j2<<endl;
546 //cout<<" "<<fC[0]<<" "<<fC[1]<<" "<<fC[6]<<" "<<fC[10]<<" "<<fC[4]<<" "<<fC[5]<<" "<<fC[8]<<" "<<fC[12]<<endl;
d54804bf 547 fC[0]*= j0*j0;
548 fC[1]*= j0;
693d2443 549 fC[3]*= j0;
d54804bf 550 fC[6]*= j0;
551 fC[10]*= j0;
552
693d2443 553 fC[3]*= j2;
d54804bf 554 fC[4]*= j2;
555 fC[5]*= j2*j2;
556 fC[8]*= j2;
557 fC[12]*= j2;
eb30eb49 558 //cout<<" "<<fC[0]<<" "<<fC[1]<<" "<<fC[6]<<" "<<fC[10]<<" "<<fC[4]<<" "<<fC[5]<<" "<<fC[8]<<" "<<fC[12]<<endl;
559 return 1;
d54804bf 560}
561
15d2e9cf 562GPUd() Bool_t AliHLTTPCCATrackParam::Rotate( Float_t alpha, AliHLTTPCCATrackLinearisation &t0, Float_t maxSinPhi )
693d2443 563{
564 //* Rotate the coordinate system in XY on the angle alpha
565
566 Float_t cA = CAMath::Cos( alpha );
567 Float_t sA = CAMath::Sin( alpha );
15d2e9cf 568 Float_t x0 = X(), y0= Y(), sP= t0.SinPhi(), cP= t0.CosPhi();
693d2443 569 Float_t cosPhi = cP*cA + sP*sA;
570 Float_t sinPhi =-cP*sA + sP*cA;
571
572 if( CAMath::Abs(sinPhi)> maxSinPhi || CAMath::Abs(cosPhi)<1.e-2 || CAMath::Abs(cP)<1.e-2 ) return 0;
573
15d2e9cf 574 //Float_t J[5][5] = { { j0, 0, 0, 0, 0 }, // Y
575 // { 0, 1, 0, 0, 0 }, // Z
576 // { 0, 0, j2, 0, 0 }, // SinPhi
577 // { 0, 0, 0, 1, 0 }, // DzDs
578 // { 0, 0, 0, 0, 1 } }; // Kappa
579
693d2443 580 Float_t j0 = cP/cosPhi;
581 Float_t j2 = cosPhi/cP;
582 Float_t d[2] = {Y() - y0, SinPhi() - sP};
583
15d2e9cf 584 SetX( x0*cA + y0*sA );
585 SetY(-x0*sA + y0*cA + j0*d[0] );
693d2443 586 t0.SetCosPhi( cosPhi );
587 t0.SetSinPhi( sinPhi );
588
693d2443 589 SetSinPhi( sinPhi + j2*d[1] );
590
591 fC[0]*= j0*j0;
592 fC[1]*= j0;
593 fC[3]*= j0;
594 fC[6]*= j0;
595 fC[10]*= j0;
596
597 fC[3]*= j2;
598 fC[4]*= j2;
599 fC[5]*= j2*j2;
600 fC[8]*= j2;
601 fC[12]*= j2;
602
603 return 1;
604}
605
606
15d2e9cf 607GPUd() Bool_t AliHLTTPCCATrackParam::Filter( Float_t y, Float_t z, Float_t err2Y, Float_t err2Z, Float_t maxSinPhi )
4687b8fc 608{
609 //* Add the y,z measurement with the Kalman filter
610
611 Float_t
612 c00 = fC[ 0],
613 c11 = fC[ 2],
614 c20 = fC[ 3],
615 c31 = fC[ 7],
616 c40 = fC[10];
617
618 err2Y+=c00;
619 err2Z+=c11;
620
621 Float_t
622 z0 = y-fP[0],
623 z1 = z-fP[1];
624
625 if( err2Y < 1.e-8 || err2Z<1.e-8 ) return 0;
626
627 Float_t mS0 = 1./err2Y;
628 Float_t mS2 = 1./err2Z;
629
630 // K = CHtS
631
632 Float_t k00, k11, k20, k31, k40;
633
634 k00 = c00*mS0;
635 k20 = c20*mS0;
636 k40 = c40*mS0;
637
638 k11 = c11*mS2;
639 k31 = c31*mS2;
640
641 Float_t sinPhi = fP[2] + k20*z0 ;
693d2443 642
15d2e9cf 643 if( maxSinPhi>0 && CAMath::Abs(sinPhi)>= maxSinPhi ) return 0;
693d2443 644
693d2443 645 fNDF += 2;
646 fChi2 += mS0*z0*z0 + mS2*z1*z1 ;
647
648 fP[ 0]+= k00*z0 ;
649 fP[ 1]+= k11*z1 ;
650 fP[ 2] = sinPhi ;
651 fP[ 3]+= k31*z1 ;
652 fP[ 4]+= k40*z0 ;
653
654 fC[ 0]-= k00*c00 ;
655 fC[ 3]-= k20*c00 ;
656 fC[ 5]-= k20*c20 ;
657 fC[10]-= k40*c00 ;
658 fC[12]-= k40*c20 ;
659 fC[14]-= k40*c40 ;
660
661 fC[ 2]-= k11*c11 ;
662 fC[ 7]-= k31*c11 ;
663 fC[ 9]-= k31*c31 ;
664
693d2443 665 return 1;
666}
667
668
eb30eb49 669
eb30eb49 670
00d07bcd 671#if !defined(HLTCA_GPUCODE)
4687b8fc 672#include <iostream>
00d07bcd 673#endif
674
675GPUd() void AliHLTTPCCATrackParam::Print() const
eb30eb49 676{
00d07bcd 677 //* print parameters
678
679#if !defined(HLTCA_GPUCODE)
15d2e9cf 680 std::cout<<"track: x="<<GetX()<<" c="<<GetSignCosPhi()<<", P= "<<GetY()<<" "<<GetZ()<<" "<<GetSinPhi()<<" "<<GetDzDs()<<" "<<GetQPt()<<std::endl;
681 std::cout<<"errs2: "<<GetErr2Y()<<" "<<GetErr2Z()<<" "<<GetErr2SinPhi()<<" "<<GetErr2DzDs()<<" "<<GetErr2QPt()<<std::endl;
00d07bcd 682#endif
d54804bf 683}
15d2e9cf 684