HLT/TPCLib/merger-ca/AliHLTTPCGMSliceTrack.cxx

   1 // $Id: AliHLTTPCGMSliceTrack.cxx 41769 2010-06-16 13:58:00Z sgorbuno $
   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 //                  for The ALICE HLT Project.                              *
   8 //                                                                          *
   9 // Permission to use, copy, modify and distribute this software and its     *
  10 // documentation strictly for non-commercial purposes is hereby granted     *
  11 // without fee, provided that the above copyright notice appears in all     *
  12 // copies and that both the copyright notice and this permission notice     *
  13 // appear in the supporting documentation. The authors make no claims       *
  14 // about the suitability of this software for any purpose. It is            *
  15 // provided "as is" without express or implied warranty.                    *
  16 //                                                                          *
  17 //***************************************************************************
  18
  19 #include "AliHLTTPCGMSliceTrack.h"
  20 #include "AliHLTTPCCAMath.h"
  21 #include "AliHLTTPCGMBorderTrack.h"
  22 #include "AliHLTTPCGMTrackLinearisation.h"
  23 #include "Riostream.h"
  24 #include "AliHLTTPCCAParam.h"
  25 #include <cmath>
  26
  27
  28
  29 bool AliHLTTPCGMSliceTrack::FilterErrors( AliHLTTPCCAParam &param, float maxSinPhi )
  30 {
  31   float lastX = fOrigTrack->Cluster(fOrigTrack->NClusters()-1 ).GetX();
  32
  33   const int N = 3;
  34   const float *cy = param.GetParamS0Par(0,0);
  35   const float *cz = param.GetParamS0Par(1,0);
  36
  37   float bz = -param.ConstBz();
  38   const float kZLength = 250.f - 0.275f;
  39
  40   float trDzDs2 = fDzDs*fDzDs;
  41   float k  = fQPt*bz;
  42   float dx = .33*(lastX - fX);
  43   float kdx = k*dx;
  44   float dxBz = dx * bz;
  45   float kdx205 = 2.f+kdx*kdx*0.5f;
  46
  47   {
  48     float secPhi2 = fSecPhi*fSecPhi;
  49     float zz = fabs( kZLength - fabs(fZ) );
  50     float zz2 = zz*zz;
  51     float angleY2 = secPhi2 - 1.f;
  52     float angleZ2 = trDzDs2 * secPhi2 ;
  53
  54     float cy0 = cy[0] + cy[1]*zz + cy[3]*zz2;
  55     float cy1 = cy[2] + cy[5]*zz;
  56     float cy2 = cy[4];
  57     float cz0 = cz[0] + cz[1]*zz + cz[3]*zz2;
  58     float cz1 = cz[2] + cz[5]*zz;
  59     float cz2 = cz[4];
  60
  61     fC0 = fabs( cy0 + angleY2 * ( cy1 + angleY2*cy2 ) );
  62     fC2 = fabs( cz0 + angleZ2 * ( cz1 + angleZ2*cz2 ) );
  63
  64     fC3 = 0;
  65     fC5 = 1;
  66     fC7 = 0;
  67     fC9 = 1;
  68     fC10 = 0;
  69     fC12 = 0;
  70     fC14 = 10;
  71   }
  72
  73   for( int iStep=0; iStep<N; iStep++ ){
  74
  75     float err2Y, err2Z;
  76
  77     { // transport block
  78
  79       float ex = fCosPhi;
  80       float ey = fSinPhi;
  81       float ey1 = kdx + ey;
  82       if( fabs( ey1 ) > maxSinPhi ) return 0;
  83
  84       float ss = ey + ey1;
  85       float ex1 = sqrt(1.f - ey1*ey1);
  86
  87       float cc = ex + ex1;
  88       float dxcci = dx / cc;
  89
  90       float dy = dxcci * ss;
  91       float norm2 = 1.f + ey*ey1 + ex*ex1;
  92       float dl = dxcci * sqrt( norm2 + norm2 );
  93
  94       float dS;
  95       {
  96         float dSin = 0.5f*k*dl;
  97         float a = dSin*dSin;
  98         const float k2 = 1.f/6.f;
  99         const float k4 = 3.f/40.f;
 100         dS = dl + dl*a*(k2 + a*(k4 ));//+ k6*a) );
 101       }
 102
 103       float dz = dS * fDzDs;
 104       float ex1i =1.f/ex1;
 105       float z = fZ+ dz;
 106       {
 107         float secPhi2 = ex1i*ex1i;
 108         float zz = fabs( kZLength - fabs(z) );
 109         float zz2 = zz*zz;
 110         float angleY2 = secPhi2 - 1.f;
 111         float angleZ2 = trDzDs2 * secPhi2 ;
 112
 113         float cy0 = cy[0] + cy[1]*zz + cy[3]*zz2;
 114         float cy1 = cy[2] + cy[5]*zz;
 115         float cy2 = cy[4];
 116         float cz0 = cz[0] + cz[1]*zz + cz[3]*zz2;
 117         float cz1 = cz[2] + cz[5]*zz;
 118         float cz2 = cz[4];
 119
 120         err2Y = fabs( cy0 + angleY2 * ( cy1 + angleY2*cy2 ) );
 121         err2Z = fabs( cz0 + angleZ2 * ( cz1 + angleZ2*cz2 ) );
 122       }
 123
 124       float hh = kdx205 * dxcci*ex1i;
 125       float h2 = hh * fSecPhi;
 126
 127       fX+=dx;
 128       fY+= dy;
 129       fZ+= dz;
 130       fSinPhi = ey1;
 131       fCosPhi = ex1;
 132       fSecPhi = ex1i;
 133
 134       float h4 = bz*dxcci*hh;
 135
 136       float c20 = fC3;
 137       float c22 = fC5;
 138       float c31 = fC7;
 139       float c33 = fC9;
 140       float c40 = fC10;
 141       float c42 = fC12;
 142       float c44 = fC14;
 143
 144       float c20ph4c42 =  c20 + h4*c42;
 145       float h2c22 = h2*c22;
 146       float h4c44 = h4*c44;
 147       float n7 = c31 + dS*c33;
 148       float n10 = c40 + h2*c42 + h4c44;
 149       float n12 = c42 + dxBz*c44;
 150
 151
 152       fC0+= h2*h2c22 + h4*h4c44 + 2.f*( h2*c20ph4c42  + h4*c40 );
 153
 154       fC3 = c20ph4c42 + h2c22  + dxBz*n10;
 155       fC10 = n10;
 156
 157       fC5 = c22 + dxBz*( c42 + n12 );
 158       fC12 = n12;
 159
 160       fC2+= dS*(c31 + n7);
 161       fC7 = n7;
 162
 163    } // end transport block
 164
 165
 166     // Filter block
 167
 168     float
 169       c00 = fC0,
 170       c11 = fC2,
 171       c20 = fC3,
 172       c31 = fC7,
 173       c40 = fC10;
 174
 175     float mS0 = 1.f/(err2Y + c00);
 176     float mS2 = 1.f/(err2Z + c11);
 177
 178     // K = CHtS
 179
 180     float k00, k11, k20, k31, k40;
 181
 182     k00 = c00 * mS0;
 183     k20 = c20 * mS0;
 184     k40 = c40 * mS0;
 185
 186     fC0 -= k00 * c00 ;
 187     fC5 -= k20 * c20 ;
 188     fC10 -= k00 * c40 ;
 189     fC12 -= k40 * c20 ;
 190     fC3 -= k20 * c00 ;
 191     fC14 -= k40 * c40 ;
 192
 193     k11 = c11 * mS2;
 194     k31 = c31 * mS2;
 195
 196     fC7 -= k31 * c11;
 197     fC2 -= k11 * c11;
 198     fC9 -= k31 * c31;
 199   }
 200
 201   //* Check that the track parameters and covariance matrix are reasonable
 202
 203   bool ok = 1;
 204
 205   const float *c = &fX;
 206   for ( int i = 0; i < 17; i++ ) ok = ok && finite( c[i] );
 207
 208   if ( fC0 <= 0.f || fC2 <= 0.f || fC5 <= 0.f || fC9 <= 0.f || fC14 <= 0.f
 209        || fC0 > 5.f || fC2 > 5.f || fC5 > 2.f || fC9 > 2.f             ) ok = 0;
 210
 211   if( ok ){
 212     ok = ok
 213       && ( fC3*fC3<=fC5*fC0 )
 214       && ( fC7*fC7<=fC9*fC2 )
 215       && ( fC10*fC10<=fC14*fC0 )
 216       && ( fC12*fC12<=fC14*fC5 );
 217   }
 218
 219   return ok;
 220 }
 221
 222
 223
 224 bool AliHLTTPCGMSliceTrack::TransportToX( float x, float Bz, AliHLTTPCGMBorderTrack &b, float maxSinPhi ) const
 225 {
 226   Bz = -Bz;
 227   float ex = fCosPhi;
 228   float ey = fSinPhi;
 229   float k  = fQPt*Bz;
 230   float dx = x - fX;
 231   float ey1 = k*dx + ey;
 232
 233   if( fabs( ey1 ) > maxSinPhi ) return 0;
 234
 235   float ex1 = sqrt( 1.f - ey1 * ey1 );
 236   float dxBz = dx * Bz;
 237
 238   float ss = ey + ey1;
 239   float cc = ex + ex1;
 240   float dxcci = dx / cc;
 241   float norm2 = 1.f + ey*ey1 + ex*ex1;
 242
 243   float dy = dxcci * ss;
 244
 245   float dS;
 246   {
 247     float dl = dxcci * sqrt( norm2 + norm2 );
 248     float dSin = 0.5f*k*dl;
 249     float a = dSin*dSin;
 250     const float k2 = 1.f/6.f;
 251     const float k4 = 3.f/40.f;
 252     //const float k6 = 5.f/112.f;
 253     dS = dl + dl*a*(k2 + a*(k4 ));//+ k6*a) );
 254   }
 255
 256   float ex1i = 1.f/ex1;
 257   float dz = dS * fDzDs;
 258
 259   float hh = dxcci*ex1i*norm2;
 260   float h2 = hh *fSecPhi;
 261   float h4 = Bz*dxcci*hh;
 262
 263   float c20 = fC3;
 264   float c22 = fC5;
 265   float c31 = fC7;
 266   float c33 = fC9;
 267   float c40 = fC10;
 268   float c42 = fC12;
 269   float c44 = fC14;
 270
 271   float c20ph4c42 =  c20 + h4*c42;
 272   float h2c22 = h2*c22;
 273   float h4c44 = h4*c44;
 274   float n7 = c31 + dS*c33;
 275
 276   b.SetPar(0, fY + dy );
 277   b.SetPar(1, fZ + dz );
 278   b.SetPar(2, ey1 );
 279   b.SetPar(3, fDzDs);
 280   b.SetPar(4, fQPt);
 281
 282   b.SetCov(0, fC0 + h2*h2c22 + h4*h4c44 + 2.f*( h2*c20ph4c42  + h4*c40 ));
 283   b.SetCov(1, fC2+ dS*(c31 + n7) );
 284   b.SetCov(2, c22 + dxBz*( c42 + c42 + dxBz*c44 ));
 285   b.SetCov(3, c33);
 286   b.SetCov(4, c44);
 287   b.SetCovD(0, c20ph4c42 + h2c22  + dxBz*(c40 + h2*c42 + h4c44) );
 288   b.SetCovD(1, n7 );
 289   return 1;
 290 }
 291
 292
 293
 294 bool AliHLTTPCGMSliceTrack::TransportToXAlpha( float newX, float sinAlpha, float cosAlpha, float Bz, AliHLTTPCGMBorderTrack &b, float maxSinPhi ) const
 295 {
 296   //*
 297
 298   float c00 = fC0;
 299   float c11 = fC2;
 300   float c20 = fC3;
 301   float c22 = fC5;
 302   float c31 = fC7;
 303   float c33 = fC9;
 304   float c40 = fC10;
 305   float c42 = fC12;
 306   float c44 = fC14;
 307
 308   float x,y;
 309   float z = fZ;
 310   float sinPhi = fSinPhi;
 311   float cosPhi = fCosPhi;
 312   float secPhi = fSecPhi;
 313   float dzds = fDzDs;
 314   float qpt = fQPt;
 315
 316   // Rotate the coordinate system in XY on the angle alpha
 317   {
 318     float sP = sinPhi, cP = cosPhi;
 319     cosPhi =  cP * cosAlpha + sP * sinAlpha;
 320     sinPhi = -cP * sinAlpha + sP * cosAlpha;
 321
 322     if ( CAMath::Abs( sinPhi ) > .999 || CAMath::Abs( cP ) < 1.e-2  ) return 0;
 323
 324     secPhi = 1./cosPhi;
 325     float j0 = cP *secPhi;
 326     float j2 = cosPhi / cP;
 327     x =   fX*cosAlpha +  fY*sinAlpha ;
 328     y =  -fX*sinAlpha +  fY*cosAlpha ;
 329
 330     c00 *= j0 * j0;
 331     c40 *= j0;
 332
 333     c22 *= j2 * j2;
 334     c42 *= j2;
 335     if( cosPhi < 0.f ){ // rotate to 180'
 336       cosPhi = -cosPhi;
 337       secPhi = -secPhi;
 338       sinPhi = -sinPhi;
 339       dzds = -dzds;
 340       qpt = -qpt;
 341       c20 = -c20;
 342       c31 = -c31;
 343       c40 = -c40;
 344    }
 345   }
 346
 347   Bz = -Bz;
 348   float ex = cosPhi;
 349   float ey = sinPhi;
 350   float k  = qpt*Bz;
 351   float dx = newX - x;
 352   float ey1 = k*dx + ey;
 353
 354   if( fabs( ey1 ) > maxSinPhi ) return 0;
 355
 356   float ex1 = sqrt( 1.f - ey1 * ey1 );
 357
 358   float dxBz = dx * Bz;
 359
 360   float ss = ey + ey1;
 361   float cc = ex + ex1;
 362   float dxcci = dx / cc;
 363   float norm2 = 1.f + ey*ey1 + ex*ex1;
 364
 365   float dy = dxcci * ss;
 366
 367   float dS;
 368   {
 369     float dl = dxcci * sqrt( norm2 + norm2 );
 370     float dSin = 0.5f*k*dl;
 371     float a = dSin*dSin;
 372     const float k2 = 1.f/6.f;
 373     const float k4 = 3.f/40.f;
 374     //const float k6 = 5.f/112.f;
 375     dS = dl + dl*a*(k2 + a*(k4 ));//+ k6*a) );
 376   }
 377
 378   float ex1i = 1.f/ex1;
 379   float dz = dS * dzds;
 380
 381   float hh = dxcci*ex1i*norm2;
 382   float h2 = hh * secPhi;
 383   float h4 = Bz*dxcci*hh;
 384
 385   float c20ph4c42 =  c20 + h4*c42;
 386   float h2c22 = h2*c22;
 387   float h4c44 = h4*c44;
 388   float n7 = c31 + dS*c33;
 389
 390   b.SetPar(0, y + dy );
 391   b.SetPar(1, z + dz );
 392   b.SetPar(2, ey1 );
 393   b.SetPar(3, dzds);
 394   b.SetPar(4, qpt);
 395
 396   b.SetCov(0, c00 + h2*h2c22 + h4*h4c44 + 2.f*( h2*c20ph4c42  + h4*c40 ));
 397   b.SetCov(1, c11 + dS*(c31 + n7) );
 398   b.SetCov(2, c22 + dxBz*( c42 + c42 + dxBz*c44 ));
 399   b.SetCov(3, c33);
 400   b.SetCov(4, c44);
 401   b.SetCovD(0, c20ph4c42 + h2c22  + dxBz*(c40 + h2*c42 + h4c44) );
 402   b.SetCovD(1, n7 );
 403
 404   return 1;
 405 }