STEER/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 AliTPCtracker, AliITStrackerV2 and AliTRDtracker
  21 //        Origin: Iouri Belikov, CERN, Jouri.Belikov@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
  31 extern TGeoManager *gGeoManager;
  32
  33 ClassImp(AliTrackerBase)
  34
  35 AliTrackerBase::AliTrackerBase():
  36   TObject(),
  37   fX(0),
  38   fY(0),
  39   fZ(0),
  40   fSigmaX(0.005),
  41   fSigmaY(0.005),
  42   fSigmaZ(0.010)
  43 {
  44   //--------------------------------------------------------------------
  45   // The default constructor.
  46   //--------------------------------------------------------------------
  47   if (!TGeoGlobalMagField::Instance()->GetField())
  48     AliWarning("Field map is not set.");
  49 }
  50
  51 //__________________________________________________________________________
  52 AliTrackerBase::AliTrackerBase(const AliTrackerBase &atr):
  53   TObject(atr),
  54   fX(atr.fX),
  55   fY(atr.fY),
  56   fZ(atr.fZ),
  57   fSigmaX(atr.fSigmaX),
  58   fSigmaY(atr.fSigmaY),
  59   fSigmaZ(atr.fSigmaZ)
  60 {
  61   //--------------------------------------------------------------------
  62   // The default constructor.
  63   //--------------------------------------------------------------------
  64   if (!TGeoGlobalMagField::Instance()->GetField())
  65     AliWarning("Field map is not set.");
  66 }
  67
  68 //__________________________________________________________________________
  69 Double_t AliTrackerBase::GetBz()
  70 {
  71   AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
  72   if (!fld) return 0.5*kAlmost0Field;
  73   Double_t bz = fld->SolenoidField();
  74   return TMath::Sign(0.5*kAlmost0Field,bz) + bz;
  75 }
  76
  77 //__________________________________________________________________________
  78 Double_t AliTrackerBase::GetBz(const Double_t *r) {
  79   //------------------------------------------------------------------
  80   // Returns Bz (kG) at the point "r" .
  81   //------------------------------------------------------------------
  82   AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
  83   if (!fld) return  0.5*kAlmost0Field;
  84   Double_t bz = fld->GetBz(r);
  85   return  TMath::Sign(0.5*kAlmost0Field,bz) + bz;
  86 }
  87
  88 //__________________________________________________________________________
  89 void AliTrackerBase::GetBxByBz(const Double_t r[3], Double_t b[3]) {
  90   //------------------------------------------------------------------
  91   // Returns Bx, By and Bz (kG) at the point "r" .
  92   //------------------------------------------------------------------
  93   AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
  94   if (!fld) {
  95      b[0] = b[1] = 0.;
  96      b[2] = 0.5*kAlmost0Field;
  97      return;
  98   }
  99
 100   if (fld->IsUniform()) {
 101      b[0] = b[1] = 0.;
 102      b[2] = fld->SolenoidField();
 103   }  else {
 104      fld->Field(r,b);
 105   }
 106   b[2] = (TMath::Sign(0.5*kAlmost0Field,b[2]) + b[2]);
 107   return;
 108 }
 109
 110 Double_t AliTrackerBase::MeanMaterialBudget(const Double_t *start, const Double_t *end, Double_t *mparam)
 111 {
 112   //
 113   // Calculate mean material budget and material properties between
 114   //    the points "start" and "end".
 115   //
 116   // "mparam" - parameters used for the energy and multiple scattering
 117   //  corrections:
 118   //
 119   // mparam[0] - mean density: sum(x_i*rho_i)/sum(x_i) [g/cm3]
 120   // mparam[1] - equivalent rad length fraction: sum(x_i/X0_i) [adimensional]
 121   // mparam[2] - mean A: sum(x_i*A_i)/sum(x_i) [adimensional]
 122   // mparam[3] - mean Z: sum(x_i*Z_i)/sum(x_i) [adimensional]
 123   // mparam[4] - length: sum(x_i) [cm]
 124   // mparam[5] - Z/A mean: sum(x_i*Z_i/A_i)/sum(x_i) [adimensional]
 125   // mparam[6] - number of boundary crosses
 126   //
 127   //  Origin:  Marian Ivanov, Marian.Ivanov@cern.ch
 128   //
 129   //  Corrections and improvements by
 130   //        Andrea Dainese, Andrea.Dainese@lnl.infn.it,
 131   //        Andrei Gheata,  Andrei.Gheata@cern.ch
 132   //
 133
 134   mparam[0]=0; mparam[1]=1; mparam[2] =0; mparam[3] =0;
 135   mparam[4]=0; mparam[5]=0; mparam[6]=0;
 136   //
 137   Double_t bparam[6]; // total parameters
 138   Double_t lparam[6]; // local parameters
 139
 140   for (Int_t i=0;i<6;i++) bparam[i]=0;
 141
 142   if (!gGeoManager) {
 143     AliErrorClass("No TGeo\n");
 144     return 0.;
 145   }
 146   //
 147   Double_t length;
 148   Double_t dir[3];
 149   length = TMath::Sqrt((end[0]-start[0])*(end[0]-start[0])+
 150                        (end[1]-start[1])*(end[1]-start[1])+
 151                        (end[2]-start[2])*(end[2]-start[2]));
 152   mparam[4]=length;
 153   if (length<TGeoShape::Tolerance()) return 0.0;
 154   Double_t invlen = 1./length;
 155   dir[0] = (end[0]-start[0])*invlen;
 156   dir[1] = (end[1]-start[1])*invlen;
 157   dir[2] = (end[2]-start[2])*invlen;
 158
 159   // Initialize start point and direction
 160   TGeoNode *currentnode = 0;
 161   TGeoNode *startnode = gGeoManager->InitTrack(start, dir);
 162   if (!startnode) {
 163     AliErrorClass(Form("start point out of geometry: x %f, y %f, z %f",
 164                   start[0],start[1],start[2]));
 165     return 0.0;
 166   }
 167   TGeoMaterial *material = startnode->GetVolume()->GetMedium()->GetMaterial();
 168   lparam[0]   = material->GetDensity();
 169   lparam[1]   = material->GetRadLen();
 170   lparam[2]   = material->GetA();
 171   lparam[3]   = material->GetZ();
 172   lparam[4]   = length;
 173   lparam[5]   = lparam[3]/lparam[2];
 174   if (material->IsMixture()) {
 175     TGeoMixture * mixture = (TGeoMixture*)material;
 176     lparam[5] =0;
 177     Double_t sum =0;
 178     for (Int_t iel=0;iel<mixture->GetNelements();iel++){
 179       sum  += mixture->GetWmixt()[iel];
 180       lparam[5]+= mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel];
 181     }
 182     lparam[5]/=sum;
 183   }
 184
 185   // Locate next boundary within length without computing safety.
 186   // Propagate either with length (if no boundary found) or just cross boundary
 187   gGeoManager->FindNextBoundaryAndStep(length, kFALSE);
 188   Double_t step = 0.0; // Step made
 189   Double_t snext = gGeoManager->GetStep();
 190   // If no boundary within proposed length, return current density
 191   if (!gGeoManager->IsOnBoundary()) {
 192     mparam[0] = lparam[0];
 193     mparam[1] = lparam[4]/lparam[1];
 194     mparam[2] = lparam[2];
 195     mparam[3] = lparam[3];
 196     mparam[4] = lparam[4];
 197     return lparam[0];
 198   }
 199   // Try to cross the boundary and see what is next
 200   Int_t nzero = 0;
 201   while (length>TGeoShape::Tolerance()) {
 202     currentnode = gGeoManager->GetCurrentNode();
 203     if (snext<2.*TGeoShape::Tolerance()) nzero++;
 204     else nzero = 0;
 205     if (nzero>3) {
 206       // This means navigation has problems on one boundary
 207       // Try to cross by making a small step
 208       AliErrorClass("Cannot cross boundary\n");
 209       mparam[0] = bparam[0]/step;
 210       mparam[1] = bparam[1];
 211       mparam[2] = bparam[2]/step;
 212       mparam[3] = bparam[3]/step;
 213       mparam[5] = bparam[5]/step;
 214       mparam[4] = step;
 215       mparam[0] = 0.;             // if crash of navigation take mean density 0
 216       mparam[1] = 1000000;        // and infinite rad length
 217       return bparam[0]/step;
 218     }
 219     mparam[6]+=1.;
 220     step += snext;
 221     bparam[1]    += snext/lparam[1];
 222     bparam[2]    += snext*lparam[2];
 223     bparam[3]    += snext*lparam[3];
 224     bparam[5]    += snext*lparam[5];
 225     bparam[0]    += snext*lparam[0];
 226
 227     if (snext>=length) break;
 228     if (!currentnode) break;
 229     length -= snext;
 230     material = currentnode->GetVolume()->GetMedium()->GetMaterial();
 231     lparam[0] = material->GetDensity();
 232     lparam[1]  = material->GetRadLen();
 233     lparam[2]  = material->GetA();
 234     lparam[3]  = material->GetZ();
 235     lparam[5]   = lparam[3]/lparam[2];
 236     if (material->IsMixture()) {
 237       TGeoMixture * mixture = (TGeoMixture*)material;
 238       lparam[5]=0;
 239       Double_t sum =0;
 240       for (Int_t iel=0;iel<mixture->GetNelements();iel++){
 241         sum+= mixture->GetWmixt()[iel];
 242         lparam[5]+= mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel];
 243       }
 244       lparam[5]/=sum;
 245     }
 246     gGeoManager->FindNextBoundaryAndStep(length, kFALSE);
 247     snext = gGeoManager->GetStep();
 248   }
 249   mparam[0] = bparam[0]/step;
 250   mparam[1] = bparam[1];
 251   mparam[2] = bparam[2]/step;
 252   mparam[3] = bparam[3]/step;
 253   mparam[5] = bparam[5]/step;
 254   return bparam[0]/step;
 255 }
 256
 257
 258 Bool_t
 259 AliTrackerBase::PropagateTrackTo(AliExternalTrackParam *track, Double_t xToGo,
 260                              Double_t mass, Double_t maxStep, Bool_t rotateTo, Double_t maxSnp, Double_t sign){
 261   //----------------------------------------------------------------
 262   //
 263   // Propagates the track to the plane X=xk (cm) using the magnetic field map
 264   // and correcting for the crossed material.
 265   //
 266   // mass     - mass used in propagation - used for energy loss correction
 267   // maxStep  - maximal step for propagation
 268   //
 269   //  Origin: Marian Ivanov,  Marian.Ivanov@cern.ch
 270   //
 271   //----------------------------------------------------------------
 272   const Double_t kEpsilon = 0.00001;
 273   Double_t xpos     = track->GetX();
 274   Double_t dir      = (xpos<xToGo) ? 1.:-1.;
 275   //
 276   while ( (xToGo-xpos)*dir > kEpsilon){
 277     Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
 278     Double_t x    = xpos+step;
 279     Double_t xyz0[3],xyz1[3],param[7];
 280     track->GetXYZ(xyz0);   //starting global position
 281
 282     Double_t bz=GetBz(xyz0); // getting the local Bz
 283
 284     if (!track->GetXYZAt(x,bz,xyz1)) return kFALSE;   // no prolongation
 285     xyz1[2]+=kEpsilon; // waiting for bug correction in geo
 286
 287     if (TMath::Abs(track->GetSnpAt(x,bz)) >= maxSnp) return kFALSE;
 288     if (!track->PropagateTo(x,bz))  return kFALSE;
 289
 290     MeanMaterialBudget(xyz0,xyz1,param);
 291     Double_t xrho=param[0]*param[4]*sign, xx0=param[1];
 292
 293     if (!track->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE;
 294     if (rotateTo){
 295       if (TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE;
 296       track->GetXYZ(xyz0);   // global position
 297       Double_t alphan = TMath::ATan2(xyz0[1], xyz0[0]);
 298       //
 299       Double_t ca=TMath::Cos(alphan-track->GetAlpha()),
 300                sa=TMath::Sin(alphan-track->GetAlpha());
 301       Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf));
 302       Double_t sinNew =  sf*ca - cf*sa;
 303       if (TMath::Abs(sinNew) >= maxSnp) return kFALSE;
 304       if (!track->Rotate(alphan)) return kFALSE;
 305     }
 306     xpos = track->GetX();
 307   }
 308   return kTRUE;
 309 }
 310
 311 Bool_t
 312 AliTrackerBase::PropagateTrackToBxByBz(AliExternalTrackParam *track,
 313 Double_t xToGo,
 314                                    Double_t mass, Double_t maxStep, Bool_t rotateTo, Double_t maxSnp,Double_t sign){
 315   //----------------------------------------------------------------
 316   //
 317   // Propagates the track to the plane X=xk (cm)
 318   // taking into account all the three components of the magnetic field
 319   // and correcting for the crossed material.
 320   //
 321   // mass     - mass used in propagation - used for energy loss correction
 322   // maxStep  - maximal step for propagation
 323   //
 324   //  Origin: Marian Ivanov,  Marian.Ivanov@cern.ch
 325   //
 326   //----------------------------------------------------------------
 327   const Double_t kEpsilon = 0.00001;
 328   Double_t xpos     = track->GetX();
 329   Double_t dir      = (xpos<xToGo) ? 1.:-1.;
 330   //
 331   while ( (xToGo-xpos)*dir > kEpsilon){
 332     Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
 333     Double_t x    = xpos+step;
 334     Double_t xyz0[3],xyz1[3],param[7];
 335     track->GetXYZ(xyz0);   //starting global position
 336
 337     Double_t b[3]; GetBxByBz(xyz0,b); // getting the local Bx, By and Bz
 338
 339     if (!track->GetXYZAt(x,b[2],xyz1)) return kFALSE;   // no prolongation
 340     xyz1[2]+=kEpsilon; // waiting for bug correction in geo
 341
 342     if (TMath::Abs(track->GetSnpAt(x,b[2])) >= maxSnp) return kFALSE;
 343     if (!track->PropagateToBxByBz(x,b))  return kFALSE;
 344
 345     MeanMaterialBudget(xyz0,xyz1,param);
 346     Double_t xrho=param[0]*param[4]*sign, xx0=param[1];
 347
 348     if (!track->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE;
 349     if (rotateTo){
 350       if (TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE;
 351       track->GetXYZ(xyz0);   // global position
 352       Double_t alphan = TMath::ATan2(xyz0[1], xyz0[0]);
 353       //
 354       Double_t ca=TMath::Cos(alphan-track->GetAlpha()),
 355                sa=TMath::Sin(alphan-track->GetAlpha());
 356       Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf));
 357       Double_t sinNew =  sf*ca - cf*sa;
 358       if (TMath::Abs(sinNew) >= maxSnp) return kFALSE;
 359       if (!track->Rotate(alphan)) return kFALSE;
 360     }
 361     xpos = track->GetX();
 362   }
 363   return kTRUE;
 364 }
 365
 366 Double_t AliTrackerBase::GetTrackPredictedChi2(AliExternalTrackParam *track,
 367                                            Double_t mass, Double_t step,
 368                                      const AliExternalTrackParam *backup) {
 369   //
 370   // This function brings the "track" with particle "mass" [GeV]
 371   // to the same local coord. system and the same reference plane as
 372   // of the "backup", doing it in "steps" [cm].
 373   // Then, it calculates the 5D predicted Chi2 for these two tracks
 374   //
 375   Double_t chi2=kVeryBig;
 376   Double_t alpha=backup->GetAlpha();
 377   if (!track->Rotate(alpha)) return chi2;
 378
 379   Double_t xb=backup->GetX();
 380   Double_t sign=(xb < track->GetX()) ? 1. : -1.;
 381   if (!PropagateTrackTo(track,xb,mass,step,kFALSE,kAlmost1,sign)) return chi2;
 382
 383   chi2=track->GetPredictedChi2(backup);
 384
 385   return chi2;
 386 }