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