]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/ESD/AliTrackerBase.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / STEER / ESD / AliTrackerBase.cxx
CommitLineData
4811a3f4 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
075f4221 20// that is the base for the AliTracker class
21// Origin: Marian.Ivanov@cern.ch
4811a3f4 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"
075f4221 30#include "AliTrackPointArray.h"
363db7c3 31#include "TVectorD.h"
4811a3f4 32
33extern TGeoManager *gGeoManager;
34
35ClassImp(AliTrackerBase)
36
37AliTrackerBase::AliTrackerBase():
38 TObject(),
39 fX(0),
40 fY(0),
41 fZ(0),
42 fSigmaX(0.005),
43 fSigmaY(0.005),
44 fSigmaZ(0.010)
45{
46 //--------------------------------------------------------------------
47 // The default constructor.
48 //--------------------------------------------------------------------
49 if (!TGeoGlobalMagField::Instance()->GetField())
50 AliWarning("Field map is not set.");
51}
52
53//__________________________________________________________________________
54AliTrackerBase::AliTrackerBase(const AliTrackerBase &atr):
55 TObject(atr),
56 fX(atr.fX),
57 fY(atr.fY),
58 fZ(atr.fZ),
59 fSigmaX(atr.fSigmaX),
60 fSigmaY(atr.fSigmaY),
61 fSigmaZ(atr.fSigmaZ)
62{
63 //--------------------------------------------------------------------
64 // The default constructor.
65 //--------------------------------------------------------------------
66 if (!TGeoGlobalMagField::Instance()->GetField())
67 AliWarning("Field map is not set.");
68}
69
70//__________________________________________________________________________
71Double_t AliTrackerBase::GetBz()
72{
73 AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
93ba142d 74 if (!fld) {
75 AliFatalClass("Field is not loaded");
76 //if (!fld)
77 return 0.5*kAlmost0Field;
78 }
4811a3f4 79 Double_t bz = fld->SolenoidField();
80 return TMath::Sign(0.5*kAlmost0Field,bz) + bz;
81}
82
83//__________________________________________________________________________
84Double_t AliTrackerBase::GetBz(const Double_t *r) {
85 //------------------------------------------------------------------
86 // Returns Bz (kG) at the point "r" .
87 //------------------------------------------------------------------
88 AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
93ba142d 89 if (!fld) {
90 AliFatalClass("Field is not loaded");
91 // if (!fld)
92 return 0.5*kAlmost0Field;
93 }
4811a3f4 94 Double_t bz = fld->GetBz(r);
95 return TMath::Sign(0.5*kAlmost0Field,bz) + bz;
96}
97
98//__________________________________________________________________________
99void AliTrackerBase::GetBxByBz(const Double_t r[3], Double_t b[3]) {
100 //------------------------------------------------------------------
101 // Returns Bx, By and Bz (kG) at the point "r" .
102 //------------------------------------------------------------------
103 AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
104 if (!fld) {
93ba142d 105 AliFatalClass("Field is not loaded");
106 // b[0] = b[1] = 0.;
107 // b[2] = 0.5*kAlmost0Field;
108 return;
4811a3f4 109 }
110
111 if (fld->IsUniform()) {
112 b[0] = b[1] = 0.;
113 b[2] = fld->SolenoidField();
114 } else {
115 fld->Field(r,b);
116 }
117 b[2] = (TMath::Sign(0.5*kAlmost0Field,b[2]) + b[2]);
118 return;
119}
120
121Double_t AliTrackerBase::MeanMaterialBudget(const Double_t *start, const Double_t *end, Double_t *mparam)
122{
123 //
124 // Calculate mean material budget and material properties between
125 // the points "start" and "end".
126 //
127 // "mparam" - parameters used for the energy and multiple scattering
128 // corrections:
129 //
130 // mparam[0] - mean density: sum(x_i*rho_i)/sum(x_i) [g/cm3]
131 // mparam[1] - equivalent rad length fraction: sum(x_i/X0_i) [adimensional]
132 // mparam[2] - mean A: sum(x_i*A_i)/sum(x_i) [adimensional]
133 // mparam[3] - mean Z: sum(x_i*Z_i)/sum(x_i) [adimensional]
134 // mparam[4] - length: sum(x_i) [cm]
135 // mparam[5] - Z/A mean: sum(x_i*Z_i/A_i)/sum(x_i) [adimensional]
136 // mparam[6] - number of boundary crosses
137 //
138 // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
139 //
140 // Corrections and improvements by
141 // Andrea Dainese, Andrea.Dainese@lnl.infn.it,
142 // Andrei Gheata, Andrei.Gheata@cern.ch
143 //
144
145 mparam[0]=0; mparam[1]=1; mparam[2] =0; mparam[3] =0;
146 mparam[4]=0; mparam[5]=0; mparam[6]=0;
147 //
148 Double_t bparam[6]; // total parameters
149 Double_t lparam[6]; // local parameters
150
151 for (Int_t i=0;i<6;i++) bparam[i]=0;
152
153 if (!gGeoManager) {
2ad9a021 154 AliFatalClass("No TGeo\n");
4811a3f4 155 return 0.;
156 }
157 //
158 Double_t length;
159 Double_t dir[3];
160 length = TMath::Sqrt((end[0]-start[0])*(end[0]-start[0])+
161 (end[1]-start[1])*(end[1]-start[1])+
162 (end[2]-start[2])*(end[2]-start[2]));
163 mparam[4]=length;
164 if (length<TGeoShape::Tolerance()) return 0.0;
165 Double_t invlen = 1./length;
166 dir[0] = (end[0]-start[0])*invlen;
167 dir[1] = (end[1]-start[1])*invlen;
168 dir[2] = (end[2]-start[2])*invlen;
169
170 // Initialize start point and direction
171 TGeoNode *currentnode = 0;
172 TGeoNode *startnode = gGeoManager->InitTrack(start, dir);
173 if (!startnode) {
174 AliErrorClass(Form("start point out of geometry: x %f, y %f, z %f",
175 start[0],start[1],start[2]));
176 return 0.0;
177 }
178 TGeoMaterial *material = startnode->GetVolume()->GetMedium()->GetMaterial();
179 lparam[0] = material->GetDensity();
180 lparam[1] = material->GetRadLen();
181 lparam[2] = material->GetA();
182 lparam[3] = material->GetZ();
183 lparam[4] = length;
184 lparam[5] = lparam[3]/lparam[2];
185 if (material->IsMixture()) {
186 TGeoMixture * mixture = (TGeoMixture*)material;
187 lparam[5] =0;
188 Double_t sum =0;
189 for (Int_t iel=0;iel<mixture->GetNelements();iel++){
190 sum += mixture->GetWmixt()[iel];
191 lparam[5]+= mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel];
192 }
193 lparam[5]/=sum;
194 }
195
196 // Locate next boundary within length without computing safety.
197 // Propagate either with length (if no boundary found) or just cross boundary
198 gGeoManager->FindNextBoundaryAndStep(length, kFALSE);
199 Double_t step = 0.0; // Step made
200 Double_t snext = gGeoManager->GetStep();
201 // If no boundary within proposed length, return current density
202 if (!gGeoManager->IsOnBoundary()) {
203 mparam[0] = lparam[0];
204 mparam[1] = lparam[4]/lparam[1];
205 mparam[2] = lparam[2];
206 mparam[3] = lparam[3];
207 mparam[4] = lparam[4];
208 return lparam[0];
209 }
210 // Try to cross the boundary and see what is next
211 Int_t nzero = 0;
212 while (length>TGeoShape::Tolerance()) {
213 currentnode = gGeoManager->GetCurrentNode();
214 if (snext<2.*TGeoShape::Tolerance()) nzero++;
215 else nzero = 0;
216 if (nzero>3) {
217 // This means navigation has problems on one boundary
218 // Try to cross by making a small step
219 AliErrorClass("Cannot cross boundary\n");
220 mparam[0] = bparam[0]/step;
221 mparam[1] = bparam[1];
222 mparam[2] = bparam[2]/step;
223 mparam[3] = bparam[3]/step;
224 mparam[5] = bparam[5]/step;
225 mparam[4] = step;
226 mparam[0] = 0.; // if crash of navigation take mean density 0
227 mparam[1] = 1000000; // and infinite rad length
228 return bparam[0]/step;
229 }
230 mparam[6]+=1.;
231 step += snext;
232 bparam[1] += snext/lparam[1];
233 bparam[2] += snext*lparam[2];
234 bparam[3] += snext*lparam[3];
235 bparam[5] += snext*lparam[5];
236 bparam[0] += snext*lparam[0];
237
238 if (snext>=length) break;
239 if (!currentnode) break;
240 length -= snext;
241 material = currentnode->GetVolume()->GetMedium()->GetMaterial();
242 lparam[0] = material->GetDensity();
243 lparam[1] = material->GetRadLen();
244 lparam[2] = material->GetA();
245 lparam[3] = material->GetZ();
246 lparam[5] = lparam[3]/lparam[2];
247 if (material->IsMixture()) {
248 TGeoMixture * mixture = (TGeoMixture*)material;
249 lparam[5]=0;
250 Double_t sum =0;
251 for (Int_t iel=0;iel<mixture->GetNelements();iel++){
252 sum+= mixture->GetWmixt()[iel];
253 lparam[5]+= mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel];
254 }
255 lparam[5]/=sum;
256 }
257 gGeoManager->FindNextBoundaryAndStep(length, kFALSE);
258 snext = gGeoManager->GetStep();
259 }
260 mparam[0] = bparam[0]/step;
261 mparam[1] = bparam[1];
262 mparam[2] = bparam[2]/step;
263 mparam[3] = bparam[3]/step;
264 mparam[5] = bparam[5]/step;
265 return bparam[0]/step;
266}
267
268
269Bool_t
270AliTrackerBase::PropagateTrackTo(AliExternalTrackParam *track, Double_t xToGo,
cd5c09f1 271 Double_t mass, Double_t maxStep, Bool_t rotateTo, Double_t maxSnp, Int_t sign, Bool_t addTimeStep, Bool_t correctMaterialBudget){
4811a3f4 272 //----------------------------------------------------------------
273 //
274 // Propagates the track to the plane X=xk (cm) using the magnetic field map
275 // and correcting for the crossed material.
276 //
1d26da6d 277 // mass - mass used in propagation - used for energy loss correction (if <0 then q = 2)
4811a3f4 278 // maxStep - maximal step for propagation
279 //
280 // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
281 //
282 //----------------------------------------------------------------
283 const Double_t kEpsilon = 0.00001;
284 Double_t xpos = track->GetX();
13545423 285 Int_t dir = (xpos<xToGo) ? 1:-1;
4811a3f4 286 //
287 while ( (xToGo-xpos)*dir > kEpsilon){
288 Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
289 Double_t x = xpos+step;
290 Double_t xyz0[3],xyz1[3],param[7];
291 track->GetXYZ(xyz0); //starting global position
292
293 Double_t bz=GetBz(xyz0); // getting the local Bz
4811a3f4 294 if (!track->GetXYZAt(x,bz,xyz1)) return kFALSE; // no prolongation
295 xyz1[2]+=kEpsilon; // waiting for bug correction in geo
cd5c09f1 296
7391b629 297 if (maxSnp>0 && TMath::Abs(track->GetSnpAt(x,bz)) >= maxSnp) return kFALSE;
4811a3f4 298 if (!track->PropagateTo(x,bz)) return kFALSE;
299
cd5c09f1 300 if (correctMaterialBudget){
301 MeanMaterialBudget(xyz0,xyz1,param);
302 Double_t xrho=param[0]*param[4], xx0=param[1];
303 if (sign) {if (sign<0) xrho = -xrho;} // sign is imposed
304 else { // determine automatically the sign from direction
305 if (dir>0) xrho = -xrho; // outward should be negative
306 }
307 //
308 if (!track->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE;
13545423 309 }
cd5c09f1 310
4811a3f4 311 if (rotateTo){
dc981549 312 track->GetXYZ(xyz1); // global position
313 Double_t alphan = TMath::ATan2(xyz1[1], xyz1[0]);
7391b629 314 if (maxSnp>0) {
315 if (TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE;
cd5c09f1 316
7391b629 317 //
318 Double_t ca=TMath::Cos(alphan-track->GetAlpha()), sa=TMath::Sin(alphan-track->GetAlpha());
319 Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf));
320 Double_t sinNew = sf*ca - cf*sa;
cd5c09f1 321 if (TMath::Abs(sinNew) >= maxSnp) return kFALSE;
322
7391b629 323 }
dc981549 324 if (!track->AliExternalTrackParam::Rotate(alphan)) return kFALSE;
cd5c09f1 325
4811a3f4 326 }
327 xpos = track->GetX();
13545423 328 if (addTimeStep && track->IsStartedTimeIntegral()) {
dc981549 329 if (!rotateTo) track->GetXYZ(xyz1); // if rotateTo==kTRUE, then xyz1 is already extracted
13545423 330 Double_t dX=xyz0[0]-xyz1[0],dY=xyz0[1]-xyz1[1],dZ=xyz0[2]-xyz1[2];
331 Double_t d=TMath::Sqrt(dX*dX + dY*dY + dZ*dZ);
332 if (sign) {if (sign>0) d = -d;} // step sign is imposed, positive means inward direction
333 else { // determine automatically the sign from direction
334 if (dir<0) d = -d;
335 }
336 track->AddTimeStep(d);
337 }
4811a3f4 338 }
339 return kTRUE;
340}
341
69749404 342Int_t AliTrackerBase::PropagateTrackTo2(AliExternalTrackParam *track, Double_t xToGo,
343 Double_t mass, Double_t maxStep, Bool_t rotateTo, Double_t maxSnp, Int_t sign, Bool_t addTimeStep, Bool_t correctMaterialBudget){
344 //----------------------------------------------------------------
345 //
346 // Propagates the track to the plane X=xk (cm) using the magnetic field map
347 // and correcting for the crossed material.
348 //
349 // mass - mass used in propagation - used for energy loss correction
350 // maxStep - maximal step for propagation
351 //
352 // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
353 //
354 //----------------------------------------------------------------
355 const Double_t kEpsilon = 0.00001;
356 Double_t xpos = track->GetX();
357 Int_t dir = (xpos<xToGo) ? 1:-1;
358 //
359 while ( (xToGo-xpos)*dir > kEpsilon){
360 Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
361 Double_t x = xpos+step;
362 Double_t xyz0[3],xyz1[3],param[7];
363 track->GetXYZ(xyz0); //starting global position
364
365 Double_t bz=GetBz(xyz0); // getting the local Bz
366 if (!track->GetXYZAt(x,bz,xyz1)) return -1; // no prolongation
367 xyz1[2]+=kEpsilon; // waiting for bug correction in geo
368
369 if (maxSnp>0 && TMath::Abs(track->GetSnpAt(x,bz)) >= maxSnp) return -2;
370 if (!track->PropagateTo(x,bz)) return -3;
371
372 if (correctMaterialBudget){
373 MeanMaterialBudget(xyz0,xyz1,param);
374 Double_t xrho=param[0]*param[4], xx0=param[1];
375 if (sign) {if (sign<0) xrho = -xrho;} // sign is imposed
376 else { // determine automatically the sign from direction
377 if (dir>0) xrho = -xrho; // outward should be negative
378 }
379 //
380 if (!track->CorrectForMeanMaterial(xx0,xrho,mass)) return -4;
381 }
382
383 if (rotateTo){
384 track->GetXYZ(xyz1); // global position
385 Double_t alphan = TMath::ATan2(xyz1[1], xyz1[0]);
386 if (maxSnp>0) {
387 if (TMath::Abs(track->GetSnp()) >= maxSnp) return -5;
388
389 //
390 Double_t ca=TMath::Cos(alphan-track->GetAlpha()), sa=TMath::Sin(alphan-track->GetAlpha());
391 Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf));
392 Double_t sinNew = sf*ca - cf*sa;
393 if (TMath::Abs(sinNew) >= maxSnp) return -6;
394
395 }
396 if (!track->AliExternalTrackParam::Rotate(alphan)) return -7;
397
398 }
399 xpos = track->GetX();
400 if (addTimeStep && track->IsStartedTimeIntegral()) {
401 if (!rotateTo) track->GetXYZ(xyz1); // if rotateTo==kTRUE, then xyz1 is already extracted
402 Double_t dX=xyz0[0]-xyz1[0],dY=xyz0[1]-xyz1[1],dZ=xyz0[2]-xyz1[2];
403 Double_t d=TMath::Sqrt(dX*dX + dY*dY + dZ*dZ);
404 if (sign) {if (sign>0) d = -d;} // step sign is imposed, positive means inward direction
405 else { // determine automatically the sign from direction
406 if (dir<0) d = -d;
407 }
408 track->AddTimeStep(d);
409 }
410 }
411 return 1;
412}
413
4811a3f4 414Bool_t
415AliTrackerBase::PropagateTrackToBxByBz(AliExternalTrackParam *track,
13545423 416 Double_t xToGo,Double_t mass, Double_t maxStep, Bool_t rotateTo, Double_t maxSnp,Int_t sign, Bool_t addTimeStep){
4811a3f4 417 //----------------------------------------------------------------
418 //
419 // Propagates the track to the plane X=xk (cm)
420 // taking into account all the three components of the magnetic field
421 // and correcting for the crossed material.
422 //
1d26da6d 423 // mass - mass used in propagation - used for energy loss correction (if <0 then q=2)
4811a3f4 424 // maxStep - maximal step for propagation
425 //
426 // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
427 //
428 //----------------------------------------------------------------
429 const Double_t kEpsilon = 0.00001;
430 Double_t xpos = track->GetX();
13545423 431 Int_t dir = (xpos<xToGo) ? 1:-1;
4811a3f4 432 //
433 while ( (xToGo-xpos)*dir > kEpsilon){
434 Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
435 Double_t x = xpos+step;
436 Double_t xyz0[3],xyz1[3],param[7];
437 track->GetXYZ(xyz0); //starting global position
438
439 Double_t b[3]; GetBxByBz(xyz0,b); // getting the local Bx, By and Bz
440
441 if (!track->GetXYZAt(x,b[2],xyz1)) return kFALSE; // no prolongation
442 xyz1[2]+=kEpsilon; // waiting for bug correction in geo
443
7391b629 444 if (maxSnp>0 && TMath::Abs(track->GetSnpAt(x,b[2])) >= maxSnp) return kFALSE;
4811a3f4 445 if (!track->PropagateToBxByBz(x,b)) return kFALSE;
446
13545423 447 MeanMaterialBudget(xyz0,xyz1,param);
448 Double_t xrho=param[0]*param[4], xx0=param[1];
449 if (sign) {if (sign<0) xrho = -xrho;} // sign is imposed
450 else { // determine automatically the sign from direction
451 if (dir>0) xrho = -xrho; // outward should be negative
452 }
453 //
4811a3f4 454 if (!track->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE;
455 if (rotateTo){
dc981549 456 track->GetXYZ(xyz1); // global position
457 Double_t alphan = TMath::ATan2(xyz1[1], xyz1[0]);
7391b629 458 if (maxSnp>0) {
459 if (TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE;
460 Double_t ca=TMath::Cos(alphan-track->GetAlpha()), sa=TMath::Sin(alphan-track->GetAlpha());
461 Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf));
462 Double_t sinNew = sf*ca - cf*sa;
463 if (TMath::Abs(sinNew) >= maxSnp) return kFALSE;
464 }
dc981549 465 if (!track->AliExternalTrackParam::Rotate(alphan)) return kFALSE;
4811a3f4 466 }
13545423 467 xpos = track->GetX();
468 if (addTimeStep && track->IsStartedTimeIntegral()) {
dc981549 469 if (!rotateTo) track->GetXYZ(xyz1); // if rotateTo==kTRUE, then xyz1 is already extracted
13545423 470 Double_t dX=xyz0[0]-xyz1[0],dY=xyz0[1]-xyz1[1],dZ=xyz0[2]-xyz1[2];
471 Double_t d=TMath::Sqrt(dX*dX + dY*dY + dZ*dZ);
472 if (sign) {if (sign>0) d = -d;} // step sign is imposed, positive means inward direction
473 else { // determine automatically the sign from direction
474 if (dir<0) d = -d;
475 }
476 track->AddTimeStep(d);
477 }
4811a3f4 478 }
479 return kTRUE;
480}
481
482Double_t AliTrackerBase::GetTrackPredictedChi2(AliExternalTrackParam *track,
483 Double_t mass, Double_t step,
484 const AliExternalTrackParam *backup) {
485 //
486 // This function brings the "track" with particle "mass" [GeV]
487 // to the same local coord. system and the same reference plane as
488 // of the "backup", doing it in "steps" [cm].
489 // Then, it calculates the 5D predicted Chi2 for these two tracks
490 //
491 Double_t chi2=kVeryBig;
492 Double_t alpha=backup->GetAlpha();
493 if (!track->Rotate(alpha)) return chi2;
494
495 Double_t xb=backup->GetX();
496 Double_t sign=(xb < track->GetX()) ? 1. : -1.;
497 if (!PropagateTrackTo(track,xb,mass,step,kFALSE,kAlmost1,sign)) return chi2;
498
499 chi2=track->GetPredictedChi2(backup);
500
501 return chi2;
502}
075f4221 503
504
505
506
507Double_t AliTrackerBase::MakeC(Double_t x1,Double_t y1,
508 Double_t x2,Double_t y2,
509 Double_t x3,Double_t y3)
510{
511 //-----------------------------------------------------------------
512 // Initial approximation of the track curvature
513 //-----------------------------------------------------------------
514 x3 -=x1;
515 x2 -=x1;
516 y3 -=y1;
517 y2 -=y1;
518 //
519 Double_t det = x3*y2-x2*y3;
520 if (TMath::Abs(det)<1e-10) {
521 return 0;
522 }
523 //
524 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
525 Double_t x0 = x3*0.5-y3*u;
526 Double_t y0 = y3*0.5+x3*u;
527 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
528 if (det>0) c2*=-1;
529 return c2;
530}
531
532Double_t AliTrackerBase::MakeSnp(Double_t x1,Double_t y1,
533 Double_t x2,Double_t y2,
534 Double_t x3,Double_t y3)
535{
536 //-----------------------------------------------------------------
537 // Initial approximation of the track snp
538 //-----------------------------------------------------------------
539 x3 -=x1;
540 x2 -=x1;
541 y3 -=y1;
542 y2 -=y1;
543 //
544 Double_t det = x3*y2-x2*y3;
545 if (TMath::Abs(det)<1e-10) {
546 return 0;
547 }
548 //
549 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
550 Double_t x0 = x3*0.5-y3*u;
551 Double_t y0 = y3*0.5+x3*u;
552 Double_t c2 = 1./TMath::Sqrt(x0*x0+y0*y0);
553 x0*=c2;
554 x0=TMath::Abs(x0);
555 if (y2*x2<0.) x0*=-1;
556 return x0;
557}
558
559Double_t AliTrackerBase::MakeTgl(Double_t x1,Double_t y1,
560 Double_t x2,Double_t y2,
561 Double_t z1,Double_t z2, Double_t c)
562{
563 //-----------------------------------------------------------------
564 // Initial approximation of the tangent of the track dip angle
565 //-----------------------------------------------------------------
566 //
e94b22cc 567 const Double_t kEpsilon =0.00001;
075f4221 568 x2-=x1;
569 y2-=y1;
570 z2-=z1;
571 Double_t d = TMath::Sqrt(x2*x2+y2*y2); // distance straight line
572 if (TMath::Abs(d*c*0.5)>1) return 0;
573 Double_t angle2 = TMath::ASin(d*c*0.5);
e94b22cc 574 if (TMath::Abs(angle2)>kEpsilon) {
575 angle2 = z2*TMath::Abs(c/(angle2*2.));
576 }else{
577 angle2=z2/d;
578 }
075f4221 579 return angle2;
580}
581
582
583Double_t AliTrackerBase::MakeTgl(Double_t x1,Double_t y1,
584 Double_t x2,Double_t y2,
585 Double_t z1,Double_t z2)
586{
587 //-----------------------------------------------------------------
588 // Initial approximation of the tangent of the track dip angle
589 //-----------------------------------------------------------------
590 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
591}
592
593
594AliExternalTrackParam * AliTrackerBase::MakeSeed( AliTrackPoint &point0, AliTrackPoint &point1, AliTrackPoint &point2){
595 //
596 // Make Seed - AliExternalTrackParam from input 3 points
597 // returning seed in local frame of point0
598 //
599 Double_t xyz0[3]={0,0,0};
600 Double_t xyz1[3]={0,0,0};
601 Double_t xyz2[3]={0,0,0};
602 Double_t alpha=point0.GetAngle();
603 Double_t xyz[3]={point0.GetX(),point0.GetY(),point0.GetZ()};
604 Double_t bxyz[3]; GetBxByBz(xyz,bxyz);
605 Double_t bz = bxyz[2];
606 //
607 // get points in frame of point 0
608 //
609 AliTrackPoint p0r = point0.Rotate(alpha);
610 AliTrackPoint p1r = point1.Rotate(alpha);
611 AliTrackPoint p2r = point2.Rotate(alpha);
612 xyz0[0]=p0r.GetX();
613 xyz0[1]=p0r.GetY();
614 xyz0[2]=p0r.GetZ();
615 xyz1[0]=p1r.GetX();
616 xyz1[1]=p1r.GetY();
617 xyz1[2]=p1r.GetZ();
618 xyz2[0]=p2r.GetX();
619 xyz2[1]=p2r.GetY();
620 xyz2[2]=p2r.GetZ();
621 //
622 // make covariance estimate
623 //
624 Double_t covar[15];
625 Double_t param[5]={0,0,0,0,0};
626 for (Int_t m=0; m<15; m++) covar[m]=0;
627 //
628 // calculate intitial param
629 param[0]=xyz0[1];
630 param[1]=xyz0[2];
631 param[2]=MakeSnp(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]);
632 param[4]=MakeC(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]);
633 param[3]=MakeTgl(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz0[2],xyz1[2],param[4]);
634
635 //covariance matrix - only diagonal elements
636 //Double_t dist=p0r.GetX()-p2r.GetX();
637 Double_t deltaP=0;
638 covar[0]= p0r.GetCov()[3];
639 covar[2]= p0r.GetCov()[5];
640 //sigma snp
641 deltaP= (MakeSnp(xyz0[0],xyz0[1]+TMath::Sqrt(p0r.GetCov()[3]),xyz1[0],xyz1[1],xyz2[0],xyz2[1])-param[2]);
642 covar[5]+= deltaP*deltaP;
643 deltaP= (MakeSnp(xyz0[0],xyz0[1],xyz1[0],xyz1[1]+TMath::Sqrt(p1r.GetCov()[3]),xyz2[0],xyz2[1])-param[2]);
644 covar[5]+= deltaP*deltaP;
645 deltaP= (MakeSnp(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]+TMath::Sqrt(p1r.GetCov()[3]))-param[2]);
646 covar[5]+= deltaP*deltaP;
647 //sigma tgl
648 //
649 deltaP=MakeTgl(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz0[2]+TMath::Sqrt(p1r.GetCov()[5]),xyz1[2],param[4])-param[3];
650 covar[9]+= deltaP*deltaP;
651 deltaP=MakeTgl(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz0[2],xyz1[2]+TMath::Sqrt(p1r.GetCov()[5]),param[4])-param[3];
652 covar[9]+= deltaP*deltaP;
653 //
654
655 deltaP=MakeC(xyz0[0],xyz0[1]+TMath::Sqrt(p0r.GetCov()[3]),xyz1[0],xyz1[1],xyz2[0],xyz2[1])-param[4];
656 covar[14]+= deltaP*deltaP;
657 deltaP=MakeC(xyz0[0],xyz0[1],xyz1[0],xyz1[1]+TMath::Sqrt(p1r.GetCov()[3]),xyz2[0],xyz2[1])-param[4];
658 covar[14]+= deltaP*deltaP;
659 deltaP=MakeC(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]+TMath::Sqrt(p2r.GetCov()[3]))-param[4];
660 covar[14]+= deltaP*deltaP;
661
662 covar[14]/=(bz*kB2C)*(bz*kB2C);
663 param[4]/=(bz*kB2C); // transform to 1/pt
664 AliExternalTrackParam * trackParam = new AliExternalTrackParam(xyz0[0],alpha,param, covar);
665 if (0) {
666 // consistency check -to put warnings here
667 // small disagrement once Track extrapolation used
668 // nice agreement in seeds with MC track parameters - problem in extrapoloation - to be fixed
669 // to check later
670 Double_t y1,y2,z1,z2;
671 trackParam->GetYAt(xyz1[0],bz,y1);
672 trackParam->GetZAt(xyz1[0],bz,z1);
673 trackParam->GetYAt(xyz2[0],bz,y2);
674 trackParam->GetZAt(xyz2[0],bz,z2);
675 if (TMath::Abs(y1-xyz1[1])> TMath::Sqrt(p1r.GetCov()[3]*5)){
676 AliWarningClass("Seeding problem y1\n");
677 }
678 if (TMath::Abs(y2-xyz2[1])> TMath::Sqrt(p2r.GetCov()[3]*5)){
679 AliWarningClass("Seeding problem y2\n");
680 }
681 if (TMath::Abs(z1-xyz1[2])> TMath::Sqrt(p1r.GetCov()[5]*5)){
682 AliWarningClass("Seeding problem z1\n");
683 }
684 }
685 return trackParam;
686}
687
688Double_t AliTrackerBase::FitTrack(AliExternalTrackParam * trackParam, AliTrackPointArray *pointArray, Double_t mass, Double_t maxStep){
689 //
690 // refit the track - trackParam using the points in point array
691 //
692 const Double_t kMaxSnp=0.99;
693 if (!trackParam) return 0;
694 Int_t npoints=pointArray->GetNPoints();
695 AliTrackPoint point,point2;
696 Double_t pointPos[2]={0,0};
697 Double_t pointCov[3]={0,0,0};
698 // choose coordinate frame
699 // in standard way the coordinate frame should be changed point by point
700 // Some problems with rotation observed
701 // rotate method of AliExternalTrackParam should be revisited
702 pointArray->GetPoint(point,0);
703 pointArray->GetPoint(point2,npoints-1);
704 Double_t alpha=TMath::ATan2(point.GetY()-point2.GetY(), point.GetX()-point2.GetX());
705
706 for (Int_t ipoint=npoints-1; ipoint>0; ipoint-=1){
707 pointArray->GetPoint(point,ipoint);
708 AliTrackPoint pr = point.Rotate(alpha);
709 trackParam->Rotate(alpha);
710 Bool_t status = PropagateTrackTo(trackParam,pr.GetX(),mass,maxStep,kFALSE,kMaxSnp);
711 if(!status){
712 AliWarningClass("Problem to propagate\n");
713 break;
714 }
715 if (TMath::Abs(trackParam->GetSnp())>kMaxSnp){
716 AliWarningClass("sin(phi) > kMaxSnp \n");
717 break;
718 }
719 pointPos[0]=pr.GetY();//local y
720 pointPos[1]=pr.GetZ();//local z
721 pointCov[0]=pr.GetCov()[3];//simay^2
722 pointCov[1]=pr.GetCov()[4];//sigmayz
723 pointCov[2]=pr.GetCov()[5];//sigmaz^2
724 trackParam->Update(pointPos,pointCov);
725 }
726 return 0;
727}
363db7c3 728
729
730
731void AliTrackerBase::UpdateTrack(AliExternalTrackParam &track1, const AliExternalTrackParam &track2){
732 //
733 // Update track 1 with track 2
734 //
735 //
736 //
737 TMatrixD vecXk(5,1); // X vector
738 TMatrixD covXk(5,5); // X covariance
739 TMatrixD matHk(5,5); // vector to mesurement
740 TMatrixD measR(5,5); // measurement error
741 TMatrixD vecZk(5,1); // measurement
742 //
743 TMatrixD vecYk(5,1); // Innovation or measurement residual
744 TMatrixD matHkT(5,5);
745 TMatrixD matSk(5,5); // Innovation (or residual) covariance
746 TMatrixD matKk(5,5); // Optimal Kalman gain
747 TMatrixD mat1(5,5); // update covariance matrix
748 TMatrixD covXk2(5,5); //
749 TMatrixD covOut(5,5);
750 //
751 Double_t *param1=(Double_t*) track1.GetParameter();
752 Double_t *covar1=(Double_t*) track1.GetCovariance();
753 Double_t *param2=(Double_t*) track2.GetParameter();
754 Double_t *covar2=(Double_t*) track2.GetCovariance();
755 //
756 // copy data to the matrix
757 for (Int_t ipar=0; ipar<5; ipar++){
758 for (Int_t jpar=0; jpar<5; jpar++){
759 covXk(ipar,jpar) = covar1[track1.GetIndex(ipar, jpar)];
760 measR(ipar,jpar) = covar2[track2.GetIndex(ipar, jpar)];
761 matHk(ipar,jpar)=0;
762 mat1(ipar,jpar)=0;
763 }
764 vecXk(ipar,0) = param1[ipar];
765 vecZk(ipar,0) = param2[ipar];
766 matHk(ipar,ipar)=1;
767 mat1(ipar,ipar)=1;
768 }
769 //
770 //
771 //
772 //
773 //
774 vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
775 matHkT=matHk.T(); matHk.T();
776 matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance
777 matSk.Invert();
778 matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
779 vecXk += matKk*vecYk; // updated vector
780 covXk2 = (mat1-(matKk*matHk));
781 covOut = covXk2*covXk;
782 //
783 //
784 //
785 // copy from matrix to parameters
786 if (0) {
787 vecXk.Print();
788 vecZk.Print();
789 //
790 measR.Print();
791 covXk.Print();
792 covOut.Print();
793 //
794 track1.Print();
795 track2.Print();
796 }
797
798 for (Int_t ipar=0; ipar<5; ipar++){
799 param1[ipar]= vecXk(ipar,0) ;
800 for (Int_t jpar=0; jpar<5; jpar++){
801 covar1[track1.GetIndex(ipar, jpar)]=covOut(ipar,jpar);
802 }
803 }
804}
805
806