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