STEER/ESD/AliTrackerBase.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: AliTrackerBase.cxx 38069 2009-12-24 16:56:18Z belikov $ */
  17
  18 //-------------------------------------------------------------------------
  19 //               Implementation of the AliTrackerBase class
  20 //                that is the base for the AliTracker class
  21 //                     Origin: Marian.Ivanov@cern.ch
  22 //-------------------------------------------------------------------------
  23 #include <TClass.h>
  24 #include <TMath.h>
  25 #include <TGeoManager.h>
  26
  27 #include "AliLog.h"
  28 #include "AliTrackerBase.h"
  29 #include "AliExternalTrackParam.h"
  30 #include "AliTrackPointArray.h"
  31 #include "TVectorD.h"
  32
  33 extern TGeoManager *gGeoManager;
  34
  35 ClassImp(AliTrackerBase)
  36
  37 AliTrackerBase::AliTrackerBase():
  38   TObject(),
  39   fX(0),
  40   fY(0),
  41   fZ(0),
  42   fSigmaX(0.005),
  43   fSigmaY(0.005),
  44   fSigmaZ(0.010)
  45 {
  46   //--------------------------------------------------------------------
  47   // The default constructor.
  48   //--------------------------------------------------------------------
  49   if (!TGeoGlobalMagField::Instance()->GetField())
  50     AliWarning("Field map is not set.");
  51 }
  52
  53 //__________________________________________________________________________
  54 AliTrackerBase::AliTrackerBase(const AliTrackerBase &atr):
  55   TObject(atr),
  56   fX(atr.fX),
  57   fY(atr.fY),
  58   fZ(atr.fZ),
  59   fSigmaX(atr.fSigmaX),
  60   fSigmaY(atr.fSigmaY),
  61   fSigmaZ(atr.fSigmaZ)
  62 {
  63   //--------------------------------------------------------------------
  64   // The default constructor.
  65   //--------------------------------------------------------------------
  66   if (!TGeoGlobalMagField::Instance()->GetField())
  67     AliWarning("Field map is not set.");
  68 }
  69
  70 //__________________________________________________________________________
  71 Double_t AliTrackerBase::GetBz()
  72 {
  73   AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
  74   if (!fld) return 0.5*kAlmost0Field;
  75   Double_t bz = fld->SolenoidField();
  76   return TMath::Sign(0.5*kAlmost0Field,bz) + bz;
  77 }
  78
  79 //__________________________________________________________________________
  80 Double_t AliTrackerBase::GetBz(const Double_t *r) {
  81   //------------------------------------------------------------------
  82   // Returns Bz (kG) at the point "r" .
  83   //------------------------------------------------------------------
  84   AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
  85   if (!fld) return  0.5*kAlmost0Field;
  86   Double_t bz = fld->GetBz(r);
  87   return  TMath::Sign(0.5*kAlmost0Field,bz) + bz;
  88 }
  89
  90 //__________________________________________________________________________
  91 void AliTrackerBase::GetBxByBz(const Double_t r[3], Double_t b[3]) {
  92   //------------------------------------------------------------------
  93   // Returns Bx, By and Bz (kG) at the point "r" .
  94   //------------------------------------------------------------------
  95   AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
  96   if (!fld) {
  97      b[0] = b[1] = 0.;
  98      b[2] = 0.5*kAlmost0Field;
  99      return;
 100   }
 101
 102   if (fld->IsUniform()) {
 103      b[0] = b[1] = 0.;
 104      b[2] = fld->SolenoidField();
 105   }  else {
 106      fld->Field(r,b);
 107   }
 108   b[2] = (TMath::Sign(0.5*kAlmost0Field,b[2]) + b[2]);
 109   return;
 110 }
 111
 112 Double_t AliTrackerBase::MeanMaterialBudget(const Double_t *start, const Double_t *end, Double_t *mparam)
 113 {
 114   //
 115   // Calculate mean material budget and material properties between
 116   //    the points "start" and "end".
 117   //
 118   // "mparam" - parameters used for the energy and multiple scattering
 119   //  corrections:
 120   //
 121   // mparam[0] - mean density: sum(x_i*rho_i)/sum(x_i) [g/cm3]
 122   // mparam[1] - equivalent rad length fraction: sum(x_i/X0_i) [adimensional]
 123   // mparam[2] - mean A: sum(x_i*A_i)/sum(x_i) [adimensional]
 124   // mparam[3] - mean Z: sum(x_i*Z_i)/sum(x_i) [adimensional]
 125   // mparam[4] - length: sum(x_i) [cm]
 126   // mparam[5] - Z/A mean: sum(x_i*Z_i/A_i)/sum(x_i) [adimensional]
 127   // mparam[6] - number of boundary crosses
 128   //
 129   //  Origin:  Marian Ivanov, Marian.Ivanov@cern.ch
 130   //
 131   //  Corrections and improvements by
 132   //        Andrea Dainese, Andrea.Dainese@lnl.infn.it,
 133   //        Andrei Gheata,  Andrei.Gheata@cern.ch
 134   //
 135
 136   mparam[0]=0; mparam[1]=1; mparam[2] =0; mparam[3] =0;
 137   mparam[4]=0; mparam[5]=0; mparam[6]=0;
 138   //
 139   Double_t bparam[6]; // total parameters
 140   Double_t lparam[6]; // local parameters
 141
 142   for (Int_t i=0;i<6;i++) bparam[i]=0;
 143
 144   if (!gGeoManager) {
 145     AliErrorClass("No TGeo\n");
 146     return 0.;
 147   }
 148   //
 149   Double_t length;
 150   Double_t dir[3];
 151   length = TMath::Sqrt((end[0]-start[0])*(end[0]-start[0])+
 152                        (end[1]-start[1])*(end[1]-start[1])+
 153                        (end[2]-start[2])*(end[2]-start[2]));
 154   mparam[4]=length;
 155   if (length<TGeoShape::Tolerance()) return 0.0;
 156   Double_t invlen = 1./length;
 157   dir[0] = (end[0]-start[0])*invlen;
 158   dir[1] = (end[1]-start[1])*invlen;
 159   dir[2] = (end[2]-start[2])*invlen;
 160
 161   // Initialize start point and direction
 162   TGeoNode *currentnode = 0;
 163   TGeoNode *startnode = gGeoManager->InitTrack(start, dir);
 164   if (!startnode) {
 165     AliErrorClass(Form("start point out of geometry: x %f, y %f, z %f",
 166                   start[0],start[1],start[2]));
 167     return 0.0;
 168   }
 169   TGeoMaterial *material = startnode->GetVolume()->GetMedium()->GetMaterial();
 170   lparam[0]   = material->GetDensity();
 171   lparam[1]   = material->GetRadLen();
 172   lparam[2]   = material->GetA();
 173   lparam[3]   = material->GetZ();
 174   lparam[4]   = length;
 175   lparam[5]   = lparam[3]/lparam[2];
 176   if (material->IsMixture()) {
 177     TGeoMixture * mixture = (TGeoMixture*)material;
 178     lparam[5] =0;
 179     Double_t sum =0;
 180     for (Int_t iel=0;iel<mixture->GetNelements();iel++){
 181       sum  += mixture->GetWmixt()[iel];
 182       lparam[5]+= mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel];
 183     }
 184     lparam[5]/=sum;
 185   }
 186
 187   // Locate next boundary within length without computing safety.
 188   // Propagate either with length (if no boundary found) or just cross boundary
 189   gGeoManager->FindNextBoundaryAndStep(length, kFALSE);
 190   Double_t step = 0.0; // Step made
 191   Double_t snext = gGeoManager->GetStep();
 192   // If no boundary within proposed length, return current density
 193   if (!gGeoManager->IsOnBoundary()) {
 194     mparam[0] = lparam[0];
 195     mparam[1] = lparam[4]/lparam[1];
 196     mparam[2] = lparam[2];
 197     mparam[3] = lparam[3];
 198     mparam[4] = lparam[4];
 199     return lparam[0];
 200   }
 201   // Try to cross the boundary and see what is next
 202   Int_t nzero = 0;
 203   while (length>TGeoShape::Tolerance()) {
 204     currentnode = gGeoManager->GetCurrentNode();
 205     if (snext<2.*TGeoShape::Tolerance()) nzero++;
 206     else nzero = 0;
 207     if (nzero>3) {
 208       // This means navigation has problems on one boundary
 209       // Try to cross by making a small step
 210       AliErrorClass("Cannot cross boundary\n");
 211       mparam[0] = bparam[0]/step;
 212       mparam[1] = bparam[1];
 213       mparam[2] = bparam[2]/step;
 214       mparam[3] = bparam[3]/step;
 215       mparam[5] = bparam[5]/step;
 216       mparam[4] = step;
 217       mparam[0] = 0.;             // if crash of navigation take mean density 0
 218       mparam[1] = 1000000;        // and infinite rad length
 219       return bparam[0]/step;
 220     }
 221     mparam[6]+=1.;
 222     step += snext;
 223     bparam[1]    += snext/lparam[1];
 224     bparam[2]    += snext*lparam[2];
 225     bparam[3]    += snext*lparam[3];
 226     bparam[5]    += snext*lparam[5];
 227     bparam[0]    += snext*lparam[0];
 228
 229     if (snext>=length) break;
 230     if (!currentnode) break;
 231     length -= snext;
 232     material = currentnode->GetVolume()->GetMedium()->GetMaterial();
 233     lparam[0] = material->GetDensity();
 234     lparam[1]  = material->GetRadLen();
 235     lparam[2]  = material->GetA();
 236     lparam[3]  = material->GetZ();
 237     lparam[5]   = lparam[3]/lparam[2];
 238     if (material->IsMixture()) {
 239       TGeoMixture * mixture = (TGeoMixture*)material;
 240       lparam[5]=0;
 241       Double_t sum =0;
 242       for (Int_t iel=0;iel<mixture->GetNelements();iel++){
 243         sum+= mixture->GetWmixt()[iel];
 244         lparam[5]+= mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel];
 245       }
 246       lparam[5]/=sum;
 247     }
 248     gGeoManager->FindNextBoundaryAndStep(length, kFALSE);
 249     snext = gGeoManager->GetStep();
 250   }
 251   mparam[0] = bparam[0]/step;
 252   mparam[1] = bparam[1];
 253   mparam[2] = bparam[2]/step;
 254   mparam[3] = bparam[3]/step;
 255   mparam[5] = bparam[5]/step;
 256   return bparam[0]/step;
 257 }
 258
 259
 260 Bool_t
 261 AliTrackerBase::PropagateTrackTo(AliExternalTrackParam *track, Double_t xToGo,
 262                                  Double_t mass, Double_t maxStep, Bool_t rotateTo, Double_t maxSnp, Int_t sign, Bool_t addTimeStep, Bool_t correctMaterialBudget){
 263   //----------------------------------------------------------------
 264   //
 265   // Propagates the track to the plane X=xk (cm) using the magnetic field map
 266   // and correcting for the crossed material.
 267   //
 268   // mass     - mass used in propagation - used for energy loss correction
 269   // maxStep  - maximal step for propagation
 270   //
 271   //  Origin: Marian Ivanov,  Marian.Ivanov@cern.ch
 272   //
 273   //----------------------------------------------------------------
 274   const Double_t kEpsilon = 0.00001;
 275   Double_t xpos     = track->GetX();
 276   Int_t dir         = (xpos<xToGo) ? 1:-1;
 277   //
 278   while ( (xToGo-xpos)*dir > kEpsilon){
 279     Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
 280     Double_t x    = xpos+step;
 281     Double_t xyz0[3],xyz1[3],param[7];
 282     track->GetXYZ(xyz0);   //starting global position
 283
 284     Double_t bz=GetBz(xyz0); // getting the local Bz
 285     if (!track->GetXYZAt(x,bz,xyz1)) return kFALSE;   // no prolongation
 286     xyz1[2]+=kEpsilon; // waiting for bug correction in geo
 287
 288     if (maxSnp>0 && TMath::Abs(track->GetSnpAt(x,bz)) >= maxSnp) return kFALSE;
 289     if (!track->PropagateTo(x,bz))  return kFALSE;
 290
 291     if (correctMaterialBudget){
 292       MeanMaterialBudget(xyz0,xyz1,param);
 293       Double_t xrho=param[0]*param[4], xx0=param[1];
 294       if (sign) {if (sign<0) xrho = -xrho;}  // sign is imposed
 295       else { // determine automatically the sign from direction
 296         if (dir>0) xrho = -xrho; // outward should be negative
 297       }
 298       //
 299       if (!track->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE;
 300     }
 301
 302     if (rotateTo){
 303       track->GetXYZ(xyz1);   // global position
 304       Double_t alphan = TMath::ATan2(xyz1[1], xyz1[0]);
 305       if (maxSnp>0) {
 306         if (TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE;
 307
 308         //
 309         Double_t ca=TMath::Cos(alphan-track->GetAlpha()), sa=TMath::Sin(alphan-track->GetAlpha());
 310         Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf));
 311         Double_t sinNew =  sf*ca - cf*sa;
 312         if (TMath::Abs(sinNew) >= maxSnp) return kFALSE;
 313
 314       }
 315       if (!track->AliExternalTrackParam::Rotate(alphan)) return kFALSE;
 316
 317     }
 318     xpos = track->GetX();
 319     if (addTimeStep && track->IsStartedTimeIntegral()) {
 320       if (!rotateTo) track->GetXYZ(xyz1); // if rotateTo==kTRUE, then xyz1 is already extracted
 321       Double_t dX=xyz0[0]-xyz1[0],dY=xyz0[1]-xyz1[1],dZ=xyz0[2]-xyz1[2];
 322       Double_t d=TMath::Sqrt(dX*dX + dY*dY + dZ*dZ);
 323       if (sign) {if (sign>0) d = -d;}  // step sign is imposed, positive means inward direction
 324       else { // determine automatically the sign from direction
 325         if (dir<0) d = -d;
 326       }
 327       track->AddTimeStep(d);
 328     }
 329   }
 330   return kTRUE;
 331 }
 332
 333 Int_t AliTrackerBase::PropagateTrackTo2(AliExternalTrackParam *track, Double_t xToGo,
 334                                         Double_t mass, Double_t maxStep, Bool_t rotateTo, Double_t maxSnp, Int_t sign, Bool_t addTimeStep, Bool_t correctMaterialBudget){
 335   //----------------------------------------------------------------
 336   //
 337   // Propagates the track to the plane X=xk (cm) using the magnetic field map
 338   // and correcting for the crossed material.
 339   //
 340   // mass     - mass used in propagation - used for energy loss correction
 341   // maxStep  - maximal step for propagation
 342   //
 343   //  Origin: Marian Ivanov,  Marian.Ivanov@cern.ch
 344   //
 345   //----------------------------------------------------------------
 346   const Double_t kEpsilon = 0.00001;
 347   Double_t xpos     = track->GetX();
 348   Int_t dir         = (xpos<xToGo) ? 1:-1;
 349   //
 350   while ( (xToGo-xpos)*dir > kEpsilon){
 351     Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
 352     Double_t x    = xpos+step;
 353     Double_t xyz0[3],xyz1[3],param[7];
 354     track->GetXYZ(xyz0);   //starting global position
 355
 356     Double_t bz=GetBz(xyz0); // getting the local Bz
 357     if (!track->GetXYZAt(x,bz,xyz1)) return -1;   // no prolongation
 358     xyz1[2]+=kEpsilon; // waiting for bug correction in geo
 359
 360     if (maxSnp>0 && TMath::Abs(track->GetSnpAt(x,bz)) >= maxSnp) return -2;
 361     if (!track->PropagateTo(x,bz))  return -3;
 362
 363     if (correctMaterialBudget){
 364       MeanMaterialBudget(xyz0,xyz1,param);
 365       Double_t xrho=param[0]*param[4], xx0=param[1];
 366       if (sign) {if (sign<0) xrho = -xrho;}  // sign is imposed
 367       else { // determine automatically the sign from direction
 368         if (dir>0) xrho = -xrho; // outward should be negative
 369       }
 370       //
 371       if (!track->CorrectForMeanMaterial(xx0,xrho,mass)) return -4;
 372     }
 373
 374     if (rotateTo){
 375       track->GetXYZ(xyz1);   // global position
 376       Double_t alphan = TMath::ATan2(xyz1[1], xyz1[0]);
 377       if (maxSnp>0) {
 378         if (TMath::Abs(track->GetSnp()) >= maxSnp) return -5;
 379
 380         //
 381         Double_t ca=TMath::Cos(alphan-track->GetAlpha()), sa=TMath::Sin(alphan-track->GetAlpha());
 382         Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf));
 383         Double_t sinNew =  sf*ca - cf*sa;
 384         if (TMath::Abs(sinNew) >= maxSnp) return -6;
 385
 386       }
 387       if (!track->AliExternalTrackParam::Rotate(alphan)) return -7;
 388
 389     }
 390     xpos = track->GetX();
 391     if (addTimeStep && track->IsStartedTimeIntegral()) {
 392       if (!rotateTo) track->GetXYZ(xyz1); // if rotateTo==kTRUE, then xyz1 is already extracted
 393       Double_t dX=xyz0[0]-xyz1[0],dY=xyz0[1]-xyz1[1],dZ=xyz0[2]-xyz1[2];
 394       Double_t d=TMath::Sqrt(dX*dX + dY*dY + dZ*dZ);
 395       if (sign) {if (sign>0) d = -d;}  // step sign is imposed, positive means inward direction
 396       else { // determine automatically the sign from direction
 397         if (dir<0) d = -d;
 398       }
 399       track->AddTimeStep(d);
 400     }
 401   }
 402   return 1;
 403 }
 404
 405 Bool_t
 406 AliTrackerBase::PropagateTrackToBxByBz(AliExternalTrackParam *track,
 407                                        Double_t xToGo,Double_t mass, Double_t maxStep, Bool_t rotateTo, Double_t maxSnp,Int_t sign, Bool_t addTimeStep){
 408   //----------------------------------------------------------------
 409   //
 410   // Propagates the track to the plane X=xk (cm)
 411   // taking into account all the three components of the magnetic field
 412   // and correcting for the crossed material.
 413   //
 414   // mass     - mass used in propagation - used for energy loss correction
 415   // maxStep  - maximal step for propagation
 416   //
 417   //  Origin: Marian Ivanov,  Marian.Ivanov@cern.ch
 418   //
 419   //----------------------------------------------------------------
 420   const Double_t kEpsilon = 0.00001;
 421   Double_t xpos     = track->GetX();
 422   Int_t dir         = (xpos<xToGo) ? 1:-1;
 423   //
 424   while ( (xToGo-xpos)*dir > kEpsilon){
 425     Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
 426     Double_t x    = xpos+step;
 427     Double_t xyz0[3],xyz1[3],param[7];
 428     track->GetXYZ(xyz0);   //starting global position
 429
 430     Double_t b[3]; GetBxByBz(xyz0,b); // getting the local Bx, By and Bz
 431
 432     if (!track->GetXYZAt(x,b[2],xyz1)) return kFALSE;   // no prolongation
 433     xyz1[2]+=kEpsilon; // waiting for bug correction in geo
 434
 435     if (maxSnp>0 && TMath::Abs(track->GetSnpAt(x,b[2])) >= maxSnp) return kFALSE;
 436     if (!track->PropagateToBxByBz(x,b))  return kFALSE;
 437
 438     MeanMaterialBudget(xyz0,xyz1,param);
 439     Double_t xrho=param[0]*param[4], xx0=param[1];
 440     if (sign) {if (sign<0) xrho = -xrho;}  // sign is imposed
 441     else { // determine automatically the sign from direction
 442       if (dir>0) xrho = -xrho; // outward should be negative
 443     }
 444     //
 445     if (!track->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE;
 446     if (rotateTo){
 447       track->GetXYZ(xyz1);   // global position
 448       Double_t alphan = TMath::ATan2(xyz1[1], xyz1[0]);
 449       if (maxSnp>0) {
 450         if (TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE;
 451         Double_t ca=TMath::Cos(alphan-track->GetAlpha()), sa=TMath::Sin(alphan-track->GetAlpha());
 452         Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf));
 453         Double_t sinNew =  sf*ca - cf*sa;
 454         if (TMath::Abs(sinNew) >= maxSnp) return kFALSE;
 455       }
 456       if (!track->AliExternalTrackParam::Rotate(alphan)) return kFALSE;
 457     }
 458     xpos = track->GetX();
 459     if (addTimeStep && track->IsStartedTimeIntegral()) {
 460       if (!rotateTo) track->GetXYZ(xyz1); // if rotateTo==kTRUE, then xyz1 is already extracted
 461       Double_t dX=xyz0[0]-xyz1[0],dY=xyz0[1]-xyz1[1],dZ=xyz0[2]-xyz1[2];
 462       Double_t d=TMath::Sqrt(dX*dX + dY*dY + dZ*dZ);
 463       if (sign) {if (sign>0) d = -d;}  // step sign is imposed, positive means inward direction
 464       else { // determine automatically the sign from direction
 465         if (dir<0) d = -d;
 466       }
 467       track->AddTimeStep(d);
 468     }
 469   }
 470   return kTRUE;
 471 }
 472
 473 Double_t AliTrackerBase::GetTrackPredictedChi2(AliExternalTrackParam *track,
 474                                            Double_t mass, Double_t step,
 475                                      const AliExternalTrackParam *backup) {
 476   //
 477   // This function brings the "track" with particle "mass" [GeV]
 478   // to the same local coord. system and the same reference plane as
 479   // of the "backup", doing it in "steps" [cm].
 480   // Then, it calculates the 5D predicted Chi2 for these two tracks
 481   //
 482   Double_t chi2=kVeryBig;
 483   Double_t alpha=backup->GetAlpha();
 484   if (!track->Rotate(alpha)) return chi2;
 485
 486   Double_t xb=backup->GetX();
 487   Double_t sign=(xb < track->GetX()) ? 1. : -1.;
 488   if (!PropagateTrackTo(track,xb,mass,step,kFALSE,kAlmost1,sign)) return chi2;
 489
 490   chi2=track->GetPredictedChi2(backup);
 491
 492   return chi2;
 493 }
 494
 495
 496
 497
 498 Double_t AliTrackerBase::MakeC(Double_t x1,Double_t y1,
 499                    Double_t x2,Double_t y2,
 500                    Double_t x3,Double_t y3)
 501 {
 502   //-----------------------------------------------------------------
 503   // Initial approximation of the track curvature
 504   //-----------------------------------------------------------------
 505   x3 -=x1;
 506   x2 -=x1;
 507   y3 -=y1;
 508   y2 -=y1;
 509   //
 510   Double_t det = x3*y2-x2*y3;
 511   if (TMath::Abs(det)<1e-10) {
 512     return 0;
 513   }
 514   //
 515   Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
 516   Double_t x0 = x3*0.5-y3*u;
 517   Double_t y0 = y3*0.5+x3*u;
 518   Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
 519   if (det>0) c2*=-1;
 520   return c2;
 521 }
 522
 523 Double_t AliTrackerBase::MakeSnp(Double_t x1,Double_t y1,
 524                    Double_t x2,Double_t y2,
 525                    Double_t x3,Double_t y3)
 526 {
 527   //-----------------------------------------------------------------
 528   // Initial approximation of the track snp
 529   //-----------------------------------------------------------------
 530   x3 -=x1;
 531   x2 -=x1;
 532   y3 -=y1;
 533   y2 -=y1;
 534   //
 535   Double_t det = x3*y2-x2*y3;
 536   if (TMath::Abs(det)<1e-10) {
 537     return 0;
 538   }
 539   //
 540   Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
 541   Double_t x0 = x3*0.5-y3*u;
 542   Double_t y0 = y3*0.5+x3*u;
 543   Double_t c2 = 1./TMath::Sqrt(x0*x0+y0*y0);
 544   x0*=c2;
 545   x0=TMath::Abs(x0);
 546   if (y2*x2<0.) x0*=-1;
 547   return x0;
 548 }
 549
 550 Double_t AliTrackerBase::MakeTgl(Double_t x1,Double_t y1,
 551                    Double_t x2,Double_t y2,
 552                    Double_t z1,Double_t z2, Double_t c)
 553 {
 554   //-----------------------------------------------------------------
 555   // Initial approximation of the tangent of the track dip angle
 556   //-----------------------------------------------------------------
 557   //
 558   x2-=x1;
 559   y2-=y1;
 560   z2-=z1;
 561   Double_t d  =  TMath::Sqrt(x2*x2+y2*y2);  // distance  straight line
 562   if (TMath::Abs(d*c*0.5)>1) return 0;
 563   Double_t   angle2    = TMath::ASin(d*c*0.5);
 564   angle2  = z2*TMath::Abs(c/(angle2*2.));
 565   return angle2;
 566 }
 567
 568
 569 Double_t AliTrackerBase::MakeTgl(Double_t x1,Double_t y1,
 570                    Double_t x2,Double_t y2,
 571                    Double_t z1,Double_t z2)
 572 {
 573   //-----------------------------------------------------------------
 574   // Initial approximation of the tangent of the track dip angle
 575   //-----------------------------------------------------------------
 576   return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
 577 }
 578
 579
 580 AliExternalTrackParam * AliTrackerBase::MakeSeed( AliTrackPoint &point0, AliTrackPoint &point1, AliTrackPoint &point2){
 581   //
 582   // Make Seed  - AliExternalTrackParam from input 3 points
 583   // returning seed in local frame of point0
 584   //
 585   Double_t xyz0[3]={0,0,0};
 586   Double_t xyz1[3]={0,0,0};
 587   Double_t xyz2[3]={0,0,0};
 588   Double_t alpha=point0.GetAngle();
 589   Double_t xyz[3]={point0.GetX(),point0.GetY(),point0.GetZ()};
 590   Double_t bxyz[3]; GetBxByBz(xyz,bxyz);
 591   Double_t bz = bxyz[2];
 592   //
 593   // get points in frame of point 0
 594   //
 595   AliTrackPoint p0r = point0.Rotate(alpha);
 596   AliTrackPoint p1r = point1.Rotate(alpha);
 597   AliTrackPoint p2r = point2.Rotate(alpha);
 598   xyz0[0]=p0r.GetX();
 599   xyz0[1]=p0r.GetY();
 600   xyz0[2]=p0r.GetZ();
 601   xyz1[0]=p1r.GetX();
 602   xyz1[1]=p1r.GetY();
 603   xyz1[2]=p1r.GetZ();
 604   xyz2[0]=p2r.GetX();
 605   xyz2[1]=p2r.GetY();
 606   xyz2[2]=p2r.GetZ();
 607   //
 608   // make covariance estimate
 609   //
 610   Double_t covar[15];
 611   Double_t param[5]={0,0,0,0,0};
 612   for (Int_t m=0; m<15; m++) covar[m]=0;
 613   //
 614   // calculate intitial param
 615   param[0]=xyz0[1];
 616   param[1]=xyz0[2];
 617   param[2]=MakeSnp(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]);
 618   param[4]=MakeC(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]);
 619   param[3]=MakeTgl(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz0[2],xyz1[2],param[4]);
 620
 621   //covariance matrix - only diagonal elements
 622   //Double_t dist=p0r.GetX()-p2r.GetX();
 623   Double_t deltaP=0;
 624   covar[0]= p0r.GetCov()[3];
 625   covar[2]= p0r.GetCov()[5];
 626   //sigma snp
 627   deltaP= (MakeSnp(xyz0[0],xyz0[1]+TMath::Sqrt(p0r.GetCov()[3]),xyz1[0],xyz1[1],xyz2[0],xyz2[1])-param[2]);
 628   covar[5]+= deltaP*deltaP;
 629   deltaP= (MakeSnp(xyz0[0],xyz0[1],xyz1[0],xyz1[1]+TMath::Sqrt(p1r.GetCov()[3]),xyz2[0],xyz2[1])-param[2]);
 630   covar[5]+= deltaP*deltaP;
 631   deltaP= (MakeSnp(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]+TMath::Sqrt(p1r.GetCov()[3]))-param[2]);
 632   covar[5]+= deltaP*deltaP;
 633   //sigma tgl
 634   //
 635   deltaP=MakeTgl(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz0[2]+TMath::Sqrt(p1r.GetCov()[5]),xyz1[2],param[4])-param[3];
 636   covar[9]+= deltaP*deltaP;
 637   deltaP=MakeTgl(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz0[2],xyz1[2]+TMath::Sqrt(p1r.GetCov()[5]),param[4])-param[3];
 638   covar[9]+= deltaP*deltaP;
 639   //
 640
 641   deltaP=MakeC(xyz0[0],xyz0[1]+TMath::Sqrt(p0r.GetCov()[3]),xyz1[0],xyz1[1],xyz2[0],xyz2[1])-param[4];
 642   covar[14]+= deltaP*deltaP;
 643   deltaP=MakeC(xyz0[0],xyz0[1],xyz1[0],xyz1[1]+TMath::Sqrt(p1r.GetCov()[3]),xyz2[0],xyz2[1])-param[4];
 644   covar[14]+= deltaP*deltaP;
 645   deltaP=MakeC(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]+TMath::Sqrt(p2r.GetCov()[3]))-param[4];
 646   covar[14]+= deltaP*deltaP;
 647
 648   covar[14]/=(bz*kB2C)*(bz*kB2C);
 649   param[4]/=(bz*kB2C); // transform to 1/pt
 650   AliExternalTrackParam * trackParam = new AliExternalTrackParam(xyz0[0],alpha,param, covar);
 651   if (0) {
 652     // consistency check  -to put warnings here
 653     // small disagrement once Track extrapolation used
 654     // nice agreement in seeds with MC track parameters - problem in extrapoloation - to be fixed
 655     // to check later
 656     Double_t y1,y2,z1,z2;
 657     trackParam->GetYAt(xyz1[0],bz,y1);
 658     trackParam->GetZAt(xyz1[0],bz,z1);
 659     trackParam->GetYAt(xyz2[0],bz,y2);
 660     trackParam->GetZAt(xyz2[0],bz,z2);
 661     if (TMath::Abs(y1-xyz1[1])> TMath::Sqrt(p1r.GetCov()[3]*5)){
 662       AliWarningClass("Seeding problem y1\n");
 663     }
 664     if (TMath::Abs(y2-xyz2[1])> TMath::Sqrt(p2r.GetCov()[3]*5)){
 665       AliWarningClass("Seeding problem y2\n");
 666     }
 667     if (TMath::Abs(z1-xyz1[2])> TMath::Sqrt(p1r.GetCov()[5]*5)){
 668       AliWarningClass("Seeding problem z1\n");
 669     }
 670   }
 671   return trackParam;
 672 }
 673
 674 Double_t  AliTrackerBase::FitTrack(AliExternalTrackParam * trackParam, AliTrackPointArray *pointArray, Double_t mass, Double_t maxStep){
 675   //
 676   // refit the track  - trackParam using the points in point array
 677   //
 678   const Double_t kMaxSnp=0.99;
 679   if (!trackParam) return 0;
 680   Int_t  npoints=pointArray->GetNPoints();
 681   AliTrackPoint point,point2;
 682   Double_t pointPos[2]={0,0};
 683   Double_t pointCov[3]={0,0,0};
 684   // choose coordinate frame
 685   // in standard way the coordinate frame should be changed point by point
 686   // Some problems with rotation observed
 687   // rotate method of AliExternalTrackParam should be revisited
 688   pointArray->GetPoint(point,0);
 689   pointArray->GetPoint(point2,npoints-1);
 690   Double_t alpha=TMath::ATan2(point.GetY()-point2.GetY(), point.GetX()-point2.GetX());
 691
 692   for (Int_t ipoint=npoints-1; ipoint>0; ipoint-=1){
 693     pointArray->GetPoint(point,ipoint);
 694     AliTrackPoint pr = point.Rotate(alpha);
 695     trackParam->Rotate(alpha);
 696     Bool_t status = PropagateTrackTo(trackParam,pr.GetX(),mass,maxStep,kFALSE,kMaxSnp);
 697     if(!status){
 698       AliWarningClass("Problem to propagate\n");
 699       break;
 700     }
 701     if (TMath::Abs(trackParam->GetSnp())>kMaxSnp){
 702       AliWarningClass("sin(phi) > kMaxSnp \n");
 703       break;
 704     }
 705     pointPos[0]=pr.GetY();//local y
 706     pointPos[1]=pr.GetZ();//local z
 707     pointCov[0]=pr.GetCov()[3];//simay^2
 708     pointCov[1]=pr.GetCov()[4];//sigmayz
 709     pointCov[2]=pr.GetCov()[5];//sigmaz^2
 710     trackParam->Update(pointPos,pointCov);
 711   }
 712   return 0;
 713 }
 714
 715
 716
 717 void AliTrackerBase::UpdateTrack(AliExternalTrackParam &track1, const AliExternalTrackParam &track2){
 718   //
 719   // Update track 1 with track 2
 720   //
 721   //
 722   //
 723   TMatrixD vecXk(5,1);    // X vector
 724   TMatrixD covXk(5,5);    // X covariance
 725   TMatrixD matHk(5,5);    // vector to mesurement
 726   TMatrixD measR(5,5);    // measurement error
 727   TMatrixD vecZk(5,1);    // measurement
 728   //
 729   TMatrixD vecYk(5,1);    // Innovation or measurement residual
 730   TMatrixD matHkT(5,5);
 731   TMatrixD matSk(5,5);    // Innovation (or residual) covariance
 732   TMatrixD matKk(5,5);    // Optimal Kalman gain
 733   TMatrixD mat1(5,5);     // update covariance matrix
 734   TMatrixD covXk2(5,5);   //
 735   TMatrixD covOut(5,5);
 736   //
 737   Double_t *param1=(Double_t*) track1.GetParameter();
 738   Double_t *covar1=(Double_t*) track1.GetCovariance();
 739   Double_t *param2=(Double_t*) track2.GetParameter();
 740   Double_t *covar2=(Double_t*) track2.GetCovariance();
 741   //
 742   // copy data to the matrix
 743   for (Int_t ipar=0; ipar<5; ipar++){
 744     for (Int_t jpar=0; jpar<5; jpar++){
 745       covXk(ipar,jpar) = covar1[track1.GetIndex(ipar, jpar)];
 746       measR(ipar,jpar) = covar2[track2.GetIndex(ipar, jpar)];
 747       matHk(ipar,jpar)=0;
 748       mat1(ipar,jpar)=0;
 749     }
 750     vecXk(ipar,0) = param1[ipar];
 751     vecZk(ipar,0) = param2[ipar];
 752     matHk(ipar,ipar)=1;
 753     mat1(ipar,ipar)=1;
 754   }
 755   //
 756   //
 757   //
 758   //
 759   //
 760   vecYk = vecZk-matHk*vecXk;                 // Innovation or measurement residual
 761   matHkT=matHk.T(); matHk.T();
 762   matSk = (matHk*(covXk*matHkT))+measR;      // Innovation (or residual) covariance
 763   matSk.Invert();
 764   matKk = (covXk*matHkT)*matSk;              //  Optimal Kalman gain
 765   vecXk += matKk*vecYk;                      //  updated vector
 766   covXk2 = (mat1-(matKk*matHk));
 767   covOut =  covXk2*covXk;
 768   //
 769   //
 770   //
 771   // copy from matrix to parameters
 772   if (0) {
 773     vecXk.Print();
 774     vecZk.Print();
 775     //
 776     measR.Print();
 777     covXk.Print();
 778     covOut.Print();
 779     //
 780     track1.Print();
 781     track2.Print();
 782   }
 783
 784   for (Int_t ipar=0; ipar<5; ipar++){
 785     param1[ipar]= vecXk(ipar,0) ;
 786     for (Int_t jpar=0; jpar<5; jpar++){
 787       covar1[track1.GetIndex(ipar, jpar)]=covOut(ipar,jpar);
 788     }
 789   }
 790 }
 791
 792