]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/AliTrackerBase.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Bug Fix: Forgot to call the initialization of
[u/mrichter/AliRoot.git] / STEER / 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"
4811a3f4 31
32extern TGeoManager *gGeoManager;
33
34ClassImp(AliTrackerBase)
35
36AliTrackerBase::AliTrackerBase():
37 TObject(),
38 fX(0),
39 fY(0),
40 fZ(0),
41 fSigmaX(0.005),
42 fSigmaY(0.005),
43 fSigmaZ(0.010)
44{
45 //--------------------------------------------------------------------
46 // The default constructor.
47 //--------------------------------------------------------------------
48 if (!TGeoGlobalMagField::Instance()->GetField())
49 AliWarning("Field map is not set.");
50}
51
52//__________________________________________________________________________
53AliTrackerBase::AliTrackerBase(const AliTrackerBase &atr):
54 TObject(atr),
55 fX(atr.fX),
56 fY(atr.fY),
57 fZ(atr.fZ),
58 fSigmaX(atr.fSigmaX),
59 fSigmaY(atr.fSigmaY),
60 fSigmaZ(atr.fSigmaZ)
61{
62 //--------------------------------------------------------------------
63 // The default constructor.
64 //--------------------------------------------------------------------
65 if (!TGeoGlobalMagField::Instance()->GetField())
66 AliWarning("Field map is not set.");
67}
68
69//__________________________________________________________________________
70Double_t AliTrackerBase::GetBz()
71{
72 AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
73 if (!fld) return 0.5*kAlmost0Field;
74 Double_t bz = fld->SolenoidField();
75 return TMath::Sign(0.5*kAlmost0Field,bz) + bz;
76}
77
78//__________________________________________________________________________
79Double_t AliTrackerBase::GetBz(const Double_t *r) {
80 //------------------------------------------------------------------
81 // Returns Bz (kG) at the point "r" .
82 //------------------------------------------------------------------
83 AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
84 if (!fld) return 0.5*kAlmost0Field;
85 Double_t bz = fld->GetBz(r);
86 return TMath::Sign(0.5*kAlmost0Field,bz) + bz;
87}
88
89//__________________________________________________________________________
90void AliTrackerBase::GetBxByBz(const Double_t r[3], Double_t b[3]) {
91 //------------------------------------------------------------------
92 // Returns Bx, By and Bz (kG) at the point "r" .
93 //------------------------------------------------------------------
94 AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
95 if (!fld) {
96 b[0] = b[1] = 0.;
97 b[2] = 0.5*kAlmost0Field;
98 return;
99 }
100
101 if (fld->IsUniform()) {
102 b[0] = b[1] = 0.;
103 b[2] = fld->SolenoidField();
104 } else {
105 fld->Field(r,b);
106 }
107 b[2] = (TMath::Sign(0.5*kAlmost0Field,b[2]) + b[2]);
108 return;
109}
110
111Double_t AliTrackerBase::MeanMaterialBudget(const Double_t *start, const Double_t *end, Double_t *mparam)
112{
113 //
114 // Calculate mean material budget and material properties between
115 // the points "start" and "end".
116 //
117 // "mparam" - parameters used for the energy and multiple scattering
118 // corrections:
119 //
120 // mparam[0] - mean density: sum(x_i*rho_i)/sum(x_i) [g/cm3]
121 // mparam[1] - equivalent rad length fraction: sum(x_i/X0_i) [adimensional]
122 // mparam[2] - mean A: sum(x_i*A_i)/sum(x_i) [adimensional]
123 // mparam[3] - mean Z: sum(x_i*Z_i)/sum(x_i) [adimensional]
124 // mparam[4] - length: sum(x_i) [cm]
125 // mparam[5] - Z/A mean: sum(x_i*Z_i/A_i)/sum(x_i) [adimensional]
126 // mparam[6] - number of boundary crosses
127 //
128 // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
129 //
130 // Corrections and improvements by
131 // Andrea Dainese, Andrea.Dainese@lnl.infn.it,
132 // Andrei Gheata, Andrei.Gheata@cern.ch
133 //
134
135 mparam[0]=0; mparam[1]=1; mparam[2] =0; mparam[3] =0;
136 mparam[4]=0; mparam[5]=0; mparam[6]=0;
137 //
138 Double_t bparam[6]; // total parameters
139 Double_t lparam[6]; // local parameters
140
141 for (Int_t i=0;i<6;i++) bparam[i]=0;
142
143 if (!gGeoManager) {
144 AliErrorClass("No TGeo\n");
145 return 0.;
146 }
147 //
148 Double_t length;
149 Double_t dir[3];
150 length = TMath::Sqrt((end[0]-start[0])*(end[0]-start[0])+
151 (end[1]-start[1])*(end[1]-start[1])+
152 (end[2]-start[2])*(end[2]-start[2]));
153 mparam[4]=length;
154 if (length<TGeoShape::Tolerance()) return 0.0;
155 Double_t invlen = 1./length;
156 dir[0] = (end[0]-start[0])*invlen;
157 dir[1] = (end[1]-start[1])*invlen;
158 dir[2] = (end[2]-start[2])*invlen;
159
160 // Initialize start point and direction
161 TGeoNode *currentnode = 0;
162 TGeoNode *startnode = gGeoManager->InitTrack(start, dir);
163 if (!startnode) {
164 AliErrorClass(Form("start point out of geometry: x %f, y %f, z %f",
165 start[0],start[1],start[2]));
166 return 0.0;
167 }
168 TGeoMaterial *material = startnode->GetVolume()->GetMedium()->GetMaterial();
169 lparam[0] = material->GetDensity();
170 lparam[1] = material->GetRadLen();
171 lparam[2] = material->GetA();
172 lparam[3] = material->GetZ();
173 lparam[4] = length;
174 lparam[5] = lparam[3]/lparam[2];
175 if (material->IsMixture()) {
176 TGeoMixture * mixture = (TGeoMixture*)material;
177 lparam[5] =0;
178 Double_t sum =0;
179 for (Int_t iel=0;iel<mixture->GetNelements();iel++){
180 sum += mixture->GetWmixt()[iel];
181 lparam[5]+= mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel];
182 }
183 lparam[5]/=sum;
184 }
185
186 // Locate next boundary within length without computing safety.
187 // Propagate either with length (if no boundary found) or just cross boundary
188 gGeoManager->FindNextBoundaryAndStep(length, kFALSE);
189 Double_t step = 0.0; // Step made
190 Double_t snext = gGeoManager->GetStep();
191 // If no boundary within proposed length, return current density
192 if (!gGeoManager->IsOnBoundary()) {
193 mparam[0] = lparam[0];
194 mparam[1] = lparam[4]/lparam[1];
195 mparam[2] = lparam[2];
196 mparam[3] = lparam[3];
197 mparam[4] = lparam[4];
198 return lparam[0];
199 }
200 // Try to cross the boundary and see what is next
201 Int_t nzero = 0;
202 while (length>TGeoShape::Tolerance()) {
203 currentnode = gGeoManager->GetCurrentNode();
204 if (snext<2.*TGeoShape::Tolerance()) nzero++;
205 else nzero = 0;
206 if (nzero>3) {
207 // This means navigation has problems on one boundary
208 // Try to cross by making a small step
209 AliErrorClass("Cannot cross boundary\n");
210 mparam[0] = bparam[0]/step;
211 mparam[1] = bparam[1];
212 mparam[2] = bparam[2]/step;
213 mparam[3] = bparam[3]/step;
214 mparam[5] = bparam[5]/step;
215 mparam[4] = step;
216 mparam[0] = 0.; // if crash of navigation take mean density 0
217 mparam[1] = 1000000; // and infinite rad length
218 return bparam[0]/step;
219 }
220 mparam[6]+=1.;
221 step += snext;
222 bparam[1] += snext/lparam[1];
223 bparam[2] += snext*lparam[2];
224 bparam[3] += snext*lparam[3];
225 bparam[5] += snext*lparam[5];
226 bparam[0] += snext*lparam[0];
227
228 if (snext>=length) break;
229 if (!currentnode) break;
230 length -= snext;
231 material = currentnode->GetVolume()->GetMedium()->GetMaterial();
232 lparam[0] = material->GetDensity();
233 lparam[1] = material->GetRadLen();
234 lparam[2] = material->GetA();
235 lparam[3] = material->GetZ();
236 lparam[5] = lparam[3]/lparam[2];
237 if (material->IsMixture()) {
238 TGeoMixture * mixture = (TGeoMixture*)material;
239 lparam[5]=0;
240 Double_t sum =0;
241 for (Int_t iel=0;iel<mixture->GetNelements();iel++){
242 sum+= mixture->GetWmixt()[iel];
243 lparam[5]+= mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel];
244 }
245 lparam[5]/=sum;
246 }
247 gGeoManager->FindNextBoundaryAndStep(length, kFALSE);
248 snext = gGeoManager->GetStep();
249 }
250 mparam[0] = bparam[0]/step;
251 mparam[1] = bparam[1];
252 mparam[2] = bparam[2]/step;
253 mparam[3] = bparam[3]/step;
254 mparam[5] = bparam[5]/step;
255 return bparam[0]/step;
256}
257
258
259Bool_t
260AliTrackerBase::PropagateTrackTo(AliExternalTrackParam *track, Double_t xToGo,
261 Double_t mass, Double_t maxStep, Bool_t rotateTo, Double_t maxSnp, Double_t sign){
262 //----------------------------------------------------------------
263 //
264 // Propagates the track to the plane X=xk (cm) using the magnetic field map
265 // and correcting for the crossed material.
266 //
267 // mass - mass used in propagation - used for energy loss correction
268 // maxStep - maximal step for propagation
269 //
270 // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
271 //
272 //----------------------------------------------------------------
273 const Double_t kEpsilon = 0.00001;
274 Double_t xpos = track->GetX();
275 Double_t dir = (xpos<xToGo) ? 1.:-1.;
276 //
277 while ( (xToGo-xpos)*dir > kEpsilon){
278 Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
279 Double_t x = xpos+step;
280 Double_t xyz0[3],xyz1[3],param[7];
281 track->GetXYZ(xyz0); //starting global position
282
283 Double_t bz=GetBz(xyz0); // getting the local Bz
284
285 if (!track->GetXYZAt(x,bz,xyz1)) return kFALSE; // no prolongation
286 xyz1[2]+=kEpsilon; // waiting for bug correction in geo
287
288 if (TMath::Abs(track->GetSnpAt(x,bz)) >= maxSnp) return kFALSE;
289 if (!track->PropagateTo(x,bz)) return kFALSE;
290
291 MeanMaterialBudget(xyz0,xyz1,param);
292 Double_t xrho=param[0]*param[4]*sign, xx0=param[1];
293
294 if (!track->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE;
295 if (rotateTo){
296 if (TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE;
297 track->GetXYZ(xyz0); // global position
298 Double_t alphan = TMath::ATan2(xyz0[1], xyz0[0]);
299 //
300 Double_t ca=TMath::Cos(alphan-track->GetAlpha()),
301 sa=TMath::Sin(alphan-track->GetAlpha());
302 Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf));
303 Double_t sinNew = sf*ca - cf*sa;
304 if (TMath::Abs(sinNew) >= maxSnp) return kFALSE;
305 if (!track->Rotate(alphan)) return kFALSE;
306 }
307 xpos = track->GetX();
308 }
309 return kTRUE;
310}
311
312Bool_t
313AliTrackerBase::PropagateTrackToBxByBz(AliExternalTrackParam *track,
314Double_t xToGo,
315 Double_t mass, Double_t maxStep, Bool_t rotateTo, Double_t maxSnp,Double_t sign){
316 //----------------------------------------------------------------
317 //
318 // Propagates the track to the plane X=xk (cm)
319 // taking into account all the three components of the magnetic field
320 // and correcting for the crossed material.
321 //
322 // mass - mass used in propagation - used for energy loss correction
323 // maxStep - maximal step for propagation
324 //
325 // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
326 //
327 //----------------------------------------------------------------
328 const Double_t kEpsilon = 0.00001;
329 Double_t xpos = track->GetX();
330 Double_t dir = (xpos<xToGo) ? 1.:-1.;
331 //
332 while ( (xToGo-xpos)*dir > kEpsilon){
333 Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
334 Double_t x = xpos+step;
335 Double_t xyz0[3],xyz1[3],param[7];
336 track->GetXYZ(xyz0); //starting global position
337
338 Double_t b[3]; GetBxByBz(xyz0,b); // getting the local Bx, By and Bz
339
340 if (!track->GetXYZAt(x,b[2],xyz1)) return kFALSE; // no prolongation
341 xyz1[2]+=kEpsilon; // waiting for bug correction in geo
342
343 if (TMath::Abs(track->GetSnpAt(x,b[2])) >= maxSnp) return kFALSE;
344 if (!track->PropagateToBxByBz(x,b)) return kFALSE;
345
346 MeanMaterialBudget(xyz0,xyz1,param);
347 Double_t xrho=param[0]*param[4]*sign, xx0=param[1];
348
349 if (!track->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE;
350 if (rotateTo){
351 if (TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE;
352 track->GetXYZ(xyz0); // global position
353 Double_t alphan = TMath::ATan2(xyz0[1], xyz0[0]);
354 //
355 Double_t ca=TMath::Cos(alphan-track->GetAlpha()),
356 sa=TMath::Sin(alphan-track->GetAlpha());
357 Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf));
358 Double_t sinNew = sf*ca - cf*sa;
359 if (TMath::Abs(sinNew) >= maxSnp) return kFALSE;
360 if (!track->Rotate(alphan)) return kFALSE;
361 }
362 xpos = track->GetX();
363 }
364 return kTRUE;
365}
366
367Double_t AliTrackerBase::GetTrackPredictedChi2(AliExternalTrackParam *track,
368 Double_t mass, Double_t step,
369 const AliExternalTrackParam *backup) {
370 //
371 // This function brings the "track" with particle "mass" [GeV]
372 // to the same local coord. system and the same reference plane as
373 // of the "backup", doing it in "steps" [cm].
374 // Then, it calculates the 5D predicted Chi2 for these two tracks
375 //
376 Double_t chi2=kVeryBig;
377 Double_t alpha=backup->GetAlpha();
378 if (!track->Rotate(alpha)) return chi2;
379
380 Double_t xb=backup->GetX();
381 Double_t sign=(xb < track->GetX()) ? 1. : -1.;
382 if (!PropagateTrackTo(track,xb,mass,step,kFALSE,kAlmost1,sign)) return chi2;
383
384 chi2=track->GetPredictedChi2(backup);
385
386 return chi2;
387}
075f4221 388
389
390
391
392Double_t AliTrackerBase::MakeC(Double_t x1,Double_t y1,
393 Double_t x2,Double_t y2,
394 Double_t x3,Double_t y3)
395{
396 //-----------------------------------------------------------------
397 // Initial approximation of the track curvature
398 //-----------------------------------------------------------------
399 x3 -=x1;
400 x2 -=x1;
401 y3 -=y1;
402 y2 -=y1;
403 //
404 Double_t det = x3*y2-x2*y3;
405 if (TMath::Abs(det)<1e-10) {
406 return 0;
407 }
408 //
409 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
410 Double_t x0 = x3*0.5-y3*u;
411 Double_t y0 = y3*0.5+x3*u;
412 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
413 if (det>0) c2*=-1;
414 return c2;
415}
416
417Double_t AliTrackerBase::MakeSnp(Double_t x1,Double_t y1,
418 Double_t x2,Double_t y2,
419 Double_t x3,Double_t y3)
420{
421 //-----------------------------------------------------------------
422 // Initial approximation of the track snp
423 //-----------------------------------------------------------------
424 x3 -=x1;
425 x2 -=x1;
426 y3 -=y1;
427 y2 -=y1;
428 //
429 Double_t det = x3*y2-x2*y3;
430 if (TMath::Abs(det)<1e-10) {
431 return 0;
432 }
433 //
434 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
435 Double_t x0 = x3*0.5-y3*u;
436 Double_t y0 = y3*0.5+x3*u;
437 Double_t c2 = 1./TMath::Sqrt(x0*x0+y0*y0);
438 x0*=c2;
439 x0=TMath::Abs(x0);
440 if (y2*x2<0.) x0*=-1;
441 return x0;
442}
443
444Double_t AliTrackerBase::MakeTgl(Double_t x1,Double_t y1,
445 Double_t x2,Double_t y2,
446 Double_t z1,Double_t z2, Double_t c)
447{
448 //-----------------------------------------------------------------
449 // Initial approximation of the tangent of the track dip angle
450 //-----------------------------------------------------------------
451 //
452 x2-=x1;
453 y2-=y1;
454 z2-=z1;
455 Double_t d = TMath::Sqrt(x2*x2+y2*y2); // distance straight line
456 if (TMath::Abs(d*c*0.5)>1) return 0;
457 Double_t angle2 = TMath::ASin(d*c*0.5);
458 angle2 = z2*TMath::Abs(c/(angle2*2.));
459 return angle2;
460}
461
462
463Double_t AliTrackerBase::MakeTgl(Double_t x1,Double_t y1,
464 Double_t x2,Double_t y2,
465 Double_t z1,Double_t z2)
466{
467 //-----------------------------------------------------------------
468 // Initial approximation of the tangent of the track dip angle
469 //-----------------------------------------------------------------
470 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
471}
472
473
474AliExternalTrackParam * AliTrackerBase::MakeSeed( AliTrackPoint &point0, AliTrackPoint &point1, AliTrackPoint &point2){
475 //
476 // Make Seed - AliExternalTrackParam from input 3 points
477 // returning seed in local frame of point0
478 //
479 Double_t xyz0[3]={0,0,0};
480 Double_t xyz1[3]={0,0,0};
481 Double_t xyz2[3]={0,0,0};
482 Double_t alpha=point0.GetAngle();
483 Double_t xyz[3]={point0.GetX(),point0.GetY(),point0.GetZ()};
484 Double_t bxyz[3]; GetBxByBz(xyz,bxyz);
485 Double_t bz = bxyz[2];
486 //
487 // get points in frame of point 0
488 //
489 AliTrackPoint p0r = point0.Rotate(alpha);
490 AliTrackPoint p1r = point1.Rotate(alpha);
491 AliTrackPoint p2r = point2.Rotate(alpha);
492 xyz0[0]=p0r.GetX();
493 xyz0[1]=p0r.GetY();
494 xyz0[2]=p0r.GetZ();
495 xyz1[0]=p1r.GetX();
496 xyz1[1]=p1r.GetY();
497 xyz1[2]=p1r.GetZ();
498 xyz2[0]=p2r.GetX();
499 xyz2[1]=p2r.GetY();
500 xyz2[2]=p2r.GetZ();
501 //
502 // make covariance estimate
503 //
504 Double_t covar[15];
505 Double_t param[5]={0,0,0,0,0};
506 for (Int_t m=0; m<15; m++) covar[m]=0;
507 //
508 // calculate intitial param
509 param[0]=xyz0[1];
510 param[1]=xyz0[2];
511 param[2]=MakeSnp(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]);
512 param[4]=MakeC(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]);
513 param[3]=MakeTgl(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz0[2],xyz1[2],param[4]);
514
515 //covariance matrix - only diagonal elements
516 //Double_t dist=p0r.GetX()-p2r.GetX();
517 Double_t deltaP=0;
518 covar[0]= p0r.GetCov()[3];
519 covar[2]= p0r.GetCov()[5];
520 //sigma snp
521 deltaP= (MakeSnp(xyz0[0],xyz0[1]+TMath::Sqrt(p0r.GetCov()[3]),xyz1[0],xyz1[1],xyz2[0],xyz2[1])-param[2]);
522 covar[5]+= deltaP*deltaP;
523 deltaP= (MakeSnp(xyz0[0],xyz0[1],xyz1[0],xyz1[1]+TMath::Sqrt(p1r.GetCov()[3]),xyz2[0],xyz2[1])-param[2]);
524 covar[5]+= deltaP*deltaP;
525 deltaP= (MakeSnp(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]+TMath::Sqrt(p1r.GetCov()[3]))-param[2]);
526 covar[5]+= deltaP*deltaP;
527 //sigma tgl
528 //
529 deltaP=MakeTgl(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz0[2]+TMath::Sqrt(p1r.GetCov()[5]),xyz1[2],param[4])-param[3];
530 covar[9]+= deltaP*deltaP;
531 deltaP=MakeTgl(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz0[2],xyz1[2]+TMath::Sqrt(p1r.GetCov()[5]),param[4])-param[3];
532 covar[9]+= deltaP*deltaP;
533 //
534
535 deltaP=MakeC(xyz0[0],xyz0[1]+TMath::Sqrt(p0r.GetCov()[3]),xyz1[0],xyz1[1],xyz2[0],xyz2[1])-param[4];
536 covar[14]+= deltaP*deltaP;
537 deltaP=MakeC(xyz0[0],xyz0[1],xyz1[0],xyz1[1]+TMath::Sqrt(p1r.GetCov()[3]),xyz2[0],xyz2[1])-param[4];
538 covar[14]+= deltaP*deltaP;
539 deltaP=MakeC(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]+TMath::Sqrt(p2r.GetCov()[3]))-param[4];
540 covar[14]+= deltaP*deltaP;
541
542 covar[14]/=(bz*kB2C)*(bz*kB2C);
543 param[4]/=(bz*kB2C); // transform to 1/pt
544 AliExternalTrackParam * trackParam = new AliExternalTrackParam(xyz0[0],alpha,param, covar);
545 if (0) {
546 // consistency check -to put warnings here
547 // small disagrement once Track extrapolation used
548 // nice agreement in seeds with MC track parameters - problem in extrapoloation - to be fixed
549 // to check later
550 Double_t y1,y2,z1,z2;
551 trackParam->GetYAt(xyz1[0],bz,y1);
552 trackParam->GetZAt(xyz1[0],bz,z1);
553 trackParam->GetYAt(xyz2[0],bz,y2);
554 trackParam->GetZAt(xyz2[0],bz,z2);
555 if (TMath::Abs(y1-xyz1[1])> TMath::Sqrt(p1r.GetCov()[3]*5)){
556 AliWarningClass("Seeding problem y1\n");
557 }
558 if (TMath::Abs(y2-xyz2[1])> TMath::Sqrt(p2r.GetCov()[3]*5)){
559 AliWarningClass("Seeding problem y2\n");
560 }
561 if (TMath::Abs(z1-xyz1[2])> TMath::Sqrt(p1r.GetCov()[5]*5)){
562 AliWarningClass("Seeding problem z1\n");
563 }
564 }
565 return trackParam;
566}
567
568Double_t AliTrackerBase::FitTrack(AliExternalTrackParam * trackParam, AliTrackPointArray *pointArray, Double_t mass, Double_t maxStep){
569 //
570 // refit the track - trackParam using the points in point array
571 //
572 const Double_t kMaxSnp=0.99;
573 if (!trackParam) return 0;
574 Int_t npoints=pointArray->GetNPoints();
575 AliTrackPoint point,point2;
576 Double_t pointPos[2]={0,0};
577 Double_t pointCov[3]={0,0,0};
578 // choose coordinate frame
579 // in standard way the coordinate frame should be changed point by point
580 // Some problems with rotation observed
581 // rotate method of AliExternalTrackParam should be revisited
582 pointArray->GetPoint(point,0);
583 pointArray->GetPoint(point2,npoints-1);
584 Double_t alpha=TMath::ATan2(point.GetY()-point2.GetY(), point.GetX()-point2.GetX());
585
586 for (Int_t ipoint=npoints-1; ipoint>0; ipoint-=1){
587 pointArray->GetPoint(point,ipoint);
588 AliTrackPoint pr = point.Rotate(alpha);
589 trackParam->Rotate(alpha);
590 Bool_t status = PropagateTrackTo(trackParam,pr.GetX(),mass,maxStep,kFALSE,kMaxSnp);
591 if(!status){
592 AliWarningClass("Problem to propagate\n");
593 break;
594 }
595 if (TMath::Abs(trackParam->GetSnp())>kMaxSnp){
596 AliWarningClass("sin(phi) > kMaxSnp \n");
597 break;
598 }
599 pointPos[0]=pr.GetY();//local y
600 pointPos[1]=pr.GetZ();//local z
601 pointCov[0]=pr.GetCov()[3];//simay^2
602 pointCov[1]=pr.GetCov()[4];//sigmayz
603 pointCov[2]=pr.GetCov()[5];//sigmaz^2
604 trackParam->Update(pointPos,pointCov);
605 }
606 return 0;
607}