STEER/AliExternalTrackParam.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //                                                                           //
  20 // track parameters in "external" format                                     //
  21 //                                                                           //
  22 // The track parameters are:                                                //
  23 // - local y coordinate                                                      //
  24 // - local z coordinate                                                      //
  25 // - sin of azimutal angle                                                   //
  26 // - tan of dip angle                                                        //
  27 // - charge/pt                                                               //
  28 // The parametrisation is given at the local x coordinate fX and the         //
  29 // azimuthal angle fAlpha.                                                   //
  30 //                                                                           //
  31 // The external parametrisation can be used to exchange track parameters     //
  32 // between different detectors.                                              //
  33 //                                                                           //
  34 ///////////////////////////////////////////////////////////////////////////////
  35
  36 #include <TMath.h>
  37 #include <TVector3.h>
  38
  39 #include "AliExternalTrackParam.h"
  40 #include "AliKalmanTrack.h"
  41 #include "AliTrackReference.h"
  42 #include "AliLog.h"
  43
  44 ClassImp(AliExternalTrackParam)
  45
  46 const AliMagF *AliExternalTrackParam::fgkFieldMap=0;
  47 Double_t AliExternalTrackParam::fgConvConst=0.;
  48
  49
  50 //_____________________________________________________________________________
  51 AliExternalTrackParam::AliExternalTrackParam() :
  52   fMass(-1),
  53   fX(0),
  54   fAlpha(0),
  55   fLocalConvConst(0)
  56 {
  57   //
  58   // default constructor
  59   //
  60   for (Int_t i = 0; i < 5; i++) fParam[i] = 0;
  61   for (Int_t i = 0; i < 15; i++) fCovar[i] = 0;
  62 }
  63
  64 //_____________________________________________________________________________
  65 AliExternalTrackParam::AliExternalTrackParam(Double_t x, Double_t alpha,
  66                                              const Double_t param[5],
  67                                              const Double_t covar[15]) :
  68   fMass(-1),
  69   fX(x),
  70   fAlpha(alpha),
  71   fLocalConvConst(0)
  72 {
  73   //
  74   // create external track parameters from given arguments
  75   //
  76   for (Int_t i = 0; i < 5; i++) fParam[i] = param[i];
  77   for (Int_t i = 0; i < 15; i++) fCovar[i] = covar[i];
  78 }
  79
  80 //_____________________________________________________________________________
  81 AliExternalTrackParam::AliExternalTrackParam(const AliKalmanTrack& track) :
  82   fMass(track.GetMass()),
  83   fX(0),
  84   fAlpha(track.GetAlpha()),
  85   fLocalConvConst(0)
  86 {
  87   //
  88   //
  89   track.GetExternalParameters(fX,fParam);
  90   track.GetExternalCovariance(fCovar);
  91   SaveLocalConvConst();
  92 }
  93
  94
  95 AliExternalTrackParam * AliExternalTrackParam::MakeTrack(const AliTrackReference *ref, Double_t mass)
  96 {
  97   //
  98   // Make dummy track from the track reference
  99   // negative mass means opposite charge
 100   //
 101   Double_t xx[5];
 102   Double_t cc[15];
 103   for (Int_t i=0;i<15;i++) cc[i]=0;
 104   Double_t x = ref->X(), y = ref->Y(), z = ref->Z();
 105   Double_t alpha = TMath::ATan2(y,x);
 106   Double_t xr = TMath::Sqrt(x*x+y*y);
 107   xx[0] = 0;
 108   xx[1] = z;
 109   xx[3] = ref->Pz()/ref->Pt();
 110   Float_t b[3];
 111   Float_t xyz[3]={x,y,z};
 112   Float_t convConst = 0;
 113   (AliKalmanTrack::GetFieldMap())->Field(xyz,b);
 114   convConst=1000/0.299792458/(1e-13 - b[2]);
 115   xx[4] = 1./(convConst*ref->Pt()); // curvature rpresentation
 116   if (mass<0) xx[4]*=-1.;  // negative mass - negative direction
 117   Double_t alphap = TMath::ATan2(ref->Py(),ref->Px())-alpha;
 118   if (alphap> TMath::Pi()) alphap-=TMath::Pi();
 119   if (alphap<-TMath::Pi()) alphap+=TMath::Pi();
 120   xx[2] = TMath::Sin(alphap);
 121   xx[4]*=convConst;   // 1/pt representation
 122   //  AliExternalTrackParam * track = new  AliExternalTrackParam(xx,cc,xr,alpha);
 123   AliExternalTrackParam * track = new  AliExternalTrackParam(xr,alpha,xx,cc);
 124   track->fMass = TMath::Abs(mass);
 125   //track->StartTimeIntegral();
 126   track->SaveLocalConvConst();
 127   return track;
 128 }
 129
 130
 131
 132 //_____________________________________________________________________________
 133 const Double_t* AliExternalTrackParam::GetParameter() const
 134 {
 135 // get a pointer to the array of track parameters
 136
 137   return fParam;
 138 }
 139
 140 //_____________________________________________________________________________
 141 const Double_t* AliExternalTrackParam::GetCovariance() const
 142 {
 143 // get a pointer to the array of the track parameter covariance matrix
 144
 145   return fCovar;
 146 }
 147
 148 //_____________________________________________________________________________
 149 AliExternalTrackParam* AliExternalTrackParam::CreateExternalParam() const
 150 {
 151 // copy this instance
 152
 153   return new AliExternalTrackParam(fX, fAlpha, fParam, fCovar);
 154 }
 155
 156 //_____________________________________________________________________________
 157 void AliExternalTrackParam::ResetCovariance(Double_t factor,
 158                                             Bool_t clearOffDiagonal)
 159 {
 160 // reset the covariance matrix ("forget" track history)
 161
 162   Int_t k = 0;
 163   for (Int_t i = 0; i < 5; i++) {
 164     for (Int_t j = 0; j < i; j++) {  // off diagonal elements
 165       if (clearOffDiagonal) {
 166         fCovar[k++] = 0;
 167       } else {
 168         fCovar[k++] *= factor;
 169       }
 170     }
 171     fCovar[k++] *= factor;     // diagonal elements
 172   }
 173 }
 174
 175
 176 //_____________________________________________________________________________
 177 Bool_t AliExternalTrackParam::PropagateTo(Double_t xk, Double_t x0, Double_t rho)
 178 {
 179   //
 180   // Propagate the track parameters to the given x coordinate assuming vacuum.
 181   // If length is not NULL, the change of track length is added to it.
 182   //
 183
 184   Double_t lcc=GetLocalConvConst();
 185   Double_t cur = fParam[4]/lcc;
 186   Double_t x1=fX, x2=xk, dx=x2-x1;
 187   Double_t f1=fParam[2], f2=f1 + cur*dx;
 188   if (TMath::Abs(f2) >= 0.98) {
 189     // MI change  - don't propagate highly inclined tracks
 190     //              covariance matrix distorted
 191     return kFALSE;
 192   }
 193
 194   // old position [SR, GSI, 17.02.2003]
 195   Double_t oldX = fX, oldY = fParam[0], oldZ = fParam[1];
 196   Double_t r1=sqrt(1.- f1*f1), r2=sqrt(1.- f2*f2);
 197   fParam[0] += dx*(f1+f2)/(r1+r2);
 198   fParam[1] += dx*(f1+f2)/(f1*r2 + f2*r1)*fParam[3];
 199   fParam[2] += dx*cur;
 200   // transform error matrix to the curvature
 201   fCovar[10]/=lcc;
 202   fCovar[11]/=lcc;
 203   fCovar[12]/=lcc;
 204   fCovar[13]/=lcc;
 205   fCovar[14]/=lcc*lcc;
 206
 207   //f = F - 1
 208
 209   Double_t f02=    dx/(r1*r1*r1);
 210   Double_t f04=0.5*dx*dx/(r1*r1*r1);
 211   Double_t f12=    dx*fParam[3]*f1/(r1*r1*r1);
 212   Double_t f14=0.5*dx*dx*fParam[3]*f1/(r1*r1*r1);
 213   Double_t f13=    dx/r1;
 214   Double_t f24=    dx;
 215
 216   //b = C*ft
 217   Double_t b00=f02*fCovar[3] + f04*fCovar[10], b01=f12*fCovar[3] + f14*fCovar[10] + f13*fCovar[6];
 218   Double_t b02=f24*fCovar[10];
 219   Double_t b10=f02*fCovar[4] + f04*fCovar[11], b11=f12*fCovar[4] + f14*fCovar[11] + f13*fCovar[7];
 220   Double_t b12=f24*fCovar[11];
 221   Double_t b20=f02*fCovar[5] + f04*fCovar[12], b21=f12*fCovar[5] + f14*fCovar[12] + f13*fCovar[8];
 222   Double_t b22=f24*fCovar[12];
 223   Double_t b40=f02*fCovar[12] + f04*fCovar[14], b41=f12*fCovar[12] + f14*fCovar[14] + f13*fCovar[13];
 224   Double_t b42=f24*fCovar[14];
 225   Double_t b30=f02*fCovar[8] + f04*fCovar[13], b31=f12*fCovar[8] + f14*fCovar[13] + f13*fCovar[9];
 226   Double_t b32=f24*fCovar[13];
 227
 228   //a = f*b = f*C*ft
 229   Double_t a00=f02*b20+f04*b40,a01=f02*b21+f04*b41,a02=f02*b22+f04*b42;
 230   Double_t a11=f12*b21+f14*b41+f13*b31,a12=f12*b22+f14*b42+f13*b32;
 231   Double_t a22=f24*b42;
 232
 233   //F*C*Ft = C + (b + bt + a)
 234   fCovar[0] += b00 + b00 + a00;
 235   fCovar[1] += b10 + b01 + a01;
 236   fCovar[2] += b11 + b11 + a11;
 237   fCovar[3] += b20 + b02 + a02;
 238   fCovar[4] += b21 + b12 + a12;
 239   fCovar[5] += b22 + b22 + a22;
 240   fCovar[6] += b30;
 241   fCovar[7] += b31;
 242   fCovar[8] += b32;
 243   fCovar[10] += b40;
 244   fCovar[11] += b41;
 245   fCovar[12] += b42;
 246
 247   fX=x2;
 248
 249   //Change of the magnetic field *************
 250   SaveLocalConvConst();
 251   // transform back error matrix from curvature to the 1/pt
 252   fCovar[10]*=lcc;
 253   fCovar[11]*=lcc;
 254   fCovar[12]*=lcc;
 255   fCovar[13]*=lcc;
 256   fCovar[14]*=lcc*lcc;
 257
 258   Double_t dist = TMath::Sqrt((fX-oldX)*(fX-oldX)+(fParam[0]-oldY)*(fParam[0]-oldY)+
 259                               (fParam[1]-oldZ)*(fParam[1]-oldZ));
 260   if (!CorrectForMaterial(dist,x0,rho)) return 0;
 261
 262   // Integrated Time [SR, GSI, 17.02.2003]
 263  //  if (IsStartedTimeIntegral() && fX>oldX) {
 264 //     Double_t l2 = (fX-oldX)*(fX-oldX)+(fParam[0]-oldY)*(fParam[0]-oldY)+
 265 //                   (fParam[1]-oldZ)*(fParam[1]-oldZ);
 266 //     AddTimeStep(TMath::Sqrt(l2));
 267 //   }
 268   //
 269
 270   return kTRUE;
 271 }
 272
 273 Bool_t     AliExternalTrackParam::PropagateToDCA(Double_t xd, Double_t yd,  Double_t x0, Double_t rho){
 274   //
 275   // Propagate the track parameters to the nearest point of given xv yv coordinate
 276   //
 277   Double_t a=fAlpha;
 278   Double_t cs=TMath::Cos(a),sn=TMath::Sin(a);
 279
 280   Double_t xv= xd*cs + yd*sn;
 281   Double_t yv=-xd*sn + yd*cs;   // vertex position in local frame
 282   //
 283   Double_t c=fParam[4]/GetLocalConvConst(), snp=fParam[2];
 284   //
 285   Double_t x=fX, y=fParam[1];
 286   Double_t tgfv=-(c*(x-xv)-snp)/(c*(y-yv) + TMath::Sqrt(1.-snp*snp));
 287   Double_t fv=TMath::ATan(tgfv);
 288   cs=TMath::Cos(fv); sn=TMath::Sin(fv);
 289   x = xv*cs + yv*sn;
 290   yv=-xv*sn + yv*cs; xv=x;
 291   RotateTo(fv+a);
 292   return PropagateTo(xv,x0,rho);
 293 }
 294
 295
 296 //_____________________________________________________________________________
 297 Bool_t AliExternalTrackParam::RotateTo(Double_t alp)
 298 {
 299   // Rotate the reference axis for the parametrisation to the given angle.
 300   //
 301   Double_t  x=fX;
 302   Double_t p0=fParam[0];
 303   //
 304   if      (alp < -TMath::Pi()) alp += 2*TMath::Pi();
 305   else if (alp >= TMath::Pi()) alp -= 2*TMath::Pi();
 306   Double_t ca=TMath::Cos(alp-fAlpha), sa=TMath::Sin(alp-fAlpha);
 307   Double_t sf=fParam[2], cf=TMath::Sqrt(1.- fParam[2]*fParam[2]);
 308   // **** rotation **********************
 309
 310   fAlpha = alp;
 311   fX =  x*ca + p0*sa;
 312   fParam[0]= -x*sa + p0*ca;
 313   fParam[2]=  sf*ca - cf*sa;
 314   Double_t rr=(ca+sf/cf*sa);
 315   //
 316   fCovar[0] *= (ca*ca);
 317   fCovar[1] *= ca;
 318   fCovar[3] *= ca*rr;
 319   fCovar[6] *= ca;
 320   fCovar[10] *= ca;
 321   fCovar[4] *= rr;
 322   fCovar[5] *= rr*rr;
 323   fCovar[7] *= rr;
 324   fCovar[11] *= rr;
 325   return kTRUE;
 326 }
 327
 328 //_____________________________________________________________________________
 329 Bool_t AliExternalTrackParam::CorrectForMaterial(Double_t d, Double_t x0, Double_t rho)
 330 {
 331   //
 332   // Take into account material effects assuming:
 333   // x0  - mean rad length
 334   // rho - mean density
 335
 336   //
 337   // multiple scattering
 338   //
 339   if (fMass<=0) {
 340     AliError("Non-positive mass");
 341     return kFALSE;
 342   }
 343   Double_t p2=(1.+ fParam[3]*fParam[3])/(fParam[4]*fParam[4]);
 344   Double_t beta2=p2/(p2 + fMass*fMass);
 345   Double_t theta2=14.1*14.1/(beta2*p2*1e6)*d/x0*rho;
 346   //
 347   fCovar[5] += theta2*(1.- fParam[2]*fParam[2])*(1. + fParam[3]*fParam[3]);
 348   fCovar[9] += theta2*(1. + fParam[3]*fParam[3])*(1. + fParam[3]*fParam[3]);
 349   fCovar[13] += theta2*fParam[3]*fParam[4]*(1. + fParam[3]*fParam[3]);
 350   fCovar[14] += theta2*fParam[3]*fParam[4]*fParam[3]*fParam[4];
 351
 352   Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2+1e-10)) - beta2)*d*rho;
 353   fParam[4] *=(1.- TMath::Sqrt(p2+fMass*fMass)/p2*dE);
 354   //
 355   Double_t sigmade = 0.02*TMath::Sqrt(TMath::Abs(dE));   // energy loss fluctuation
 356   Double_t sigmac2 = sigmade*sigmade*fParam[4]*fParam[4]*(p2+fMass*fMass)/(p2*p2);
 357   fCovar[14] += sigmac2;
 358   //
 359   //
 360
 361
 362
 363   return kTRUE;
 364 }
 365
 366 //_____________________________________________________________________________
 367 Bool_t AliExternalTrackParam::GetProlongationAt(Double_t xk,
 368                                                 Double_t& y,
 369                                                 Double_t& z) const
 370 {
 371   //
 372   // Get the local y and z coordinates at the given x value
 373   //
 374   Double_t lcc=GetLocalConvConst();
 375   Double_t cur = fParam[4]/lcc;
 376   Double_t x1=fX, x2=xk, dx=x2-x1;
 377   Double_t f1=fParam[2], f2=f1 + cur*dx;
 378   //
 379   if (TMath::Abs(f2) >= 0.98) {
 380     // MI change  - don't propagate highly inclined tracks
 381     //              covariance matrix distorted
 382     return kFALSE;
 383   }
 384   Double_t r1=sqrt(1.- f1*f1), r2=sqrt(1.- f2*f2);
 385   y = fParam[0] + dx*(f1+f2)/(r1+r2);
 386   z = fParam[1] + dx*(f1+f2)/(f1*r2 + f2*r1)*fParam[3];
 387   return kTRUE;
 388 }
 389
 390 //_____________________________________________________________________________
 391 Double_t AliExternalTrackParam::GetXAtVertex(Double_t /*x*/,
 392                                              Double_t /*y*/) const
 393 {
 394 // Get the x coordinate at the given vertex (x,y)
 395 //
 396 // NOT IMPLEMENTED for this class
 397
 398   return 0;
 399 }
 400
 401
 402 // //_____________________________________________________________________________
 403 // Double_t AliExternalTrackParam::GetPredictedChi2(const AliCluster* /*cluster*/)
 404 // {
 405 // // calculate the chi2 contribution of the given cluster
 406 // //
 407 // // NOT IMPLEMENTED for this class
 408
 409 //   return -1;
 410 // }
 411
 412 // //_____________________________________________________________________________
 413 // Bool_t AliExternalTrackParam::Update(const AliCluster* /*cluster*/)
 414 // {
 415 // // update the track parameters using the position and error
 416 // // of the given cluster
 417 // //
 418 // // NOT IMPLEMENTED for this class
 419
 420 //   return kFALSE;
 421 // }
 422
 423
 424 //_____________________________________________________________________________
 425 Double_t AliExternalTrackParam::SigmaPhi() const
 426 {
 427 // get the error of the azimuthal angle
 428
 429   return TMath::Sqrt(TMath::Abs(fCovar[5] / (1. - fParam[2]*fParam[2])));
 430 }
 431
 432 //_____________________________________________________________________________
 433 Double_t AliExternalTrackParam::SigmaTheta() const
 434 {
 435 // get the error of the polar angle
 436
 437   return TMath::Sqrt(TMath::Abs(fCovar[9])) / (1. + fParam[3]*fParam[3]);
 438 }
 439
 440 //_____________________________________________________________________________
 441 Double_t AliExternalTrackParam::SigmaPt() const
 442 {
 443 // get the error of the transversal component of the momentum
 444
 445   return TMath::Sqrt(fCovar[14]) / TMath::Abs(fParam[4]);
 446 }
 447
 448 //_____________________________________________________________________________
 449 TVector3 AliExternalTrackParam::Momentum() const
 450 {
 451 // get the momentum vector
 452
 453   Double_t phi = TMath::ASin(fParam[2]) + fAlpha;
 454   Double_t pt = 1. / TMath::Abs(fParam[4]);
 455   return TVector3(pt * TMath::Cos(phi),
 456                   pt * TMath::Sin(phi),
 457                   pt * fParam[3]);
 458 }
 459
 460 //_____________________________________________________________________________
 461 TVector3 AliExternalTrackParam::Position() const
 462 {
 463 // get the current spatial position in global coordinates
 464
 465   return TVector3(fX * TMath::Cos(fAlpha) - fParam[0] * TMath::Sin(fAlpha),
 466                   fX * TMath::Sin(fAlpha) + fParam[0] * TMath::Cos(fAlpha),
 467                   fParam[1]);
 468 }
 469
 470
 471 //_____________________________________________________________________________
 472 void AliExternalTrackParam::Print(Option_t* /*option*/) const
 473 {
 474 // print the parameters and the covariance matrix
 475
 476   printf("AliExternalTrackParam: x = %-12g  alpha = %-12g\n", fX, fAlpha);
 477   printf("  parameters: %12g %12g %12g %12g %12g\n",
 478          fParam[0], fParam[1], fParam[2], fParam[3], fParam[4]);
 479   printf("  covariance: %12g\n", fCovar[0]);
 480   printf("              %12g %12g\n", fCovar[1], fCovar[2]);
 481   printf("              %12g %12g %12g\n", fCovar[3], fCovar[4], fCovar[5]);
 482   printf("              %12g %12g %12g %12g\n",
 483          fCovar[6], fCovar[7], fCovar[8], fCovar[9]);
 484   printf("              %12g %12g %12g %12g %12g\n",
 485          fCovar[10], fCovar[11], fCovar[12], fCovar[13], fCovar[14]);
 486 }