STEER/AliKalmanTrack.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 //                Implementation of the AliKalmanTrack class
  20 //   that is the base for AliTPCtrack, AliITStrackV2 and AliTRDtrack
  21 //        Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
  22 //-------------------------------------------------------------------------
  23 #include "AliTracker.h"
  24 #include "AliKalmanTrack.h"
  25 #include "TGeoManager.h"
  26
  27 ClassImp(AliKalmanTrack)
  28
  29 //_______________________________________________________________________
  30 AliKalmanTrack::AliKalmanTrack():
  31   fLab(-3141593),
  32   fFakeRatio(0),
  33   fChi2(0),
  34   fMass(AliPID::ParticleMass(AliPID::kPion)),
  35   fN(0),
  36   fLocalConvConst(0),
  37   fStartTimeIntegral(kFALSE),
  38   fIntegratedLength(0)
  39 {
  40   //
  41   // Default constructor
  42   //
  43   if (AliTracker::GetFieldMap()==0) {
  44       AliFatal("The magnetic field has not been set!");
  45   }
  46
  47   for(Int_t i=0; i<AliPID::kSPECIES; i++) fIntegratedTime[i] = 0;
  48 }
  49
  50 //_______________________________________________________________________
  51 AliKalmanTrack::AliKalmanTrack(const AliKalmanTrack &t):
  52   TObject(t),
  53   fLab(t.fLab),
  54   fFakeRatio(t.fFakeRatio),
  55   fChi2(t.fChi2),
  56   fMass(t.fMass),
  57   fN(t.fN),
  58   fLocalConvConst(t.fLocalConvConst),
  59   fStartTimeIntegral(t.fStartTimeIntegral),
  60   fIntegratedLength(t.fIntegratedLength)
  61 {
  62   //
  63   // Copy constructor
  64   //
  65   if (AliTracker::GetFieldMap()==0) {
  66     AliFatal("The magnetic field has not been set!");
  67   }
  68
  69   for (Int_t i=0; i<AliPID::kSPECIES; i++)
  70       fIntegratedTime[i] = t.fIntegratedTime[i];
  71 }
  72
  73 //_______________________________________________________________________
  74 void AliKalmanTrack::StartTimeIntegral()
  75 {
  76   // Sylwester Radomski, GSI
  77   // S.Radomski@gsi.de
  78   //
  79   // Start time integration
  80   // To be called at Vertex by ITS tracker
  81   //
  82
  83   //if (fStartTimeIntegral)
  84   //  AliWarning("Reseting Recorded Time.");
  85
  86   fStartTimeIntegral = kTRUE;
  87   for(Int_t i=0; i<AliPID::kSPECIES; i++) fIntegratedTime[i] = 0;
  88   fIntegratedLength = 0;
  89 }
  90
  91 //_______________________________________________________________________
  92 void AliKalmanTrack:: AddTimeStep(Double_t length)
  93 {
  94   //
  95   // Add step to integrated time
  96   // this method should be called by a sublasses at the end
  97   // of the PropagateTo function or by a tracker
  98   // each time step is made.
  99   //
 100   // If integration not started function does nothing
 101   //
 102   // Formula
 103   // dt = dl * sqrt(p^2 + m^2) / p
 104   // p = pT * (1 + tg^2 (lambda) )
 105   //
 106   // pt = 1/external parameter [4]
 107   // tg lambda = external parameter [3]
 108   //
 109   //
 110   // Sylwester Radomski, GSI
 111   // S.Radomski@gsi.de
 112   //
 113
 114   static const Double_t kcc = 2.99792458e-2;
 115
 116   if (!fStartTimeIntegral) return;
 117
 118   fIntegratedLength += length;
 119
 120   Double_t xr, param[5];
 121   Double_t pt, tgl;
 122
 123   GetExternalParameters(xr, param);
 124   pt =  1/param[4] ;
 125   tgl = param[3];
 126
 127   Double_t p = TMath::Abs(pt * TMath::Sqrt(1+tgl*tgl));
 128
 129   if (length > 100) return;
 130
 131   for (Int_t i=0; i<AliPID::kSPECIES; i++) {
 132
 133     Double_t mass = AliPID::ParticleMass(i);
 134     Double_t correction = TMath::Sqrt( pt*pt * (1 + tgl*tgl) + mass * mass ) / p;
 135     Double_t time = length * correction / kcc;
 136
 137     fIntegratedTime[i] += time;
 138   }
 139 }
 140
 141 //_______________________________________________________________________
 142 Double_t AliKalmanTrack::GetIntegratedTime(Int_t pdg) const
 143 {
 144   // Sylwester Radomski, GSI
 145   // S.Radomski@gsi.de
 146   //
 147   // Return integrated time hypothesis for a given particle
 148   // type assumption.
 149   //
 150   // Input parameter:
 151   // pdg - Pdg code of a particle type
 152   //
 153
 154
 155   if (!fStartTimeIntegral) {
 156     AliWarning("Time integration not started");
 157     return 0.;
 158   }
 159
 160   for (Int_t i=0; i<AliPID::kSPECIES; i++)
 161     if (AliPID::ParticleCode(i) == TMath::Abs(pdg)) return fIntegratedTime[i];
 162
 163   AliWarning(Form("Particle type [%d] not found", pdg));
 164   return 0;
 165 }
 166
 167 void AliKalmanTrack::GetIntegratedTimes(Double_t *times) const {
 168   for (Int_t i=0; i<AliPID::kSPECIES; i++) times[i]=fIntegratedTime[i];
 169 }
 170
 171 void AliKalmanTrack::SetIntegratedTimes(const Double_t *times) {
 172   for (Int_t i=0; i<AliPID::kSPECIES; i++) fIntegratedTime[i]=times[i];
 173 }
 174
 175 void AliKalmanTrack::External2Helix(Double_t helix[6]) const {
 176   //--------------------------------------------------------------------
 177   // External track parameters -> helix parameters
 178   //--------------------------------------------------------------------
 179   Double_t alpha,x,cs,sn;
 180   GetExternalParameters(x,helix); alpha=GetAlpha();
 181
 182   cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
 183   helix[5]=x*cs - helix[0]*sn;            // x0
 184   helix[0]=x*sn + helix[0]*cs;            // y0
 185 //helix[1]=                               // z0
 186   helix[2]=TMath::ASin(helix[2]) + alpha; // phi0
 187 //helix[3]=                               // tgl
 188   helix[4]=helix[4]/GetLocalConvConst();  // C
 189 }
 190
 191 Double_t AliKalmanTrack::MeanMaterialBudget(Double_t *start, Double_t *end, Double_t *mparam)
 192 {
 193   //
 194   // calculate mean material budget and material properties beween point start and end
 195   // mparam - returns parameters used for dEdx and multiple scatering
 196   //
 197   // mparam[0] - density mean
 198   // mparam[1] - rad length
 199   // mparam[2] - A mean
 200   // mparam[3] - Z mean
 201   // mparam[4] - length
 202   // mparam[5] - Z/A mean
 203   // mparam[6] - number of boundary crosses
 204   //
 205     mparam[0]=0; mparam[1]=1; mparam[2] =0; mparam[3] =0, mparam[4]=0, mparam[5]=0; mparam[6]=0;
 206   //
 207   Double_t bparam[6], lparam[6];          // bparam - total param - lparam - local parameters
 208   for (Int_t i=0;i<6;i++) bparam[i]=0;    //
 209
 210   if (!gGeoManager) {
 211     printf("ERROR: no TGeo\n");
 212     return 0.;
 213   }
 214   //
 215   Double_t length;
 216   Double_t dir[3];
 217   length = TMath::Sqrt((end[0]-start[0])*(end[0]-start[0])+
 218                        (end[1]-start[1])*(end[1]-start[1])+
 219                        (end[2]-start[2])*(end[2]-start[2]));
 220   mparam[4]=length;
 221   if (length<TGeoShape::Tolerance()) return 0.0;
 222   Double_t invlen = 1./length;
 223   dir[0] = (end[0]-start[0])*invlen;
 224   dir[1] = (end[1]-start[1])*invlen;
 225   dir[2] = (end[2]-start[2])*invlen;
 226   // Initialize start point and direction
 227   TGeoNode *currentnode = 0;
 228   TGeoNode *startnode = gGeoManager->InitTrack(start, dir);
 229   //  printf("%s length=%f\n",gGeoManager->GetPath(),length);
 230   if (!startnode) {
 231     printf("ERROR: start point out of geometry\n");
 232     return 0.0;
 233   }
 234   TGeoMaterial *material = startnode->GetVolume()->GetMedium()->GetMaterial();
 235   lparam[0] = material->GetDensity();
 236   lparam[1]   = material->GetRadLen();
 237   lparam[2]   = material->GetA();
 238   lparam[3]   = material->GetZ();
 239   lparam[4]   = length;
 240   lparam[5]   = lparam[3]/lparam[2];
 241   if (material->IsMixture()) {
 242     lparam[1]*=lparam[0];  // different normalization in the modeler for mixture
 243     TGeoMixture * mixture = (TGeoMixture*)material;
 244     lparam[5] =0;
 245     Double_t sum =0;
 246     for (Int_t iel=0;iel<mixture->GetNelements();iel++){
 247       sum  += mixture->GetWmixt()[iel];
 248       lparam[5]+= mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel];
 249     }
 250     lparam[5]/=sum;
 251   }
 252   gGeoManager->FindNextBoundary(length);
 253   Double_t snext = gGeoManager->GetStep();
 254   Double_t step = 0.0;
 255   // If no boundary within proposed length, return current density
 256   if (snext>=length) {
 257     for (Int_t ip=0;ip<5;ip++) mparam[ip] = lparam[ip];
 258     return lparam[0];
 259   }
 260   // Try to cross the boundary and see what is next
 261   while (length>TGeoShape::Tolerance()) {
 262     mparam[6]+=1.;
 263     currentnode = gGeoManager->Step();
 264     step += snext+1.E-6;
 265     bparam[1]    += snext*lparam[1];
 266     bparam[2]    += snext*lparam[2];
 267     bparam[3]    += snext*lparam[3];
 268     bparam[5]    += snext*lparam[5];
 269     bparam[0]    += snext*lparam[0];
 270
 271     if (snext>=length) break;
 272     if (!currentnode) break;
 273     //    printf("%s snext=%f  density=%f bparam[0]=%f\n", gGeoManager->GetPath(),snext,density,bparam[0]);
 274     if (!gGeoManager->IsEntering()) {
 275       gGeoManager->SetStep(1.E-3);
 276       currentnode = gGeoManager->Step();
 277       if (!gGeoManager->IsEntering() || !currentnode) {
 278         //      printf("ERROR: cannot cross boundary\n");
 279         mparam[0] = bparam[0]/step;
 280         mparam[1] = bparam[1]/step;
 281         mparam[2] = bparam[2]/step;
 282         mparam[3] = bparam[3]/step;
 283         mparam[5] = bparam[5]/step;
 284         mparam[4] = step;
 285         mparam[0] = 0.;             // if crash of navigation take mean density 0
 286         mparam[1] = 1000000;        // and infinite rad length
 287         return bparam[0]/step;
 288       }
 289       step += 1.E-3;
 290       snext += 1.E-3;
 291       bparam[0] += lparam[0]*1.E-3;
 292       bparam[1]    += lparam[1]*1.E-3;
 293       bparam[2]    += lparam[2]*1.E-3;
 294       bparam[3]    += lparam[3]*1.E-3;
 295       bparam[5]    += lparam[5]*1.E-3;
 296     }
 297     length -= snext;
 298     material = currentnode->GetVolume()->GetMedium()->GetMaterial();
 299     lparam[0] = material->GetDensity();
 300     lparam[1]  = material->GetRadLen();
 301     lparam[2]  = material->GetA();
 302     lparam[3]  = material->GetZ();
 303     lparam[5]   = lparam[3]/lparam[2];
 304     if (material->IsMixture()) {
 305       lparam[1]*=lparam[0];
 306       TGeoMixture * mixture = (TGeoMixture*)material;
 307       lparam[5]=0;
 308       Double_t sum =0;
 309       for (Int_t iel=0;iel<mixture->GetNelements();iel++){
 310         sum+= mixture->GetWmixt()[iel];
 311         lparam[5]+= mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel];
 312       }
 313       lparam[5]/=sum;
 314     }
 315     gGeoManager->FindNextBoundary(length);
 316     snext = gGeoManager->GetStep();
 317   }
 318   mparam[0] = bparam[0]/step;
 319   mparam[1] = bparam[1]/step;
 320   mparam[2] = bparam[2]/step;
 321   mparam[3] = bparam[3]/step;
 322   mparam[5] = bparam[5]/step;
 323   return bparam[0]/step;
 324
 325 }
 326
 327 Double_t AliKalmanTrack::GetConvConst() {
 328   return 1000/0.299792458/AliTracker::GetBz();
 329 }
 330
 331 void AliKalmanTrack::SaveLocalConvConst() {
 332   //---------------------------------------------------------------------
 333   // Saves local conversion constant "curvature (1/cm) -> pt (GeV/c)"
 334   //---------------------------------------------------------------------
 335   if (AliTracker::UniformField()) {
 336      fLocalConvConst=1000/0.299792458/AliTracker::GetBz();
 337   } else {
 338      Float_t r[3]; GetXYZ(r);
 339      fLocalConvConst=1000/0.299792458/AliTracker::GetBz(r);
 340   }
 341 }
 342