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