Possibility to calculate the DCA between two ESD track. The V0 and cascade vertexes...
[u/mrichter/AliRoot.git] / STEER / AliKalmanTrack.cxx
CommitLineData
87594435 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
acd84897 16/* $Id$ */
fb17acd4 17
87594435 18//-------------------------------------------------------------------------
19// Implementation of the AliKalmanTrack class
066782e8 20// that is the base for AliTPCtrack, AliITStrackV2 and AliTRDtrack
87594435 21// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
22//-------------------------------------------------------------------------
7d0f8548 23#include "AliTracker.h"
87594435 24#include "AliKalmanTrack.h"
4557b520 25#include "TGeoManager.h"
87594435 26
27ClassImp(AliKalmanTrack)
28
e2afb3b6 29//_______________________________________________________________________
30AliKalmanTrack::AliKalmanTrack():
31 fLab(-3141593),
68b8060b 32 fFakeRatio(0),
e2afb3b6 33 fChi2(0),
304864ab 34 fMass(AliPID::ParticleMass(AliPID::kPion)),
90e48c0c 35 fN(0),
c84a5e9e 36 fLocalConvConst(0),
90e48c0c 37 fStartTimeIntegral(kFALSE),
38 fIntegratedLength(0)
e2afb3b6 39{
116cbefd 40 //
41 // Default constructor
42 //
7d0f8548 43 if (AliTracker::GetFieldMap()==0) {
f37d970d 44 AliFatal("The magnetic field has not been set!");
7d0f8548 45 }
c84a5e9e 46
7d0f8548 47 for(Int_t i=0; i<AliPID::kSPECIES; i++) fIntegratedTime[i] = 0;
e2afb3b6 48}
49
50//_______________________________________________________________________
51AliKalmanTrack::AliKalmanTrack(const AliKalmanTrack &t):
52 TObject(t),
53 fLab(t.fLab),
babd135a 54 fFakeRatio(t.fFakeRatio),
e2afb3b6 55 fChi2(t.fChi2),
56 fMass(t.fMass),
90e48c0c 57 fN(t.fN),
c84a5e9e 58 fLocalConvConst(t.fLocalConvConst),
90e48c0c 59 fStartTimeIntegral(t.fStartTimeIntegral),
60 fIntegratedLength(t.fIntegratedLength)
e2afb3b6 61{
116cbefd 62 //
63 // Copy constructor
64 //
7d0f8548 65 if (AliTracker::GetFieldMap()==0) {
f37d970d 66 AliFatal("The magnetic field has not been set!");
8de97894 67 }
74f9526e 68
c84a5e9e 69 for (Int_t i=0; i<AliPID::kSPECIES; i++)
70 fIntegratedTime[i] = t.fIntegratedTime[i];
74f9526e 71}
c5507f6d 72
73//_______________________________________________________________________
74f9526e 74void AliKalmanTrack::StartTimeIntegral()
75{
49a7a79a 76 // Sylwester Radomski, GSI
77 // S.Radomski@gsi.de
74f9526e 78 //
79 // Start time integration
80 // To be called at Vertex by ITS tracker
81 //
82
83 //if (fStartTimeIntegral)
f37d970d 84 // AliWarning("Reseting Recorded Time.");
74f9526e 85
86 fStartTimeIntegral = kTRUE;
304864ab 87 for(Int_t i=0; i<AliPID::kSPECIES; i++) fIntegratedTime[i] = 0;
74f9526e 88 fIntegratedLength = 0;
89}
7d0f8548 90
74f9526e 91//_______________________________________________________________________
92void 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
5d8718b8 114 static const Double_t kcc = 2.99792458e-2;
74f9526e 115
116 if (!fStartTimeIntegral) return;
117
118 fIntegratedLength += length;
119
74f9526e 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
304864ab 131 for (Int_t i=0; i<AliPID::kSPECIES; i++) {
74f9526e 132
304864ab 133 Double_t mass = AliPID::ParticleMass(i);
74f9526e 134 Double_t correction = TMath::Sqrt( pt*pt * (1 + tgl*tgl) + mass * mass ) / p;
5d8718b8 135 Double_t time = length * correction / kcc;
74f9526e 136
74f9526e 137 fIntegratedTime[i] += time;
138 }
e2afb3b6 139}
140
74f9526e 141//_______________________________________________________________________
74f9526e 142Double_t AliKalmanTrack::GetIntegratedTime(Int_t pdg) const
143{
49a7a79a 144 // Sylwester Radomski, GSI
145 // S.Radomski@gsi.de
74f9526e 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) {
f37d970d 156 AliWarning("Time integration not started");
74f9526e 157 return 0.;
158 }
159
304864ab 160 for (Int_t i=0; i<AliPID::kSPECIES; i++)
161 if (AliPID::ParticleCode(i) == TMath::Abs(pdg)) return fIntegratedTime[i];
74f9526e 162
f37d970d 163 AliWarning(Form("Particle type [%d] not found", pdg));
74f9526e 164 return 0;
165}
ae982df3 166
167void AliKalmanTrack::GetIntegratedTimes(Double_t *times) const {
304864ab 168 for (Int_t i=0; i<AliPID::kSPECIES; i++) times[i]=fIntegratedTime[i];
ae982df3 169}
170
171void AliKalmanTrack::SetIntegratedTimes(const Double_t *times) {
304864ab 172 for (Int_t i=0; i<AliPID::kSPECIES; i++) fIntegratedTime[i]=times[i];
ae982df3 173}
174
c84a5e9e 175void AliKalmanTrack::External2Helix(Double_t helix[6]) const {
49a7a79a 176 //--------------------------------------------------------------------
177 // External track parameters -> helix parameters
178 //--------------------------------------------------------------------
179 Double_t alpha,x,cs,sn;
c84a5e9e 180 GetExternalParameters(x,helix); alpha=GetAlpha();
49a7a79a 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
c84a5e9e 188 helix[4]=helix[4]/GetLocalConvConst(); // C
49a7a79a 189}
190
4557b520 191Double_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();
5df14aca 239 lparam[4] = length;
4557b520 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;
58d96064 272 if (!currentnode) break;
4557b520 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}
7d0f8548 326
327Double_t AliKalmanTrack::GetConvConst() {
328 return 1000/0.299792458/AliTracker::GetBz();
329}
330
331void 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