AliRoot Core  edcc906 (edcc906)
AliTrackerBase.cxx
Go to the documentation of this file.
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
20 // that is the base for the AliTracker class
21 // Origin: Marian.Ivanov@cern.ch
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"
30 #include "AliTrackPointArray.h"
31 #include "TVectorD.h"
32 
33 extern TGeoManager *gGeoManager;
34 
35 ClassImp(AliTrackerBase)
36 
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 //__________________________________________________________________________
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 //__________________________________________________________________________
72 {
73  AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
74  if (!fld) {
75  AliFatalClass("Field is not loaded");
76  //if (!fld)
77  return 0.5*kAlmost0Field;
78  }
79  Double_t bz = fld->SolenoidField();
80  return TMath::Sign(0.5*kAlmost0Field,bz) + bz;
81 }
82 
83 //__________________________________________________________________________
84 Double_t AliTrackerBase::GetBz(const Double_t *r) {
85  //------------------------------------------------------------------
86  // Returns Bz (kG) at the point "r" .
87  //------------------------------------------------------------------
88  AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
89  if (!fld) {
90  AliFatalClass("Field is not loaded");
91  // if (!fld)
92  return 0.5*kAlmost0Field;
93  }
94  Double_t bz = fld->GetBz(r);
95  return TMath::Sign(0.5*kAlmost0Field,bz) + bz;
96 }
97 
98 //__________________________________________________________________________
99 void 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) {
105  AliFatalClass("Field is not loaded");
106  // b[0] = b[1] = 0.;
107  // b[2] = 0.5*kAlmost0Field;
108  return;
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 
121 Double_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) {
154  AliFatalClass("No TGeo\n");
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  AliDebugClass(1,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  static int show_error = !(getenv("HLT_ONLINE_MODE") && strcmp(getenv("HLT_ONLINE_MODE"), "on") == 0);
220  if (show_error) AliErrorClass("Cannot cross boundary");
221  mparam[0] = bparam[0]/step;
222  mparam[1] = bparam[1];
223  mparam[2] = bparam[2]/step;
224  mparam[3] = bparam[3]/step;
225  mparam[5] = bparam[5]/step;
226  mparam[4] = step;
227  mparam[0] = 0.; // if crash of navigation take mean density 0
228  mparam[1] = 1000000; // and infinite rad length
229  return bparam[0]/step;
230  }
231  mparam[6]+=1.;
232  step += snext;
233  bparam[1] += snext/lparam[1];
234  bparam[2] += snext*lparam[2];
235  bparam[3] += snext*lparam[3];
236  bparam[5] += snext*lparam[5];
237  bparam[0] += snext*lparam[0];
238 
239  if (snext>=length) break;
240  if (!currentnode) break;
241  length -= snext;
242  material = currentnode->GetVolume()->GetMedium()->GetMaterial();
243  lparam[0] = material->GetDensity();
244  lparam[1] = material->GetRadLen();
245  lparam[2] = material->GetA();
246  lparam[3] = material->GetZ();
247  lparam[5] = lparam[3]/lparam[2];
248  if (material->IsMixture()) {
249  TGeoMixture * mixture = (TGeoMixture*)material;
250  lparam[5]=0;
251  Double_t sum =0;
252  for (Int_t iel=0;iel<mixture->GetNelements();iel++){
253  sum+= mixture->GetWmixt()[iel];
254  lparam[5]+= mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel];
255  }
256  lparam[5]/=sum;
257  }
258  gGeoManager->FindNextBoundaryAndStep(length, kFALSE);
259  snext = gGeoManager->GetStep();
260  }
261  mparam[0] = bparam[0]/step;
262  mparam[1] = bparam[1];
263  mparam[2] = bparam[2]/step;
264  mparam[3] = bparam[3]/step;
265  mparam[5] = bparam[5]/step;
266  return bparam[0]/step;
267 }
268 
269 
270 Bool_t
272  Double_t mass, Double_t maxStep, Bool_t rotateTo, Double_t maxSnp, Int_t sign, Bool_t addTimeStep, Bool_t correctMaterialBudget){
273  //----------------------------------------------------------------
274  //
275  // Propagates the track to the plane X=xk (cm) using the magnetic field map
276  // and correcting for the crossed material.
277  //
278  // mass - mass used in propagation - used for energy loss correction (if <0 then q = 2)
279  // maxStep - maximal step for propagation
280  //
281  // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
282  //
283  //----------------------------------------------------------------
284  const Double_t kEpsilon = 0.00001;
285  Double_t xpos = track->GetX();
286  Int_t dir = (xpos<xToGo) ? 1:-1;
287  //
288  while ( (xToGo-xpos)*dir > kEpsilon){
289  Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
290  Double_t x = xpos+step;
291  Double_t xyz0[3],xyz1[3],param[7];
292  track->GetXYZ(xyz0); //starting global position
293 
294  Double_t bz=GetBz(xyz0); // getting the local Bz
295  if (!track->GetXYZAt(x,bz,xyz1)) return kFALSE; // no prolongation
296  xyz1[2]+=kEpsilon; // waiting for bug correction in geo
297 
298  if (maxSnp>0 && TMath::Abs(track->GetSnpAt(x,bz)) >= maxSnp) return kFALSE;
299  if (!track->PropagateTo(x,bz)) return kFALSE;
300 
301  if (correctMaterialBudget){
302  MeanMaterialBudget(xyz0,xyz1,param);
303  Double_t xrho=param[0]*param[4], xx0=param[1];
304  if (sign) {if (sign<0) xrho = -xrho;} // sign is imposed
305  else { // determine automatically the sign from direction
306  if (dir>0) xrho = -xrho; // outward should be negative
307  }
308  //
309  if (!track->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE;
310  }
311 
312  if (rotateTo){
313  track->GetXYZ(xyz1); // global position
314  Double_t alphan = TMath::ATan2(xyz1[1], xyz1[0]);
315  if (maxSnp>0) {
316  if (TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE;
317 
318  //
319  Double_t ca=TMath::Cos(alphan-track->GetAlpha()), sa=TMath::Sin(alphan-track->GetAlpha());
320  Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf));
321  Double_t sinNew = sf*ca - cf*sa;
322  if (TMath::Abs(sinNew) >= maxSnp) return kFALSE;
323 
324  }
325  if (!track->AliExternalTrackParam::Rotate(alphan)) return kFALSE;
326 
327  }
328  xpos = track->GetX();
329  if (addTimeStep && track->IsStartedTimeIntegral()) {
330  if (!rotateTo) track->GetXYZ(xyz1); // if rotateTo==kTRUE, then xyz1 is already extracted
331  Double_t dX=xyz0[0]-xyz1[0],dY=xyz0[1]-xyz1[1],dZ=xyz0[2]-xyz1[2];
332  Double_t d=TMath::Sqrt(dX*dX + dY*dY + dZ*dZ);
333  if (sign) {if (sign>0) d = -d;} // step sign is imposed, positive means inward direction
334  else { // determine automatically the sign from direction
335  if (dir<0) d = -d;
336  }
337  track->AddTimeStep(d);
338  }
339  }
340  return kTRUE;
341 }
342 
343 Bool_t AliTrackerBase::PropagateTrackParamOnlyTo(AliExternalTrackParam *track, Double_t xToGo,Double_t maxStep, Bool_t rotateTo, Double_t maxSnp)
344 {
345  //----------------------------------------------------------------
346  //
347  // Propagates in fixed step size the track params ONLY to the plane X=xk (cm) using the magnetic field map
348  // W/O correcting for the crossed material.
349  // maxStep - maximal step for propagation
350  //
351  //----------------------------------------------------------------
352  const Double_t kEpsilon = 0.00001;
353  double xpos = track->GetX();
354  Int_t dir = (xpos<xToGo) ? 1:-1;
355  //
356  double xyz[3];
357  track->GetXYZ(xyz);
358  while ( (xToGo-xpos)*dir > kEpsilon){
359  Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
360  Double_t x = track->GetX()+step;
361  Double_t bz=GetBz(xyz); // getting the local Bz
362  if (!track->PropagateParamOnlyTo(x,bz)) return kFALSE;
363  track->GetXYZ(xyz); // global position
364  if (rotateTo){
365  Double_t alphan = TMath::ATan2(xyz[1], xyz[0]);
366  if (maxSnp>0 && TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE;
367  if (!track->AliExternalTrackParam::RotateParamOnly(alphan)) return kFALSE;
368  }
369  xpos = track->GetX();
370  }
371  return kTRUE;
372 }
373 
375  Double_t mass, Double_t maxStep, Bool_t rotateTo, Double_t maxSnp, Int_t sign, Bool_t addTimeStep, Bool_t correctMaterialBudget){
376  //----------------------------------------------------------------
377  //
378  // Propagates the track to the plane X=xk (cm) using the magnetic field map
379  // and correcting for the crossed material.
380  //
381  // mass - mass used in propagation - used for energy loss correction
382  // maxStep - maximal step for propagation
383  //
384  // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
385  //
386  //----------------------------------------------------------------
387  const Double_t kEpsilon = 0.00001;
388  Double_t xpos = track->GetX();
389  Int_t dir = (xpos<xToGo) ? 1:-1;
390  //
391  while ( (xToGo-xpos)*dir > kEpsilon){
392  Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
393  Double_t x = xpos+step;
394  Double_t xyz0[3],xyz1[3],param[7];
395  track->GetXYZ(xyz0); //starting global position
396 
397  Double_t bz=GetBz(xyz0); // getting the local Bz
398  if (!track->GetXYZAt(x,bz,xyz1)) return -1; // no prolongation
399  xyz1[2]+=kEpsilon; // waiting for bug correction in geo
400 
401  if (maxSnp>0 && TMath::Abs(track->GetSnpAt(x,bz)) >= maxSnp) return -2;
402  if (!track->PropagateTo(x,bz)) return -3;
403 
404  if (correctMaterialBudget){
405  MeanMaterialBudget(xyz0,xyz1,param);
406  Double_t xrho=param[0]*param[4], xx0=param[1];
407  if (sign) {if (sign<0) xrho = -xrho;} // sign is imposed
408  else { // determine automatically the sign from direction
409  if (dir>0) xrho = -xrho; // outward should be negative
410  }
411  //
412  if (!track->CorrectForMeanMaterial(xx0,xrho,mass)) return -4;
413  }
414 
415  if (rotateTo){
416  track->GetXYZ(xyz1); // global position
417  Double_t alphan = TMath::ATan2(xyz1[1], xyz1[0]);
418  if (maxSnp>0) {
419  if (TMath::Abs(track->GetSnp()) >= maxSnp) return -5;
420 
421  //
422  Double_t ca=TMath::Cos(alphan-track->GetAlpha()), sa=TMath::Sin(alphan-track->GetAlpha());
423  Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf));
424  Double_t sinNew = sf*ca - cf*sa;
425  if (TMath::Abs(sinNew) >= maxSnp) return -6;
426 
427  }
428  if (!track->AliExternalTrackParam::Rotate(alphan)) return -7;
429 
430  }
431  xpos = track->GetX();
432  if (addTimeStep && track->IsStartedTimeIntegral()) {
433  if (!rotateTo) track->GetXYZ(xyz1); // if rotateTo==kTRUE, then xyz1 is already extracted
434  Double_t dX=xyz0[0]-xyz1[0],dY=xyz0[1]-xyz1[1],dZ=xyz0[2]-xyz1[2];
435  Double_t d=TMath::Sqrt(dX*dX + dY*dY + dZ*dZ);
436  if (sign) {if (sign>0) d = -d;} // step sign is imposed, positive means inward direction
437  else { // determine automatically the sign from direction
438  if (dir<0) d = -d;
439  }
440  track->AddTimeStep(d);
441  }
442  }
443  return 1;
444 }
445 
446 Bool_t
448  Double_t xToGo,Double_t mass, Double_t maxStep, Bool_t rotateTo, Double_t maxSnp,Int_t sign, Bool_t addTimeStep,
449  Bool_t correctMaterialBudget){
450  //----------------------------------------------------------------
451  //
452  // Propagates the track to the plane X=xk (cm)
453  // taking into account all the three components of the magnetic field
454  // and correcting for the crossed material.
455  //
456  // mass - mass used in propagation - used for energy loss correction (if <0 then q=2)
457  // maxStep - maximal step for propagation
458  //
459  // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
460  //
461  //----------------------------------------------------------------
462  const Double_t kEpsilon = 0.00001;
463  Double_t xpos = track->GetX();
464  Int_t dir = (xpos<xToGo) ? 1:-1;
465  //
466  while ( (xToGo-xpos)*dir > kEpsilon){
467  Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
468  Double_t x = xpos+step;
469  Double_t xyz0[3],xyz1[3],param[7];
470  track->GetXYZ(xyz0); //starting global position
471 
472  Double_t b[3]; GetBxByBz(xyz0,b); // getting the local Bx, By and Bz
473 
474  if (!track->GetXYZAt(x,b[2],xyz1)) return kFALSE; // no prolongation
475  xyz1[2]+=kEpsilon; // waiting for bug correction in geo
476 
477  // if (maxSnp>0 && TMath::Abs(track->GetSnpAt(x,b[2])) >= maxSnp) return kFALSE;
478  if (!track->PropagateToBxByBz(x,b)) return kFALSE;
479  if (maxSnp>0 && TMath::Abs(track->GetSnp())>=maxSnp) return kFALSE;
480 
481  if (correctMaterialBudget) {
482  MeanMaterialBudget(xyz0,xyz1,param);
483  Double_t xrho=param[0]*param[4], xx0=param[1];
484  if (sign) {if (sign<0) xrho = -xrho;} // sign is imposed
485  else { // determine automatically the sign from direction
486  if (dir>0) xrho = -xrho; // outward should be negative
487  }
488  //
489  if (!track->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE;
490  }
491  if (rotateTo){
492  track->GetXYZ(xyz1); // global position
493  Double_t alphan = TMath::ATan2(xyz1[1], xyz1[0]);
494  /*
495  if (maxSnp>0) {
496  if (TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE;
497  Double_t ca=TMath::Cos(alphan-track->GetAlpha()), sa=TMath::Sin(alphan-track->GetAlpha());
498  Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf));
499  Double_t sinNew = sf*ca - cf*sa;
500  if (TMath::Abs(sinNew) >= maxSnp) return kFALSE;
501  }
502  */
503  if (!track->AliExternalTrackParam::Rotate(alphan)) return kFALSE;
504  if (maxSnp>0 && TMath::Abs(track->GetSnp())>=maxSnp) return kFALSE;
505  }
506  xpos = track->GetX();
507  if (addTimeStep && track->IsStartedTimeIntegral()) {
508  if (!rotateTo) track->GetXYZ(xyz1); // if rotateTo==kTRUE, then xyz1 is already extracted
509  Double_t dX=xyz0[0]-xyz1[0],dY=xyz0[1]-xyz1[1],dZ=xyz0[2]-xyz1[2];
510  Double_t d=TMath::Sqrt(dX*dX + dY*dY + dZ*dZ);
511  if (sign) {if (sign>0) d = -d;} // step sign is imposed, positive means inward direction
512  else { // determine automatically the sign from direction
513  if (dir<0) d = -d;
514  }
515  track->AddTimeStep(d);
516  }
517  }
518  return kTRUE;
519 }
520 
522  Double_t xToGo,Double_t maxStep, Bool_t rotateTo, Double_t maxSnp)
523 {
524  //----------------------------------------------------------------
525  //
526  // Propagates in fixed step size the track params ONLY to the plane X=xk (cm) using the magnetic field map
527  // W/O correcting for the crossed material.
528  // maxStep - maximal step for propagation
529  //
530  //----------------------------------------------------------------
531  const Double_t kEpsilon = 0.00001;
532  Double_t xpos = track->GetX();
533  Int_t dir = (xpos<xToGo) ? 1:-1;
534  Double_t xyz[3];
535  track->GetXYZ(xyz);
536  //
537  while ( (xToGo-xpos)*dir > kEpsilon){
538  Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
539  Double_t x = xpos+step;
540  Double_t b[3]; GetBxByBz(xyz,b); // getting the local Bx, By and Bz
541  if (!track->PropagateParamOnlyBxByBzTo(x,b)) return kFALSE;
542  if (maxSnp>0 && TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE;
543  track->GetXYZ(xyz);
544  if (rotateTo){
545  Double_t alphan = TMath::ATan2(xyz[1], xyz[0]);
546  if (!track->AliExternalTrackParam::Rotate(alphan)) return kFALSE;
547  }
548  xpos = track->GetX();
549  }
550  return kTRUE;
551 }
552 
554  Double_t mass, Double_t step,
555  const AliExternalTrackParam *backup) {
556  //
557  // This function brings the "track" with particle "mass" [GeV]
558  // to the same local coord. system and the same reference plane as
559  // of the "backup", doing it in "steps" [cm].
560  // Then, it calculates the 5D predicted Chi2 for these two tracks
561  //
562  Double_t chi2=kVeryBig;
563  Double_t alpha=backup->GetAlpha();
564  if (!track->Rotate(alpha)) return chi2;
565 
566  Double_t xb=backup->GetX();
567  Double_t sign=(xb < track->GetX()) ? 1. : -1.;
568  if (!PropagateTrackTo(track,xb,mass,step,kFALSE,kAlmost1,sign)) return chi2;
569 
570  chi2=track->GetPredictedChi2(backup);
571 
572  return chi2;
573 }
574 
575 
576 
577 
578 Double_t AliTrackerBase::MakeC(Double_t x1,Double_t y1,
579  Double_t x2,Double_t y2,
580  Double_t x3,Double_t y3)
581 {
582  //-----------------------------------------------------------------
583  // Initial approximation of the track curvature
584  //-----------------------------------------------------------------
585  x3 -=x1;
586  x2 -=x1;
587  y3 -=y1;
588  y2 -=y1;
589  //
590  Double_t det = x3*y2-x2*y3;
591  if (TMath::Abs(det)<1e-10) {
592  return 0;
593  }
594  //
595  Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
596  Double_t x0 = x3*0.5-y3*u;
597  Double_t y0 = y3*0.5+x3*u;
598  Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
599  if (det>0) c2*=-1;
600  return c2;
601 }
602 
603 Double_t AliTrackerBase::MakeSnp(Double_t x1,Double_t y1,
604  Double_t x2,Double_t y2,
605  Double_t x3,Double_t y3)
606 {
607  //-----------------------------------------------------------------
608  // Initial approximation of the track snp
609  //-----------------------------------------------------------------
610  x3 -=x1;
611  x2 -=x1;
612  y3 -=y1;
613  y2 -=y1;
614  //
615  Double_t det = x3*y2-x2*y3;
616  if (TMath::Abs(det)<1e-10) {
617  return 0;
618  }
619  //
620  Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
621  Double_t x0 = x3*0.5-y3*u;
622  Double_t y0 = y3*0.5+x3*u;
623  Double_t c2 = 1./TMath::Sqrt(x0*x0+y0*y0);
624  x0*=c2;
625  x0=TMath::Abs(x0);
626  if (y2*x2<0.) x0*=-1;
627  return x0;
628 }
629 
630 Double_t AliTrackerBase::MakeTgl(Double_t x1,Double_t y1,
631  Double_t x2,Double_t y2,
632  Double_t z1,Double_t z2, Double_t c)
633 {
634  //-----------------------------------------------------------------
635  // Initial approximation of the tangent of the track dip angle
636  //-----------------------------------------------------------------
637  //
638  const Double_t kEpsilon =0.00001;
639  x2-=x1;
640  y2-=y1;
641  z2-=z1;
642  Double_t d = TMath::Sqrt(x2*x2+y2*y2); // distance straight line
643  if (TMath::Abs(d*c*0.5)>1) return 0;
644  Double_t angle2 = TMath::ASin(d*c*0.5);
645  if (TMath::Abs(angle2)>kEpsilon) {
646  angle2 = z2*TMath::Abs(c/(angle2*2.));
647  }else{
648  angle2=z2/d;
649  }
650  return angle2;
651 }
652 
653 
654 Double_t AliTrackerBase::MakeTgl(Double_t x1,Double_t y1,
655  Double_t x2,Double_t y2,
656  Double_t z1,Double_t z2)
657 {
658  //-----------------------------------------------------------------
659  // Initial approximation of the tangent of the track dip angle
660  //-----------------------------------------------------------------
661  return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
662 }
663 
664 
666  //
667  // Make Seed - AliExternalTrackParam from input 3 points
668  // returning seed in local frame of point0
669  //
670  Double_t xyz0[3]={0,0,0};
671  Double_t xyz1[3]={0,0,0};
672  Double_t xyz2[3]={0,0,0};
673  Double_t alpha=point0.GetAngle();
674  Double_t xyz[3]={point0.GetX(),point0.GetY(),point0.GetZ()};
675  Double_t bxyz[3]; GetBxByBz(xyz,bxyz);
676  Double_t bz = bxyz[2];
677  //
678  // get points in frame of point 0
679  //
680  AliTrackPoint p0r = point0.Rotate(alpha);
681  AliTrackPoint p1r = point1.Rotate(alpha);
682  AliTrackPoint p2r = point2.Rotate(alpha);
683  xyz0[0]=p0r.GetX();
684  xyz0[1]=p0r.GetY();
685  xyz0[2]=p0r.GetZ();
686  xyz1[0]=p1r.GetX();
687  xyz1[1]=p1r.GetY();
688  xyz1[2]=p1r.GetZ();
689  xyz2[0]=p2r.GetX();
690  xyz2[1]=p2r.GetY();
691  xyz2[2]=p2r.GetZ();
692  //
693  // make covariance estimate
694  //
695  Double_t covar[15];
696  Double_t param[5]={0,0,0,0,0};
697  for (Int_t m=0; m<15; m++) covar[m]=0;
698  //
699  // calculate intitial param
700  param[0]=xyz0[1];
701  param[1]=xyz0[2];
702  param[2]=MakeSnp(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]);
703  param[4]=MakeC(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]);
704  param[3]=MakeTgl(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz0[2],xyz1[2],param[4]);
705 
706  //covariance matrix - only diagonal elements
707  //Double_t dist=p0r.GetX()-p2r.GetX();
708  Double_t deltaP=0;
709  covar[0]= p0r.GetCov()[3];
710  covar[2]= p0r.GetCov()[5];
711  //sigma snp
712  deltaP= (MakeSnp(xyz0[0],xyz0[1]+TMath::Sqrt(p0r.GetCov()[3]),xyz1[0],xyz1[1],xyz2[0],xyz2[1])-param[2]);
713  covar[5]+= deltaP*deltaP;
714  deltaP= (MakeSnp(xyz0[0],xyz0[1],xyz1[0],xyz1[1]+TMath::Sqrt(p1r.GetCov()[3]),xyz2[0],xyz2[1])-param[2]);
715  covar[5]+= deltaP*deltaP;
716  deltaP= (MakeSnp(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]+TMath::Sqrt(p1r.GetCov()[3]))-param[2]);
717  covar[5]+= deltaP*deltaP;
718  //sigma tgl
719  //
720  deltaP=MakeTgl(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz0[2]+TMath::Sqrt(p1r.GetCov()[5]),xyz1[2],param[4])-param[3];
721  covar[9]+= deltaP*deltaP;
722  deltaP=MakeTgl(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz0[2],xyz1[2]+TMath::Sqrt(p1r.GetCov()[5]),param[4])-param[3];
723  covar[9]+= deltaP*deltaP;
724  //
725 
726  deltaP=MakeC(xyz0[0],xyz0[1]+TMath::Sqrt(p0r.GetCov()[3]),xyz1[0],xyz1[1],xyz2[0],xyz2[1])-param[4];
727  covar[14]+= deltaP*deltaP;
728  deltaP=MakeC(xyz0[0],xyz0[1],xyz1[0],xyz1[1]+TMath::Sqrt(p1r.GetCov()[3]),xyz2[0],xyz2[1])-param[4];
729  covar[14]+= deltaP*deltaP;
730  deltaP=MakeC(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]+TMath::Sqrt(p2r.GetCov()[3]))-param[4];
731  covar[14]+= deltaP*deltaP;
732 
733  if (TMath::Abs(bz)>kAlmost0Field) {
734  covar[14]/=(bz*kB2C)*(bz*kB2C);
735  param[4]/=(bz*kB2C); // transform to 1/pt
736  }
737  else { // assign 0.6 GeV pT
738  const double kq2pt = 1./0.6;
739  param[4] = kq2pt;
740  covar[14] = (0.5*0.5)*kq2pt;
741  }
742  AliExternalTrackParam * trackParam = new AliExternalTrackParam(xyz0[0],alpha,param, covar);
743  if (0) {
744  // consistency check -to put warnings here
745  // small disagrement once Track extrapolation used
746  // nice agreement in seeds with MC track parameters - problem in extrapoloation - to be fixed
747  // to check later
748  Double_t y1,y2,z1,z2;
749  trackParam->GetYAt(xyz1[0],bz,y1);
750  trackParam->GetZAt(xyz1[0],bz,z1);
751  trackParam->GetYAt(xyz2[0],bz,y2);
752  trackParam->GetZAt(xyz2[0],bz,z2);
753  if (TMath::Abs(y1-xyz1[1])> TMath::Sqrt(p1r.GetCov()[3]*5)){
754  AliWarningClass("Seeding problem y1\n");
755  }
756  if (TMath::Abs(y2-xyz2[1])> TMath::Sqrt(p2r.GetCov()[3]*5)){
757  AliWarningClass("Seeding problem y2\n");
758  }
759  if (TMath::Abs(z1-xyz1[2])> TMath::Sqrt(p1r.GetCov()[5]*5)){
760  AliWarningClass("Seeding problem z1\n");
761  }
762  }
763  return trackParam;
764 }
765 
766 Double_t AliTrackerBase::FitTrack(AliExternalTrackParam * trackParam, AliTrackPointArray *pointArray, Double_t mass, Double_t maxStep){
767  //
768  // refit the track - trackParam using the points in point array
769  //
770  const Double_t kMaxSnp=0.99;
771  if (!trackParam) return 0;
772  Int_t npoints=pointArray->GetNPoints();
773  AliTrackPoint point,point2;
774  Double_t pointPos[2]={0,0};
775  Double_t pointCov[3]={0,0,0};
776  // choose coordinate frame
777  // in standard way the coordinate frame should be changed point by point
778  // Some problems with rotation observed
779  // rotate method of AliExternalTrackParam should be revisited
780  pointArray->GetPoint(point,0);
781  pointArray->GetPoint(point2,npoints-1);
782  Double_t alpha=TMath::ATan2(point.GetY()-point2.GetY(), point.GetX()-point2.GetX());
783 
784  for (Int_t ipoint=npoints-1; ipoint>0; ipoint-=1){
785  pointArray->GetPoint(point,ipoint);
786  AliTrackPoint pr = point.Rotate(alpha);
787  trackParam->Rotate(alpha);
788  Bool_t status = PropagateTrackTo(trackParam,pr.GetX(),mass,maxStep,kFALSE,kMaxSnp);
789  if(!status){
790  AliWarningClass("Problem to propagate\n");
791  break;
792  }
793  if (TMath::Abs(trackParam->GetSnp())>kMaxSnp){
794  AliWarningClass("sin(phi) > kMaxSnp \n");
795  break;
796  }
797  pointPos[0]=pr.GetY();//local y
798  pointPos[1]=pr.GetZ();//local z
799  pointCov[0]=pr.GetCov()[3];//simay^2
800  pointCov[1]=pr.GetCov()[4];//sigmayz
801  pointCov[2]=pr.GetCov()[5];//sigmaz^2
802  trackParam->Update(pointPos,pointCov);
803  }
804  return 0;
805 }
806 
807 
808 
810  //
811  // Update track 1 with track 2
812  //
813  //
814  //
815  TMatrixD vecXk(5,1); // X vector
816  TMatrixD covXk(5,5); // X covariance
817  TMatrixD matHk(5,5); // vector to mesurement
818  TMatrixD measR(5,5); // measurement error
819  TMatrixD vecZk(5,1); // measurement
820  //
821  TMatrixD vecYk(5,1); // Innovation or measurement residual
822  TMatrixD matHkT(5,5);
823  TMatrixD matSk(5,5); // Innovation (or residual) covariance
824  TMatrixD matKk(5,5); // Optimal Kalman gain
825  TMatrixD mat1(5,5); // update covariance matrix
826  TMatrixD covXk2(5,5); //
827  TMatrixD covOut(5,5);
828  //
829  Double_t *param1=(Double_t*) track1.GetParameter();
830  Double_t *covar1=(Double_t*) track1.GetCovariance();
831  Double_t *param2=(Double_t*) track2.GetParameter();
832  Double_t *covar2=(Double_t*) track2.GetCovariance();
833  //
834  // copy data to the matrix
835  for (Int_t ipar=0; ipar<5; ipar++){
836  for (Int_t jpar=0; jpar<5; jpar++){
837  covXk(ipar,jpar) = covar1[track1.GetIndex(ipar, jpar)];
838  measR(ipar,jpar) = covar2[track2.GetIndex(ipar, jpar)];
839  matHk(ipar,jpar)=0;
840  mat1(ipar,jpar)=0;
841  }
842  vecXk(ipar,0) = param1[ipar];
843  vecZk(ipar,0) = param2[ipar];
844  matHk(ipar,ipar)=1;
845  mat1(ipar,ipar)=1;
846  }
847  //
848  //
849  //
850  //
851  //
852  vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
853  matHkT=matHk.T(); matHk.T();
854  matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance
855  matSk.Invert();
856  matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
857  vecXk += matKk*vecYk; // updated vector
858  covXk2 = (mat1-(matKk*matHk));
859  covOut = covXk2*covXk;
860  //
861  //
862  //
863  // copy from matrix to parameters
864  if (0) {
865  vecXk.Print();
866  vecZk.Print();
867  //
868  measR.Print();
869  covXk.Print();
870  covOut.Print();
871  //
872  track1.Print();
873  track2.Print();
874  }
875 
876  for (Int_t ipar=0; ipar<5; ipar++){
877  param1[ipar]= vecXk(ipar,0) ;
878  for (Int_t jpar=0; jpar<5; jpar++){
879  covar1[track1.GetIndex(ipar, jpar)]=covOut(ipar,jpar);
880  }
881  }
882 }
883 
884 
TBrowser b
Definition: RunAnaESD.C:12
static Double_t GetBz()
const Float_t * GetCov() const
const Double_t kEpsilon
void Print(Option_t *option="") const
Bool_t Rotate(Double_t alpha)
Double_t GetBz(const Double_t *xyz) const
Definition: AliMagF.cxx:296
Bool_t Update(const Double_t p[2], const Double_t cov[3])
Bool_t GetYAt(Double_t x, Double_t b, Double_t &y) const
Bool_t GetXYZ(Double_t *p) const
static Int_t GetIndex(Int_t i, Int_t j)
Float_t GetAngle() const
void dZ()
Definition: CalibAlign.C:517
Bool_t IsUniform() const
Definition: AliMagF.h:57
virtual Bool_t IsStartedTimeIntegral() const
virtual void Field(const Double_t *x, Double_t *b)
Definition: AliMagF.cxx:280
Int_t GetNPoints() const
Double_t GetAlpha() const
AliTPCfastTrack * track
static Double_t MakeTgl(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Double_t z1, Double_t z2)
Bool_t PropagateParamOnlyTo(Double_t xk, Double_t b)
static Double_t MakeSnp(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Double_t x3, Double_t y3)
Double_t SolenoidField() const
Definition: AliMagF.h:67
npoints
Definition: driftITSTPC.C:85
#define AliFatalClass(message)
Definition: AliLog.h:645
const Double_t * GetParameter() const
static Double_t MeanMaterialBudget(const Double_t *start, const Double_t *end, Double_t *mparam)
#define AliWarning(message)
Definition: AliLog.h:541
static void UpdateTrack(AliExternalTrackParam &, const AliExternalTrackParam &)
void sum()
static Bool_t PropagateTrackParamOnlyToBxByBz(AliExternalTrackParam *track, Double_t xToGo, Double_t maxStep, Bool_t rotateTo=kTRUE, Double_t maxSnp=0.8)
void dY()
Definition: CalibAlign.C:547
#define AliErrorClass(message)
Definition: AliLog.h:596
const Double_t kB2C
Definition: AliVParticle.h:25
Float_t GetY() const
static Bool_t PropagateTrackParamOnlyTo(AliExternalTrackParam *track, Double_t xToGo, Double_t maxStep, Bool_t rotateTo=kTRUE, Double_t maxSnp=0.8)
#define AliDebugClass(logLevel, message)
Definition: AliLog.h:313
Double_t chi2
Definition: AnalyzeLaser.C:7
Bool_t GetPoint(AliTrackPoint &p, Int_t i) const
static Int_t PropagateTrackTo2(AliExternalTrackParam *track, Double_t x, Double_t m, Double_t maxStep, Bool_t rotateTo=kTRUE, Double_t maxSnp=0.8, Int_t sign=0, Bool_t addTimeStep=kFALSE, Bool_t correctMaterialBudget=kTRUE)
Float_t GetZ() const
TGeoManager * gGeoManager
#define AliWarningClass(message)
Definition: AliLog.h:546
const Double_t * GetCovariance() const
static void GetBxByBz(const Double_t r[3], Double_t b[3])
const Double_t kVeryBig
const Double_t kAlmost0Field
Definition: AliVParticle.h:26
static Double_t FitTrack(AliExternalTrackParam *trackParam, AliTrackPointArray *pointArray, Double_t mass, Double_t maxStep)
virtual void AddTimeStep(Double_t)
static Bool_t PropagateTrackTo(AliExternalTrackParam *track, Double_t x, Double_t m, Double_t maxStep, Bool_t rotateTo=kTRUE, Double_t maxSnp=0.8, Int_t sign=0, Bool_t addTimeStep=kFALSE, Bool_t correctMaterialBudget=kTRUE)
Bool_t PropagateParamOnlyBxByBzTo(Double_t xk, const Double_t b[3])
static Bool_t PropagateTrackToBxByBz(AliExternalTrackParam *track, Double_t x, Double_t m, Double_t maxStep, Bool_t rotateTo=kTRUE, Double_t maxSnp=0.8, Int_t sign=0, Bool_t addTimeStep=kFALSE, Bool_t correctMaterialBudget=kTRUE)
Bool_t GetZAt(Double_t x, Double_t b, Double_t &z) const
static Double_t MakeC(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Double_t x3, Double_t y3)
Bool_t CorrectForMeanMaterial(Double_t xOverX0, Double_t xTimesRho, Double_t mass, Bool_t anglecorr=kFALSE, Double_t(*f)(Double_t)=AliExternalTrackParam::BetheBlochSolid)
static Double_t GetTrackPredictedChi2(AliExternalTrackParam *track, Double_t mass, Double_t step, const AliExternalTrackParam *backup)
Double_t GetPredictedChi2(const Double_t p[2], const Double_t cov[3]) const
Float_t GetX() const
static AliExternalTrackParam * MakeSeed(AliTrackPoint &point0, AliTrackPoint &point1, AliTrackPoint &point2)
Bool_t PropagateToBxByBz(Double_t x, const Double_t b[3])
Bool_t PropagateTo(Double_t p[3], Double_t covyz[3], Double_t covxyz[3], Double_t b)
Double_t GetSnpAt(Double_t x, Double_t b) const
AliTrackPoint & Rotate(Float_t alpha) const
class TMatrixT< Double_t > TMatrixD
Bool_t GetXYZAt(Double_t x, Double_t b, Double_t r[3]) const
const Double_t kAlmost1
Definition: AliVParticle.h:22