AliRoot Core  ee782a0 (ee782a0)
AliMFTTrackExtrap.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$ */
17 
18 //-----------------------------------------------------------------------------
19 // Class AliMFTTrackExtrap
20 // ------------------------
21 // Tools for track extrapolation in ALICE MFT
22 // Adapted from AliMUONTrackExtrap by R. Tieulent
23 // Original Author: Philippe Pillot
24 //-----------------------------------------------------------------------------
25 
26 #include "AliMFTTrackExtrap.h"
27 #include "AliMFTTrackParam.h"
28 #include "AliMFTConstants.h"
29 
30 #include "AliMagF.h"
31 #include "AliExternalTrackParam.h"
32 
33 #include <TGeoGlobalMagField.h>
34 #include <TGeoManager.h>
35 #include <TMath.h>
36 #include <TDatabasePDG.h>
37 
38 #include <Riostream.h>
39 
40 using std::endl;
41 using std::cout;
43 ClassImp(AliMFTTrackExtrap) // Class implementation in ROOT context
45 
46 const Double_t AliMFTTrackExtrap::fgkSimpleBPosition = 0.5 * (AliMFTConstants::DefaultPlaneZ(0) + AliMFTConstants::DefaultPlaneZ(9));
47 //const Double_t AliMFTTrackExtrap::fgkSimpleBLength = 0.5 * (AliMUONConstants::CoilL() + AliMUONConstants::YokeL());
48  Double_t AliMFTTrackExtrap::fgSimpleBValue = 0.;
49  Bool_t AliMFTTrackExtrap::fgFieldON = kFALSE;
50 const Bool_t AliMFTTrackExtrap::fgkUseHelix = kTRUE;
51 const Int_t AliMFTTrackExtrap::fgkMaxStepNumber = 5000;
52 const Double_t AliMFTTrackExtrap::fgkHelixStepLength = 0.1;
53 const Double_t AliMFTTrackExtrap::fgkRungeKuttaMaxResidue = 0.002;
54 
55 //__________________________________________________________________________
56 void AliMFTTrackExtrap::SetField()
57 {
60  const Double_t x[3] = {0.,0.,fgkSimpleBPosition};
61  Double_t b[3] = {0.,0.,0.};
62  TGeoGlobalMagField::Instance()->Field(x,b);
63  cout<<Form("Field = %e %e %e",b[0],b[1],b[2])<<endl;
64  fgSimpleBValue = b[2];
65  fgFieldON = (TMath::Abs(fgSimpleBValue) > 1.e-10) ? kTRUE : kFALSE;
66 // fgFieldON = kFALSE;
67 }
68 
69 //__________________________________________________________________________
70 Double_t AliMFTTrackExtrap::GetImpactParamFromBendingMomentum(Double_t bendingMomentum)
71 {
76 
77 // if (bendingMomentum == 0.) return 1.e10;
78 //
79 // const Double_t kCorrectionFactor = 1.1; // impact parameter is 10% underestimated
80 //
81 // return kCorrectionFactor * (-0.0003 * fgSimpleBValue * fgkSimpleBLength * fgkSimpleBPosition / bendingMomentum);
82 }
83 
84 //__________________________________________________________________________
85 Double_t
87 {
92 
93 // if (impactParam == 0.) return 1.e10;
94 //
95 // const Double_t kCorrectionFactor = 1.1; // bending momentum is 10% underestimated
96 //
97 // if (fgFieldON)
98 // {
99 // return kCorrectionFactor * (-0.0003 * fgSimpleBValue * fgkSimpleBLength * fgkSimpleBPosition / impactParam);
100 // }
101 // else
102 // {
103 // return AliMUONConstants::GetMostProbBendingMomentum();
104 // }
105 }
106 //__________________________________________________________________________
107 Double_t AliMFTTrackExtrap::Sagitta(Int_t nVal, Double_t *xVal, Double_t *yVal, Double_t &distL2, Double_t &q2)
108 {
111  cout<<"Sagitta"<<endl;
112  Double_t q0,q1,chi2;
113  // fit by a polynom of 2nd order
114  chi2 = QuadraticRegression(nVal,xVal ,yVal, q0, q1,q2);
115  cout<<Form("q param = %f %f %f chi2 = %e ",q0, q1,q2,chi2)<<endl;
116  cout<<Form("pt from 2nd order parameter %f ",0.01/q2/2.*0.3*fgSimpleBValue/10.)<<endl;
117 
118 
119  Double_t maxDist = 0.;
120  Int_t i1 = -1;
121  Int_t i2 = -1;
122  Double_t sagitta = 0.;
123  Double_t dist ;
124  Double_t distTot = 0.;
125 
126  for (int i=0; i<nVal-1; i++) {
127  distTot += TMath::Sqrt((xVal[i]-xVal[i+1]) * (xVal[i]-xVal[i+1]) + (yVal[i]-yVal[i+1]) * (yVal[i]-yVal[i+1]));
128 
129  for (int j = i+1; j<nVal; j++) {
130  Double_t dist = (xVal[i]-xVal[j]) * (xVal[i]-xVal[j]) + (yVal[i]-yVal[j]) * (yVal[i]-yVal[j]);
131  if(dist>maxDist){
132  i1 = i;
133  i2 = j;
134  maxDist = dist;
135  }
136  }
137  }
138 
139  cout<<Form("distMAx = %f distTot =%f i1 = %d i2 =%d",TMath::Sqrt(maxDist),distTot,i1,i2)<<endl;
140 
141  distL2 = 1.e-2*TMath::Sqrt((xVal[i1]-xVal[i2]) * (xVal[i1]-xVal[i2]) + (yVal[i1]-yVal[i2]) * (yVal[i1]-yVal[i2]) );
142 
143  Double_t p1 = (yVal[i1]-yVal[i2]) / (xVal[i1]-xVal[i2]);
144  Double_t p0 = yVal[i2] - xVal[i2] * p1;
145  cout<<Form("p param = %f %f ",p0, p1)<<endl;
146  q2 = TMath::Sign(q2,q1*q2*(-(p0 + p1 * (xVal[i1]+xVal[i2])/2.))*(fgSimpleBValue));
147 
148 
149  Double_t y12 = q0+q1*xVal[i1]+q2*xVal[i1]*xVal[i1];
150  Double_t y22 = q0+q1*xVal[i2]+q2*xVal[i2]*xVal[i2];
151  cout<<Form("x1, y1 %f %f ",xVal[i1], y12)<<endl;
152  cout<<Form("x2, y2 %f %f ",xVal[i2], y22)<<endl;
153 
154 // Double_t p12 = (y12-y22) / (xVal[i1]-xVal[i2]);
155 // Double_t p02 = y22 - xVal[i2] * p12;
156 // cout<<Form("p2 param = %f %f ",p02, p12)<<endl;
157 
158  Double_t maxSagitta = 0.;
159  for (int i=0; i<20; i++) {
160  Double_t step =TMath::Abs(xVal[i1]-xVal[i2])/20.;
161  Double_t xmiddle = xVal[i1] + i*step;
162 
163  Double_t y1 = p0 + p1 * xmiddle;
164  Double_t p1perp = -1./p1;
165  Double_t p0perp = p0 + xmiddle *(p1-p1perp);
166  //cout<<Form("p param perp = %f %f ",p0perp, p1perp)<<endl;
167 
168  Double_t xmax = (p1perp - q1 + TMath::Sqrt(p1perp*p1perp - 2*p1perp*q1 + q1*q1 + 4*p0perp*q2 - 4*q0*q2))/(2*q2);
169  Double_t xmax2 = -(q1 -p1perp + TMath::Sqrt(p1perp*p1perp - 2*p1perp*q1 + q1*q1 + 4*p0perp*q2 - 4*q0*q2))/(2*q2);
170  // cout<<Form("xmax = %f xmax2 %f",xmax,xmax2)<<endl;
171 
172  if (TMath::Abs(xmax2-xmiddle) < TMath::Abs(xmax-xmiddle)) xmax = xmax2;
173 
174 
175  Double_t y2 = q0 + q1 * xmax + q2 * xmax * xmax;
176 
177  sagitta = 1e-2*TMath::Sqrt((xmiddle-xmax) * (xmiddle-xmax) + (y1-y2) * (y1-y2) );
178 
179  // cout<<Form("sagitta = %f Bz = %f",sagitta,fgSimpleBValue)<<endl;
180 
181  sagitta = TMath::Sign(sagitta,q1*q2*(-y1)*(fgSimpleBValue));
182 
183  if(TMath::Abs(sagitta)>TMath::Abs(maxSagitta)) maxSagitta = sagitta;
184 
185  }
186  cout<<Form(" Max sagitta = %e => pt = %f",maxSagitta, 0.3*0.5*distL2*distL2/8./maxSagitta)<<endl;
187 // p0 = p02;
188 // p1 = p12;
189 //
190 //
191 // maxSagitta = 0.;
192 // for (int i=0; i<20; i++) {
193 // Double_t step =TMath::Abs(xVal[i1]-xVal[i2])/20.;
194 // Double_t xmiddle = xVal[i1] + i*step;
196 //
197 // Double_t y1 = p0 + p1 * xmiddle;
198 // Double_t p1perp = -1./p1;
199 // Double_t p0perp = p0 + xmiddle *(p1-p1perp);
200 // //cout<<Form("p param perp = %f %f ",p0perp, p1perp)<<endl;
201 //
202 // Double_t xmax = (p1perp - q1 + TMath::Sqrt(p1perp*p1perp - 2*p1perp*q1 + q1*q1 + 4*p0perp*q2 - 4*q0*q2))/(2*q2);
203 // Double_t xmax2 = -(q1 - p1perp + TMath::Sqrt(p1perp*p1perp - 2*p1perp*q1 + q1*q1 + 4*p0perp*q2 - 4*q0*q2))/(2*q2);
205 //
206 // if (TMath::Abs(xmax2-xmiddle) < TMath::Abs(xmax-xmiddle)) xmax = xmax2;
207 //
208 //
209 // Double_t y2 = q0 + q1 * xmax + q2 * xmax * xmax;
210 //
211 // sagitta = 1e-2*TMath::Sqrt((xmiddle-xmax) * (xmiddle-xmax) + (y1-y2) * (y1-y2) );
212 //
214 //
215 // sagitta = TMath::Sign(sagitta,q1*q2*(-y1)*(fgSimpleBValue));
216 //
218 // if(sagitta>maxSagitta) maxSagitta = sagitta;
219 //
220 // }
221 // cout<<Form(" Max sagitta 2 = %e => pt = %f",maxSagitta, 0.3*0.5*distL2*distL2/8./maxSagitta )<<endl;
222 
223  return maxSagitta;
224 }
225 
226 //__________________________________________________________________________
227 Double_t AliMFTTrackExtrap::LinearRegression(Int_t nVal, Double_t *xVal, Double_t *yVal, Double_t &p0, Double_t &p1)
228 {
231  Double_t meanx =0, meany=0.;
232 
233  for (Int_t i = 0; i< nVal; i++) {
234  meanx =(meanx*i+xVal[i])/(i+1);
235  meany =(meany*i+yVal[i])/(i+1);
236  }
237  Double_t cov_xy = 0., var_x=0., var_y=0.;
238  for (Int_t i = 0; i< nVal; i++) {
239  var_x += (xVal[i] -meanx) * (xVal[i] -meanx);
240  var_y += (yVal[i] -meany) * (yVal[i] -meany);
241  cov_xy += (xVal[i] -meanx) * (yVal[i] -meany);
242  }
243  if(var_x<1.e-6) return 0.;
244  p1 = cov_xy/var_x;
245  p0 = meany - p1 * meanx;
246 
247  Double_t chi2 = 0.;
248  Double_t yest=0.;
249  for (Int_t i = 0; i< nVal; i++) {
250  yest = xVal[i]*p1+p0;
251  chi2 += (yest-yVal[i]) * (yest-yVal[i]);
252  }
253  chi2 /= var_y;
254 
255 
256  return chi2;
257 }
258 //__________________________________________________________________________
259 Double_t AliMFTTrackExtrap::QuadraticRegression(Int_t nVal, Double_t *xVal, Double_t *yVal, Double_t &p0, Double_t &p1, Double_t &p2)
260 {
264 
265  TMatrixD y(nVal,1);
266  TMatrixD x(nVal,3);
267  TMatrixD xtrans(3,nVal);
268 
269  for (int i=0; i<nVal; i++) {
270  y(i,0) = yVal[i];
271  x(i,0) = 1.;
272  x(i,1) = xVal[i];
273  x(i,2) = xVal[i]*xVal[i];
274  xtrans(0,i) = 1.;
275  xtrans(1,i) = xVal[i];
276  xtrans(2,i) = xVal[i]*xVal[i];
277  }
278  TMatrixD tmp(xtrans,TMatrixD::kMult,x);
279  tmp.Invert();
280 
281  TMatrixD tmp2(xtrans,TMatrixD::kMult,y);
282  TMatrixD b(tmp,TMatrixD::kMult,tmp2);
283 
284  p0 = b(0,0);
285  p1 = b(1,0);
286  p2 = b(2,0);
287 
288  // chi2 = (y-xb)^t . W . (y-xb)
289  TMatrixD tmp3(x,TMatrixD::kMult,b);
290  TMatrixD tmp4(y,TMatrixD::kMinus,tmp3);
291  TMatrixD chi2(tmp4,TMatrixD::kTransposeMult,tmp4);
292 
293 
294  return chi2(0,0);
295 }
296 //__________________________________________________________________________
297 Double_t AliMFTTrackExtrap::CircleRegression(Int_t nVal, Double_t *xVal, Double_t *yVal)
298 {
302  Double_t sumxi2 =0., sumxi =0., sumxiyi =0., sumyi =0.,sumyi2 =0., sumxi3 =0., sumyi3 =0.;
303  Double_t sumxi2yi=0., sumxiyi2=0.;
304  Double_t xi,yi, ri;
305  for (int i=0; i<nVal; i++) {
306  xi = xVal[i]/100.;
307  yi = yVal[i]/100.;
308  ri = TMath::Sqrt(xi*xi + yi*yi);
309 // xi /= ri*ri;
310 // yi /= ri*ri;
311 
312  sumxi += xi;
313  sumyi += yi;
314  sumxi2 += xi*xi;
315  sumyi2 += yi*yi;
316  sumxi3 += xi*xi*xi;
317  sumyi3 += yi*yi*yi;
318  sumxiyi += xi*yi;
319  sumxi2yi += xi*xi*yi;
320  sumxiyi2 += xi*yi*yi;
321  }
322 
323  Double_t A = nVal * sumxi2 - sumxi*sumxi;
324  Double_t B = nVal * sumxiyi - sumxi*sumyi;
325  Double_t C = nVal * sumyi2 - sumyi*sumyi;
326  Double_t D = 0.5*(nVal*sumxiyi2 -sumxi*sumyi2 +nVal*sumxi3 -sumxi*sumxi2);
327  Double_t E = 0.5*(nVal*sumxi2yi -sumyi*sumxi2 +nVal*sumyi3 -sumyi*sumyi2);
328 
329  Double_t aM = (D*C-B*E) / (A*C-B*B);
330  Double_t bM = (A*E-B*D) / (A*C-B*B);
331 
332  Double_t rM = 0.;
333  Double_t rK = 0.;
334 
335  for (int i=0; i<nVal; i++) {
336  xi = xVal[i]/100.;
337  yi = yVal[i]/100.;
338 
339  rM += TMath::Sqrt( (xi-aM)*(xi-aM) + (yi-bM)*(yi-bM) );
340  rK += ((xi-aM)*(xi-aM) + (yi-bM)*(yi-bM) );
341  }
342  rM /= nVal;
343  rK = TMath::Sqrt(rK/nVal);
344 
345  cout<<Form("aM %f bM %f rM %f rK %f => pt = %f or %f ",aM,bM,rM,rK,rM*0.3*0.5, rK*0.3*0.5)<<endl;
346 
347  return (rM+rK)/2.;
348 }
349 
350 //__________________________________________________________________________
351 void AliMFTTrackExtrap::LinearExtrapToZ(AliMFTTrackParam* trackParam, Double_t zEnd)
352 {
355 
356  if (trackParam->GetZ() == zEnd) return; // nothing to be done if same z
357 
358  // Compute track parameters
359  Double_t dZ = zEnd - trackParam->GetZ();
360 
361  trackParam->SetX(trackParam->GetX() + trackParam->GetSlopeX() * dZ);
362  trackParam->SetY(trackParam->GetY() + trackParam->GetSlopeY() * dZ);
363  trackParam->SetZ(zEnd);
364 }
365 
366 //__________________________________________________________________________
367 void AliMFTTrackExtrap::LinearExtrapToZCov(AliMFTTrackParam* trackParam, Double_t zEnd, Bool_t updatePropagator)
368 {
371 
372  if (trackParam->GetZ() == zEnd) return; // nothing to be done if same z
373 
374  // No need to propagate the covariance matrix if it does not exist
375  if (!trackParam->CovariancesExist()) {
376  cout<<"W-AliMUONTrackExtrap::LinearExtrapToZCov: Covariance matrix does not exist"<<endl;
377  // Extrapolate linearly track parameters to "zEnd"
378  LinearExtrapToZ(trackParam,zEnd);
379  return;
380  }
381 
382  // Compute track parameters
383  Double_t dZ = zEnd - trackParam->GetZ();
384  trackParam->SetX(trackParam->GetX() + trackParam->GetSlopeX() * dZ);
385  trackParam->SetY(trackParam->GetY() + trackParam->GetSlopeY() * dZ);
386  trackParam->SetZ(zEnd);
387 
388  // Calculate the jacobian related to the track parameters linear extrapolation to "zEnd"
389  TMatrixD jacob(5,5);
390  jacob.UnitMatrix();
391  jacob(0,2) = dZ;
392  jacob(1,3) = dZ;
393 
394  // Extrapolate track parameter covariances to "zEnd"
395  TMatrixD tmp(trackParam->GetCovariances(),TMatrixD::kMultTranspose,jacob);
396  TMatrixD tmp2(jacob,TMatrixD::kMult,tmp);
397  trackParam->SetCovariances(tmp2);
398 
399  // Update the propagator if required
400  if (updatePropagator) trackParam->UpdatePropagator(jacob);
401 }
402 
403 //__________________________________________________________________________
404 Bool_t AliMFTTrackExtrap::ExtrapToZ(AliMFTTrackParam* trackParam, Double_t zEnd)
405 {
408 //return AliMFTTrackExtrap::ExtrapToZHelix(trackParam,zEnd);
409  if (!fgFieldON) {
410  AliMFTTrackExtrap::LinearExtrapToZ(trackParam,zEnd);
411  return kTRUE;
412  }
413  else if (fgkUseHelix) return AliMFTTrackExtrap::ExtrapToZHelix(trackParam,zEnd);
414  else return AliMFTTrackExtrap::ExtrapToZRungekutta(trackParam,zEnd);
415 }
416 
417 //__________________________________________________________________________
418 Bool_t AliMFTTrackExtrap::ExtrapToZHelix(AliMFTTrackParam* trackParam, Double_t zEnd)
419 {
420 // cout<<"I-AliMFTTrackExtrap::ExtrapToZHelix: Entering ------ "<<endl;
423  if (trackParam->GetZ() == zEnd) return kTRUE; // nothing to be done if same Z
424  Double_t forwardBackward; // +1 if forward, -1 if backward
425  if (zEnd < trackParam->GetZ()) forwardBackward = 1.0; // spectro. z<0
426  else forwardBackward = -1.0;
427  Double_t v3[7], v3New[7]; // 7 in parameter ????
428  Int_t i3, stepNumber;
429  // For safety: return kTRUE or kFALSE ????
430  // Parameter vector for calling EXTRAP_ONESTEP
431  ConvertTrackParamForExtrap(trackParam, forwardBackward, v3);
432  // sign of charge (sign of fInverseBendingMomentum if forward motion)
433  // must be changed if backward extrapolation
434  Double_t chargeExtrap = forwardBackward * TMath::Sign(Double_t(1.0), trackParam->GetInverseTransverseMomentum());
435 // cout<<"chargeExtrap = "<<chargeExtrap<<endl;
436 
437  // Extrapolation loop
438  stepNumber = 0;
439 // cout<<" (-forwardBackward * (v3[2] - zEnd)) = "<<(-forwardBackward * (v3[2] - zEnd))<<endl;
440  while (((-forwardBackward * (v3[2] - zEnd)) <= 0.0) && (stepNumber < fgkMaxStepNumber)) { // spectro. z<0
441  stepNumber++;
442  ExtrapOneStepHelix(chargeExtrap, fgkHelixStepLength, v3, v3New);
443 // cout<<" (-forwardBackward * (v3New[2] - zEnd)) = "<<(-forwardBackward * (v3New[2] - zEnd))<<endl;
444 
445  if ((-forwardBackward * (v3New[2] - zEnd)) > 0.0) break; // one is beyond Z spectro. z<0
446  // better use TArray ????
447  for (i3 = 0; i3 < 7; i3++) {
448 // cout<<" v3New["<<i3<< "] = "<<v3New[i3]<<endl;
449  v3[i3] = v3New[i3];
450  }
451  }
452  // check fgkMaxStepNumber ????
453  // Interpolation back to exact Z (2nd order)
454  // should be in function ???? using TArray ????
455  Double_t dZ12 = v3New[2] - v3[2]; // 1->2
456  if (TMath::Abs(dZ12) > 0) {
457  Double_t dZ1i = zEnd - v3[2]; // 1-i
458  Double_t dZi2 = v3New[2] - zEnd; // i->2
459  Double_t xPrime = (v3New[0] - v3[0]) / dZ12;
460  Double_t xSecond = ((v3New[3] / v3New[5]) - (v3[3] / v3[5])) / dZ12;
461  Double_t yPrime = (v3New[1] - v3[1]) / dZ12;
462  Double_t ySecond = ((v3New[4] / v3New[5]) - (v3[4] / v3[5])) / dZ12;
463  v3[0] = v3[0] + xPrime * dZ1i - 0.5 * xSecond * dZ1i * dZi2; // X
464  v3[1] = v3[1] + yPrime * dZ1i - 0.5 * ySecond * dZ1i * dZi2; // Y
465  v3[2] = zEnd; // Z
466  Double_t xPrimeI = xPrime - 0.5 * xSecond * (dZi2 - dZ1i);
467  Double_t yPrimeI = yPrime - 0.5 * ySecond * (dZi2 - dZ1i);
468  // (PX, PY, PZ)/PTOT assuming forward motion
469  v3[5] = 1.0 / TMath::Sqrt(1.0 + xPrimeI * xPrimeI + yPrimeI * yPrimeI); // PZ/PTOT
470  v3[3] = xPrimeI * v3[5]; // PX/PTOT
471  v3[4] = yPrimeI * v3[5]; // PY/PTOT
472  } else {
473  cout<<"W-AliMFTTrackExtrap::ExtrapToZHelix: Extrap. to Z not reached, Z = "<<zEnd<<endl;
474  }
475  // Recover track parameters (charge back for forward motion)
476  RecoverTrackParam(v3, chargeExtrap * forwardBackward, trackParam);
477  return kTRUE;
478 }
479 
480 //__________________________________________________________________________
482 {
485  if (trackParam->GetZ() == zEnd) return kTRUE; // nothing to be done if same Z
486  Double_t forwardBackward; // +1 if forward, -1 if backward
487  if (zEnd < trackParam->GetZ()) forwardBackward = 1.0; // spectro. z<0
488  else forwardBackward = -1.0;
489  // sign of charge (sign of fInverseBendingMomentum if forward motion)
490  // must be changed if backward extrapolation
491  Double_t chargeExtrap = forwardBackward * TMath::Sign(Double_t(1.0), trackParam->GetInverseTransverseMomentum());
492  Double_t v3[7], v3New[7];
493  Double_t dZ, step;
494  Int_t stepNumber = 0;
495 
496  // Extrapolation loop (until within tolerance or the track turn around)
497  Double_t residue = zEnd - trackParam->GetZ();
498  Bool_t uturn = kFALSE;
499  Bool_t trackingFailed = kFALSE;
500  Bool_t tooManyStep = kFALSE;
501  while (TMath::Abs(residue) > fgkRungeKuttaMaxResidue && stepNumber <= fgkMaxStepNumber) {
502 
503  dZ = zEnd - trackParam->GetZ();
504  // step lenght assuming linear trajectory
505  step = dZ * TMath::Sqrt(1.0 + trackParam->GetSlopeY()*trackParam->GetSlopeY() +
506  trackParam->GetSlopeX()*trackParam->GetSlopeX());
507  ConvertTrackParamForExtrap(trackParam, forwardBackward, v3);
508 
509  do { // reduce step lenght while zEnd oversteped
510  if (stepNumber > fgkMaxStepNumber) {
511  cout<<"W-AliMFTTrackExtrap::ExtrapToZRungekutta: Too many trials: "<<stepNumber<<endl;
512  tooManyStep = kTRUE;
513  break;
514  }
515  stepNumber ++;
516  step = TMath::Abs(step);
517  if (!AliMFTTrackExtrap::ExtrapOneStepRungekutta(chargeExtrap,step,v3,v3New)) {
518  trackingFailed = kTRUE;
519  break;
520  }
521  residue = zEnd - v3New[2];
522  step *= dZ/(v3New[2]-trackParam->GetZ());
523  } while (residue*dZ < 0 && TMath::Abs(residue) > fgkRungeKuttaMaxResidue);
524 
525  if (trackingFailed) break;
526  else if (v3New[5]*v3[5] < 0) { // the track turned around
527  cout<<"W-AliMFTTrackExtrap::ExtrapToZRungekutta: The track turned around"<<endl;
528  uturn = kTRUE;
529  break;
530  } else RecoverTrackParam(v3New, chargeExtrap * forwardBackward, trackParam);
531 
532  }
533 
534  // terminate the extropolation with a straight line up to the exact "zEnd" value
536 // if (trackingFailed || uturn) {
537 //
538 // // track ends +-100 meters away in the bending direction
539 // dZ = zEnd - v3[2];
540 // Double_t bendingSlope = TMath::Sign(1.e4,-fgSimpleBValue*trackParam->GetInverseBendingMomentum()) / dZ;
541 // Double_t pZ = TMath::Abs(1. / trackParam->GetInverseBendingMomentum()) / TMath::Sqrt(1.0 + bendingSlope * bendingSlope);
542 // Double_t nonBendingSlope = TMath::Sign(TMath::Abs(v3[3]) * v3[6] / pZ, trackParam->GetNonBendingSlope());
543 // trackParam->SetNonBendingCoor(trackParam->GetNonBendingCoor() + dZ * nonBendingSlope);
544 // trackParam->SetNonBendingSlope(nonBendingSlope);
545 // trackParam->SetBendingCoor(trackParam->GetBendingCoor() + dZ * bendingSlope);
546 // trackParam->SetBendingSlope(bendingSlope);
547 // trackParam->SetZ(zEnd);
548 //
549 // return kFALSE;
550 //
551 // } else {
552 //
553 // // track extrapolated normally
554 // trackParam->SetNonBendingCoor(trackParam->GetNonBendingCoor() + residue * trackParam->GetNonBendingSlope());
555 // trackParam->SetBendingCoor(trackParam->GetBendingCoor() + residue * trackParam->GetBendingSlope());
556 // trackParam->SetZ(zEnd);
557 //
558 // return !tooManyStep;
559 //
560 // }
561 
562 }
563 
564 //__________________________________________________________________________
565 void AliMFTTrackExtrap::ConvertTrackParamForExtrap(AliMFTTrackParam* trackParam, Double_t forwardBackward, Double_t *v3)
566 {
570  v3[0] = trackParam->GetX(); // X
571  v3[1] = trackParam->GetY(); // Y
572  v3[2] = trackParam->GetZ(); // Z
573  Double_t slopeX = trackParam->GetSlopeX() ;
574  Double_t slopeY = trackParam->GetSlopeY() ;
575  Double_t slope2 = TMath::Sqrt(1.+slopeX*slopeX +slopeY*slopeY);
576  Double_t pt = TMath::Abs(1.0 / trackParam->GetInverseTransverseMomentum());
577  v3[6] = pt*slope2/TMath::Sqrt(slopeX*slopeX +slopeY*slopeY); // PTOT
578 
579  v3[5] = -forwardBackward / slope2; // PZ/PTOT spectro. z<0
580  v3[3] = slopeX / slope2; // PX/PTOT
581  v3[4] = slopeY / slope2; // PY/PTOT
582 
583 
584 // v3[0] = trackParam->GetX(); // X
585 // v3[1] = trackParam->GetY(); // Y
586 // v3[2] = trackParam->GetZ(); // Z
587 // Double_t pYZ = TMath::Abs(1.0 / trackParam->GetInverseMomentum());
588 // Double_t pZ = pYZ / TMath::Sqrt(1.0 + trackParam->GetSlopeY() * trackParam->GetSlopeY());
589 // v3[6] = TMath::Sqrt(pYZ * pYZ + pZ * pZ * trackParam->GetSlopeX() * trackParam->GetSlopeX()); // PTOT
590 // v3[5] = -forwardBackward * pZ / v3[6]; // PZ/PTOT spectro. z<0
591 // v3[3] = trackParam->GetSlopeX() * v3[5]; // PX/PTOT
592 // v3[4] = trackParam->GetSlopeY() * v3[5]; // PY/PTOT
593 
594 }
595 
596 //__________________________________________________________________________
597 void AliMFTTrackExtrap::RecoverTrackParam(Double_t *v3, Double_t charge, AliMFTTrackParam* trackParam)
598 {
602  trackParam->SetX(v3[0]); // X
603  trackParam->SetY(v3[1]); // Y
604  trackParam->SetZ(v3[2]); // Z
605  Double_t pt = v3[6]*TMath::Sqrt(1. - v3[5]*v3[5]);
606  trackParam->SetInverseTransverseMomentum(charge/pt);
607  trackParam->SetSlopeY(v3[4]/v3[5]);
608  trackParam->SetSlopeX(v3[3]/v3[5]);
609 //
610 // trackParam->SetX(v3[0]); // X
611 // trackParam->SetY(v3[1]); // Y
612 // trackParam->SetZ(v3[2]); // Z
613 // Double_t pYZ = v3[6] * TMath::Sqrt((1.-v3[3])*(1.+v3[3]));
614 // trackParam->SetInverseMomentum(charge/pYZ);
615 // trackParam->SetSlopeY(v3[4]/v3[5]);
616 // trackParam->SetSlopeX(v3[3]/v3[5]);
617 
618 }
619 
620 //__________________________________________________________________________
621 Bool_t AliMFTTrackExtrap::ExtrapToZCov(AliMFTTrackParam* trackParam, Double_t zEnd, Bool_t updatePropagator)
622 {
625  if (trackParam->GetZ() == zEnd) return kTRUE; // nothing to be done if same z
626 
627  if (!fgFieldON) { // linear extrapolation if no magnetic field
628  AliMFTTrackExtrap::LinearExtrapToZCov(trackParam,zEnd,updatePropagator);
629  return kTRUE;
630  }
631 
632  // No need to propagate the covariance matrix if it does not exist
633  if (!trackParam->CovariancesExist()) {
634  cout<<"W-AliMFTTrackExtrap::ExtrapToZCov: Covariance matrix does not exist"<<endl;
635  // Extrapolate track parameters to "zEnd"
636  return ExtrapToZ(trackParam,zEnd);
637  }
638 
639  // Save the actual track parameters
640  AliMFTTrackParam trackParamSave(*trackParam);
641  TMatrixD paramSave(trackParamSave.GetParameters());
642  Double_t zBegin = trackParamSave.GetZ();
643 
644  // Get reference to the parameter covariance matrix
645  const TMatrixD& kParamCov = trackParam->GetCovariances();
646 
647  // Extrapolate track parameters to "zEnd"
648  // Do not update the covariance matrix if the extrapolation failed
649  if (!ExtrapToZ(trackParam,zEnd)) return kFALSE;
650 
651  // Get reference to the extrapolated parameters
652  const TMatrixD& extrapParam = trackParam->GetParameters();
653 
654  // Calculate the jacobian related to the track parameters extrapolation to "zEnd"
655  Bool_t extrapStatus = kTRUE;
656  TMatrixD jacob(5,5);
657  jacob.Zero();
658  TMatrixD dParam(5,1);
659  Double_t direction[5] = {-1.,-1.,1.,1.,-1.};
660  for (Int_t i=0; i<5; i++) {
661  // Skip jacobian calculation for parameters with no associated error
662  if (kParamCov(i,i) <= 0.) continue;
663 
664  // Small variation of parameter i only
665  for (Int_t j=0; j<5; j++) {
666  if (j==i) {
667  dParam(j,0) = TMath::Sqrt(kParamCov(i,i));
668  dParam(j,0) *= TMath::Sign(1.,direction[j]*paramSave(j,0)); // variation always in the same direction
669  } else dParam(j,0) = 0.;
670  }
671  // Set new parameters
672  trackParamSave.SetParameters(paramSave);
673  trackParamSave.AddParameters(dParam);
674  trackParamSave.SetZ(zBegin);
675  // Extrapolate new track parameters to "zEnd"
676  if (!ExtrapToZ(&trackParamSave,zEnd)) {
677  cout<<"W-AliMFTTrackExtrap::ExtrapToZCov: Bad covariance matrix"<<endl;
678  extrapStatus = kFALSE;
679  }
680 
681  // Calculate the jacobian
682  TMatrixD jacobji(trackParamSave.GetParameters(),TMatrixD::kMinus,extrapParam);
683 
684  jacobji *= 1. / dParam(i,0);
685  jacob.SetSub(0,i,jacobji);
686  }
687  cout<<"jacob"<<endl;
688  jacob.Print();
689 
690  // Extrapolate track parameter covariances to "zEnd"
691  cout<<"Initial Cov MAtrix "<<endl;
692  kParamCov.Print();
693  TMatrixD tmp(kParamCov,TMatrixD::kMultTranspose,jacob);
694  TMatrixD tmp2(jacob,TMatrixD::kMult,tmp);
695  cout<<"Extrapolated Cov MAtrix "<<endl;
696  tmp2.Print();
697 
698  trackParam->SetCovariances(tmp2);
699  // Update the propagator if required
700  if (updatePropagator) trackParam->UpdatePropagator(jacob);
701 
702 // return extrapStatus;
703  return kTRUE;
704 }
705 
707 void AliMFTTrackExtrap::AddMCSEffectInAbsorber(AliMFTTrackParam* param, Double_t signedPathLength, Double_t f0, Double_t f1, Double_t f2)
708 {
711 //
712 // // absorber related covariance parameters
713 // Double_t bendingSlope = param->GetSlopeY();
714 // Double_t nonBendingSlope = param->GetSlopeX();
715 // Double_t inverseBendingMomentum = param->GetInverseMomentum();
716 // Double_t alpha2 = 0.0136 * 0.0136 * inverseBendingMomentum * inverseBendingMomentum * (1.0 + bendingSlope * bendingSlope) /
717 // (1.0 + bendingSlope *bendingSlope + nonBendingSlope * nonBendingSlope); // velocity = 1
718 // Double_t pathLength = TMath::Abs(signedPathLength);
719 // Double_t varCoor = alpha2 * (pathLength * pathLength * f0 - 2. * pathLength * f1 + f2);
720 // Double_t covCorrSlope = TMath::Sign(1.,signedPathLength) * alpha2 * (pathLength * f0 - f1);
721 // Double_t varSlop = alpha2 * f0;
722 //
723 // // Set MCS covariance matrix
724 // TMatrixD newParamCov(param->GetCovariances());
725 // // Non bending plane
726 // newParamCov(0,0) += varCoor; newParamCov(0,1) += covCorrSlope;
727 // newParamCov(1,0) += covCorrSlope; newParamCov(1,1) += varSlop;
728 // // Bending plane
729 // newParamCov(2,2) += varCoor; newParamCov(2,3) += covCorrSlope;
730 // newParamCov(3,2) += covCorrSlope; newParamCov(3,3) += varSlop;
731 //
732 // // Set momentum related covariances if B!=0
733 // if (fgFieldON) {
734 // // compute derivative d(q/Pxy) / dSlopeX and d(q/Pxy) / dSlopeY
735 // Double_t dqPxydSlopeX = inverseBendingMomentum * nonBendingSlope / (1. + nonBendingSlope*nonBendingSlope + bendingSlope*bendingSlope);
736 // Double_t dqPxydSlopeY = - inverseBendingMomentum * nonBendingSlope*nonBendingSlope * bendingSlope /
737 // (1. + bendingSlope*bendingSlope) / (1. + nonBendingSlope*nonBendingSlope + bendingSlope*bendingSlope);
738 // // Inverse bending momentum (due to dependences with bending and non bending slopes)
739 // newParamCov(4,0) += dqPxydSlopeX * covCorrSlope; newParamCov(0,4) += dqPxydSlopeX * covCorrSlope;
740 // newParamCov(4,1) += dqPxydSlopeX * varSlop; newParamCov(1,4) += dqPxydSlopeX * varSlop;
741 // newParamCov(4,2) += dqPxydSlopeY * covCorrSlope; newParamCov(2,4) += dqPxydSlopeY * covCorrSlope;
742 // newParamCov(4,3) += dqPxydSlopeY * varSlop; newParamCov(3,4) += dqPxydSlopeY * varSlop;
743 // newParamCov(4,4) += (dqPxydSlopeX*dqPxydSlopeX + dqPxydSlopeY*dqPxydSlopeY) * varSlop;
744 // }
745 //
746 // // Set new covariances
747 // param->SetCovariances(newParamCov);
748 }
749 
750 //__________________________________________________________________________
752  Double_t xVtx, Double_t yVtx, Double_t zVtx,
753  Double_t errXVtx, Double_t errYVtx,
754  Double_t absZBeg, Double_t pathLength, Double_t f0, Double_t f1, Double_t f2)
755 {
756 // /// Correct parameters and corresponding covariances using Branson correction
757 // /// - input param are parameters and covariances at the end of absorber
758 // /// - output param are parameters and covariances at vertex
759 // /// Absorber correction parameters are supposed to be calculated at the current track z-position
760 //
761 // // Position of the Branson plane (spectro. (z<0))
762 // Double_t zB = (f1>0.) ? absZBeg - f2/f1 : 0.;
763 //
764 // // Add MCS effects to current parameter covariances (spectro. (z<0))
765 // AddMCSEffectInAbsorber(param, -pathLength, f0, f1, f2);
766 //
767 // // Get track parameters and covariances in the Branson plane corrected for magnetic field effect
768 // ExtrapToZCov(param,zVtx);
769 // LinearExtrapToZCov(param,zB);
770 //
771 // // compute track parameters at vertex
772 // TMatrixD newParam(5,1);
773 // newParam(0,0) = xVtx;
774 // newParam(1,0) = (param->GetX() - xVtx) / (zB - zVtx);
775 // newParam(2,0) = yVtx;
776 // newParam(3,0) = (param->GetY() - yVtx) / (zB - zVtx);
777 // newParam(4,0) = param->GetCharge() / param->P() *
778 // TMath::Sqrt(1.0 + newParam(1,0)*newParam(1,0) + newParam(3,0)*newParam(3,0)) /
779 // TMath::Sqrt(1.0 + newParam(3,0)*newParam(3,0));
780 //
781 // // Get covariances in (X, SlopeX, Y, SlopeY, q*PTot) coordinate system
782 // TMatrixD paramCovP(param->GetCovariances());
783 // Cov2CovP(param->GetParameters(),paramCovP);
784 //
785 // // Get the covariance matrix in the (XVtx, X, YVtx, Y, q*PTot) coordinate system
786 // TMatrixD paramCovVtx(5,5);
787 // paramCovVtx.Zero();
788 // paramCovVtx(0,0) = errXVtx * errXVtx;
789 // paramCovVtx(1,1) = paramCovP(0,0);
790 // paramCovVtx(2,2) = errYVtx * errYVtx;
791 // paramCovVtx(3,3) = paramCovP(2,2);
792 // paramCovVtx(4,4) = paramCovP(4,4);
793 // paramCovVtx(1,3) = paramCovP(0,2);
794 // paramCovVtx(3,1) = paramCovP(2,0);
795 // paramCovVtx(1,4) = paramCovP(0,4);
796 // paramCovVtx(4,1) = paramCovP(4,0);
797 // paramCovVtx(3,4) = paramCovP(2,4);
798 // paramCovVtx(4,3) = paramCovP(4,2);
799 //
800 // // Jacobian of the transformation (XVtx, X, YVtx, Y, q*PTot) -> (XVtx, SlopeXVtx, YVtx, SlopeYVtx, q*PTotVtx)
801 // TMatrixD jacob(5,5);
802 // jacob.UnitMatrix();
803 // jacob(1,0) = - 1. / (zB - zVtx);
804 // jacob(1,1) = 1. / (zB - zVtx);
805 // jacob(3,2) = - 1. / (zB - zVtx);
806 // jacob(3,3) = 1. / (zB - zVtx);
807 //
808 // // Compute covariances at vertex in the (XVtx, SlopeXVtx, YVtx, SlopeYVtx, q*PTotVtx) coordinate system
809 // TMatrixD tmp(paramCovVtx,TMatrixD::kMultTranspose,jacob);
810 // TMatrixD newParamCov(jacob,TMatrixD::kMult,tmp);
811 //
812 // // Compute covariances at vertex in the (XVtx, SlopeXVtx, YVtx, SlopeYVtx, q/PyzVtx) coordinate system
813 // CovP2Cov(newParam,newParamCov);
814 //
815 // // Set parameters and covariances at vertex
816 // param->SetParameters(newParam);
817 // param->SetZ(zVtx);
818 // param->SetCovariances(newParamCov);
819 }
820 
821 //__________________________________________________________________________
822 void AliMFTTrackExtrap::CorrectELossEffectInAbsorber(AliMFTTrackParam* param, Double_t eLoss, Double_t sigmaELoss2)
823 {
825 
826 // // Get parameter covariances in (X, SlopeX, Y, SlopeY, q*PTot) coordinate system
827 // TMatrixD newParamCov(param->GetCovariances());
828 // Cov2CovP(param->GetParameters(),newParamCov);
829 //
830 // // Compute new parameters corrected for energy loss
831 // Double_t muMass = TDatabasePDG::Instance()->GetParticle("mu-")->Mass(); // GeV
832 // Double_t p = param->P();
833 // Double_t e = TMath::Sqrt(p*p + muMass*muMass);
834 // Double_t eCorr = e + eLoss;
835 // Double_t pCorr = TMath::Sqrt(eCorr*eCorr - muMass*muMass);
836 // Double_t nonBendingSlope = param->GetSlopeX();
837 // Double_t bendingSlope = param->GetSlopeY();
838 // param->SetInverseMomentum(param->GetCharge() / pCorr *
839 // TMath::Sqrt(1.0 + nonBendingSlope*nonBendingSlope + bendingSlope*bendingSlope) /
840 // TMath::Sqrt(1.0 + bendingSlope*bendingSlope));
841 //
842 // // Add effects of energy loss fluctuation to covariances
843 // newParamCov(4,4) += eCorr * eCorr / pCorr / pCorr * sigmaELoss2;
844 //
845 // // Get new parameter covariances in (X, SlopeX, Y, SlopeY, q/Pyz) coordinate system
846 // CovP2Cov(param->GetParameters(),newParamCov);
847 //
848 // // Set new parameter covariances
849 // param->SetCovariances(newParamCov);
850 }
851 
852 //__________________________________________________________________________
853 Bool_t AliMFTTrackExtrap::GetAbsorberCorrectionParam(Double_t trackXYZIn[3], Double_t trackXYZOut[3], Double_t pTotal,
854  Double_t &pathLength, Double_t &f0, Double_t &f1, Double_t &f2,
855  Double_t &meanRho, Double_t &totalELoss, Double_t &sigmaELoss2)
856 {
859  // pathLength: path length between trackXYZIn and trackXYZOut (cm)
860  // f0: 0th moment of z calculated with the inverse radiation-length distribution
861  // f1: 1st moment of z calculated with the inverse radiation-length distribution
862  // f2: 2nd moment of z calculated with the inverse radiation-length distribution
863  // meanRho: average density of crossed material (g/cm3)
864  // totalELoss: total energy loss in absorber
865 
866  // Reset absorber's parameters
867  pathLength = 0.;
868  f0 = 0.;
869  f1 = 0.;
870  f2 = 0.;
871  meanRho = 0.;
872  totalELoss = 0.;
873  sigmaELoss2 = 0.;
874 
875  // Check whether the geometry is available
876  if (!gGeoManager) {
877  cout<<"E-AliMFTTrackExtrap::GetAbsorberCorrectionParam: no TGeo"<<endl;
878  return kFALSE;
879  }
880 
881  // Initialize starting point and direction
882  pathLength = TMath::Sqrt((trackXYZOut[0] - trackXYZIn[0])*(trackXYZOut[0] - trackXYZIn[0])+
883  (trackXYZOut[1] - trackXYZIn[1])*(trackXYZOut[1] - trackXYZIn[1])+
884  (trackXYZOut[2] - trackXYZIn[2])*(trackXYZOut[2] - trackXYZIn[2]));
885  if (pathLength < TGeoShape::Tolerance()) return kFALSE;
886  Double_t b[3];
887  b[0] = (trackXYZOut[0] - trackXYZIn[0]) / pathLength;
888  b[1] = (trackXYZOut[1] - trackXYZIn[1]) / pathLength;
889  b[2] = (trackXYZOut[2] - trackXYZIn[2]) / pathLength;
890  TGeoNode *currentnode = gGeoManager->InitTrack(trackXYZIn, b);
891  if (!currentnode) {
892  cout<<"E-AliMFTTrackExtrap::GetAbsorberCorrectionParam: start point out of geometry"<<endl;
893  return kFALSE;
894  }
895 
896  // loop over absorber slices and calculate absorber's parameters
897  Double_t rho = 0.; // material density (g/cm3)
898  Double_t x0 = 0.; // radiation-length (cm-1)
899  Double_t atomicA = 0.; // A of material
900  Double_t atomicZ = 0.; // Z of material
901  Double_t atomicZoverA = 0.; // Z/A of material
902  Double_t localPathLength = 0;
903  Double_t remainingPathLength = pathLength;
904  Double_t sigmaELoss = 0.;
905  Double_t zB = trackXYZIn[2];
906  Double_t zE, dzB, dzE;
907  do {
908  // Get material properties
909  TGeoMaterial *material = currentnode->GetVolume()->GetMedium()->GetMaterial();
910  rho = material->GetDensity();
911  x0 = material->GetRadLen();
912  atomicA = material->GetA();
913  atomicZ = material->GetZ();
914  if(material->IsMixture()){
915  TGeoMixture * mixture = (TGeoMixture*)material;
916  atomicZoverA = 0.;
917  Double_t sum = 0.;
918  for (Int_t iel=0;iel<mixture->GetNelements();iel++){
919  sum += mixture->GetWmixt()[iel];
920  atomicZoverA += mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel];
921  }
922  atomicZoverA/=sum;
923  }
924  else atomicZoverA = atomicZ/atomicA;
925 
926  // Get path length within this material
927  gGeoManager->FindNextBoundary(remainingPathLength);
928  localPathLength = gGeoManager->GetStep() + 1.e-6;
929  // Check if boundary within remaining path length. If so, make sure to cross the boundary to prepare the next step
930  if (localPathLength >= remainingPathLength) localPathLength = remainingPathLength;
931  else {
932  currentnode = gGeoManager->Step();
933  if (!currentnode) {
934  cout<<"E-AliMFTTrackExtrap::GetAbsorberCorrectionParam: navigation failed"<<endl;
935  f0 = f1 = f2 = meanRho = totalELoss = sigmaELoss2 = 0.;
936  return kFALSE;
937  }
938  if (!gGeoManager->IsEntering()) {
939  // make another small step to try to enter in new absorber slice
940  gGeoManager->SetStep(0.001);
941  currentnode = gGeoManager->Step();
942  if (!gGeoManager->IsEntering() || !currentnode) {
943  cout<<"E-AliMFTTrackExtrap::GetAbsorberCorrectionParam: navigation failed"<<endl;
944  f0 = f1 = f2 = meanRho = totalELoss = sigmaELoss2 = 0.;
945  return kFALSE;
946  }
947  localPathLength += 0.001;
948  }
949  }
950 
951  // calculate absorber's parameters
952  zE = b[2] * localPathLength + zB;
953  dzB = zB - trackXYZIn[2];
954  dzE = zE - trackXYZIn[2];
955  f0 += localPathLength / x0;
956  f1 += (dzE*dzE - dzB*dzB) / b[2] / b[2] / x0 / 2.;
957  f2 += (dzE*dzE*dzE - dzB*dzB*dzB) / b[2] / b[2] / b[2] / x0 / 3.;
958  meanRho += localPathLength * rho;
959  totalELoss += BetheBloch(pTotal, localPathLength, rho, atomicZ, atomicZoverA);
960  sigmaELoss += EnergyLossFluctuation(pTotal, localPathLength, rho, atomicZoverA);
961 
962  // prepare next step
963  zB = zE;
964  remainingPathLength -= localPathLength;
965  } while (remainingPathLength > TGeoShape::Tolerance());
966 
967  meanRho /= pathLength;
968  sigmaELoss2 = sigmaELoss*sigmaELoss;
969 
970  return kTRUE;
971 }
972 
973 //__________________________________________________________________________
974 Double_t AliMFTTrackExtrap::GetMCSAngle2(const AliMFTTrackParam& param, Double_t dZ, Double_t x0)
975 {
979  return 0.;
980 // Double_t bendingSlope = param.GetSlopeY();
981 // Double_t nonBendingSlope = param.GetSlopeX();
982 // Double_t inverseTotalMomentum2 = param.GetInverseMomentum() * param.GetInverseMomentum() *
983 // (1.0 + bendingSlope * bendingSlope) /
984 // (1.0 + bendingSlope *bendingSlope + nonBendingSlope * nonBendingSlope);
985 // // Path length in the material
986 // Double_t pathLength = TMath::Abs(dZ) * TMath::Sqrt(1.0 + bendingSlope*bendingSlope + nonBendingSlope*nonBendingSlope);
987 // // relativistic velocity
988 // Double_t velo = 1.;
989 // // Angular dispersion square of the track (variance) in a plane perpendicular to the trajectory
990 // Double_t theta02 = 0.0136 / velo * (1 + 0.038 * TMath::Log(pathLength/x0));
991 //
992 // return theta02 * theta02 * inverseTotalMomentum2 * pathLength / x0;
993 }
994 
995 
996 
997 //__________________________________________________________________________
998 void AliMFTTrackExtrap::AddMCSEffect(AliMFTTrackParam *param, Double_t dZ, Double_t x0)
999 {
1005  Double_t slopeX = param->GetSlopeX();
1006  Double_t slopeY = param->GetSlopeY();
1007  Double_t slope2 = slopeX*slopeX+slopeY*slopeY;
1008 
1009  Double_t inversePt = TMath::Abs(param->GetInverseTransverseMomentum());
1010 
1011  Double_t inverseTotalMomentum2 = inversePt*inversePt / (1. + 1./slope2 );
1012  // Path length in the material
1013  Double_t signedPathLength = dZ * TMath::Sqrt(1.0 + slope2);
1014  Double_t pathLengthOverX0 = (x0 > 0.) ? TMath::Abs(signedPathLength * x0 /dZ) : TMath::Abs(signedPathLength);
1015  // relativistic velocity
1016  Double_t velo = 1.;
1017  // Angular dispersion square of the track (variance) in a plane perpendicular to the trajectory
1018  Double_t theta02 = 0.0136 / velo * (1 + 0.038 * TMath::Log(pathLengthOverX0));
1019  theta02 *= theta02 * inverseTotalMomentum2 * pathLengthOverX0;
1020 
1021  Double_t varCoor = (x0 > 0.) ? signedPathLength * signedPathLength * theta02 / 3. : 0.;
1022  Double_t varSlop = theta02;
1023  Double_t covCorrSlope = (x0 > 0.) ? signedPathLength * theta02/ 2. : 0.;
1024 
1025  cout<<Form("theta02=%e inverseTotalMomentum2=%e signedPathLength=%e pathLengthOverX0=%e ",theta02,inverseTotalMomentum2,signedPathLength,pathLengthOverX0 )<<endl;
1026 
1027  cout<<Form("dz=%e x0=%e varCoor=%e varSlop=%e covCorrSlope=%e ",dZ,x0,varCoor,varSlop,covCorrSlope )<<endl;
1028  // Set MCS covariance matrix
1029  TMatrixD newParamCov(param->GetCovariances());
1030  cout<<"Covariance avant MCS"<<endl;
1031  newParamCov.Print();
1032  // Non bending plane
1033  newParamCov(0,0) += varCoor; newParamCov(0,2) += covCorrSlope;
1034  newParamCov(2,0) += covCorrSlope; newParamCov(2,2) += varSlop;
1035  // Bending plane
1036  newParamCov(1,1) += varCoor; newParamCov(1,3) += covCorrSlope;
1037  newParamCov(3,1) += covCorrSlope; newParamCov(3,3) += varSlop;
1038 
1039  cout<<"Covariance apres MCS"<<endl;
1040  newParamCov.Print();
1041 
1042 // // Set momentum related covariances if B!=0
1043 // if (fgFieldON) {
1044 // // compute derivative d(q/Pxy) / dSlopeX and d(q/Pxy) / dSlopeY
1045 // Double_t dqPxydSlopeX = inverseBendingMomentum * nonBendingSlope / (1. + nonBendingSlope*nonBendingSlope + bendingSlope*bendingSlope);
1046 // Double_t dqPxydSlopeY = - inverseBendingMomentum * nonBendingSlope*nonBendingSlope * bendingSlope /
1047 // (1. + bendingSlope*bendingSlope) / (1. + nonBendingSlope*nonBendingSlope + bendingSlope*bendingSlope);
1048 // // Inverse bending momentum (due to dependences with bending and non bending slopes)
1049 // newParamCov(4,0) += dqPxydSlopeX * covCorrSlope; newParamCov(0,4) += dqPxydSlopeX * covCorrSlope;
1050 // newParamCov(4,1) += dqPxydSlopeX * varSlop; newParamCov(1,4) += dqPxydSlopeX * varSlop;
1051 // newParamCov(4,2) += dqPxydSlopeY * covCorrSlope; newParamCov(2,4) += dqPxydSlopeY * covCorrSlope;
1052 // newParamCov(4,3) += dqPxydSlopeY * varSlop; newParamCov(3,4) += dqPxydSlopeY * varSlop;
1053 // newParamCov(4,4) += (dqPxydSlopeX*dqPxydSlopeX + dqPxydSlopeY*dqPxydSlopeY) * varSlop;
1054 // }
1055 // cout<<"Covariance apres"<<endl;
1056 // newParamCov.Print();
1057 
1058  // Set new covariances
1059  param->SetCovariances(newParamCov);
1060 }
1061 
1062 //__________________________________________________________________________
1064  Double_t xVtx, Double_t yVtx, Double_t zVtx,
1065  Double_t errXVtx, Double_t errYVtx,
1066  Bool_t correctForMCS, Bool_t correctForEnergyLoss)
1067 {
1075 
1076 // if (trackParam->GetZ() == zVtx) return; // nothing to be done if already at vertex
1077 //
1078 // if (trackParam->GetZ() > zVtx) { // spectro. (z<0)
1079 // cout<<"E-AliMFTTrackExtrap::ExtrapToVertex: Starting Z ("<<trackParam->GetZ()
1080 // <<") upstream the vertex (zVtx = "<<zVtx<<")"<<endl;
1081 // return;
1082 // }
1083 //
1084 // // Check the vertex position relatively to the absorber
1085 // if (zVtx < AliMUONConstants::AbsZBeg() && zVtx > AliMUONConstants::AbsZEnd()) { // spectro. (z<0)
1086 // cout<<"W-AliMFTTrackExtrap::ExtrapToVertex: Ending Z ("<<zVtx
1087 // <<") inside the front absorber ("<<AliMUONConstants::AbsZBeg()<<","<<AliMUONConstants::AbsZEnd()<<")"<<endl;
1088 // } else if (zVtx < AliMUONConstants::AbsZEnd() ) { // spectro. (z<0)
1089 // cout<<"W-AliMFTTrackExtrap::ExtrapToVertex: Ending Z ("<<zVtx
1090 // <<") downstream the front absorber (zAbsorberEnd = "<<AliMUONConstants::AbsZEnd()<<")"<<endl;
1091 // if (trackParam->CovariancesExist()) ExtrapToZCov(trackParam,zVtx);
1092 // else ExtrapToZ(trackParam,zVtx);
1093 // return;
1094 // }
1095 //
1096 // // Check the track position relatively to the absorber and extrapolate track parameters to the end of the absorber if needed
1097 // if (trackParam->GetZ() > AliMUONConstants::AbsZBeg()) { // spectro. (z<0)
1098 // cout<<"W-AliMFTTrackExtrap::ExtrapToVertex: Starting Z ("<<trackParam->GetZ()
1099 // <<") upstream the front absorber (zAbsorberBegin = "<<AliMUONConstants::AbsZBeg()<<")"<<endl;
1100 // if (trackParam->CovariancesExist()) ExtrapToZCov(trackParam,zVtx);
1101 // else ExtrapToZ(trackParam,zVtx);
1102 // return;
1103 // } else if (trackParam->GetZ() > AliMUONConstants::AbsZEnd()) { // spectro. (z<0)
1104 // cout<<"W-AliMFTTrackExtrap::ExtrapToVertex: Starting Z ("<<trackParam->GetZ()
1105 // <<") inside the front absorber ("<<AliMUONConstants::AbsZBeg()<<","<<AliMUONConstants::AbsZEnd()<<")"<<endl;
1106 // } else {
1107 // if (trackParam->CovariancesExist()) ExtrapToZCov(trackParam,AliMUONConstants::AbsZEnd());
1108 // else ExtrapToZ(trackParam,AliMUONConstants::AbsZEnd());
1109 // }
1110 //
1111 // // Get absorber correction parameters assuming linear propagation in absorber
1112 // Double_t trackXYZOut[3];
1113 // trackXYZOut[0] = trackParam->GetX();
1114 // trackXYZOut[1] = trackParam->GetY();
1115 // trackXYZOut[2] = trackParam->GetZ();
1116 // Double_t trackXYZIn[3];
1117 // if (correctForMCS) { // assume linear propagation until the vertex
1118 // trackXYZIn[2] = TMath::Min(zVtx, AliMUONConstants::AbsZBeg()); // spectro. (z<0)
1119 // trackXYZIn[0] = trackXYZOut[0] + (xVtx - trackXYZOut[0]) / (zVtx - trackXYZOut[2]) * (trackXYZIn[2] - trackXYZOut[2]);
1120 // trackXYZIn[1] = trackXYZOut[1] + (yVtx - trackXYZOut[1]) / (zVtx - trackXYZOut[2]) * (trackXYZIn[2] - trackXYZOut[2]);
1121 // } else {
1122 // AliMFTTrackParam trackParamIn(*trackParam);
1123 // ExtrapToZ(&trackParamIn, TMath::Min(zVtx, AliMUONConstants::AbsZBeg()));
1124 // trackXYZIn[0] = trackParamIn.GetX();
1125 // trackXYZIn[1] = trackParamIn.GetY();
1126 // trackXYZIn[2] = trackParamIn.GetZ();
1127 // }
1128 // Double_t pTot = trackParam->P();
1129 // Double_t pathLength, f0, f1, f2, meanRho, totalELoss, sigmaELoss2;
1130 // if (!GetAbsorberCorrectionParam(trackXYZIn,trackXYZOut,pTot,pathLength,f0,f1,f2,meanRho,totalELoss,sigmaELoss2)) {
1131 // cout<<"E-AliMFTTrackExtrap::ExtrapToVertex: Unable to take into account the absorber effects"<<endl;
1132 // if (trackParam->CovariancesExist()) ExtrapToZCov(trackParam,zVtx);
1133 // else ExtrapToZ(trackParam,zVtx);
1134 // return;
1135 // }
1136 //
1137 // // Compute track parameters and covariances at vertex according to correctForMCS and correctForEnergyLoss flags
1138 // if (correctForMCS) {
1139 //
1140 // if (correctForEnergyLoss) {
1141 //
1142 // // Correct for multiple scattering and energy loss
1143 // CorrectELossEffectInAbsorber(trackParam, 0.5*totalELoss, 0.5*sigmaELoss2);
1144 // CorrectMCSEffectInAbsorber(trackParam, xVtx, yVtx, zVtx, errXVtx, errYVtx,
1145 // trackXYZIn[2], pathLength, f0, f1, f2);
1146 // CorrectELossEffectInAbsorber(trackParam, 0.5*totalELoss, 0.5*sigmaELoss2);
1147 //
1148 // } else {
1149 //
1150 // // Correct for multiple scattering
1151 // CorrectMCSEffectInAbsorber(trackParam, xVtx, yVtx, zVtx, errXVtx, errYVtx,
1152 // trackXYZIn[2], pathLength, f0, f1, f2);
1153 // }
1154 //
1155 // } else {
1156 //
1157 // if (correctForEnergyLoss) {
1158 //
1159 // // Correct for energy loss add multiple scattering dispersion in covariance matrix
1160 // CorrectELossEffectInAbsorber(trackParam, 0.5*totalELoss, 0.5*sigmaELoss2);
1161 // AddMCSEffectInAbsorber(trackParam, -pathLength, f0, f1, f2); // (spectro. (z<0))
1162 // ExtrapToZCov(trackParam, trackXYZIn[2]);
1163 // CorrectELossEffectInAbsorber(trackParam, 0.5*totalELoss, 0.5*sigmaELoss2);
1164 // ExtrapToZCov(trackParam, zVtx);
1165 //
1166 // } else {
1167 //
1168 // // add multiple scattering dispersion in covariance matrix
1169 // AddMCSEffectInAbsorber(trackParam, -pathLength, f0, f1, f2); // (spectro. (z<0))
1170 // ExtrapToZCov(trackParam, zVtx);
1171 //
1172 // }
1173 //
1174 // }
1175 
1176 }
1177 
1178 //__________________________________________________________________________
1180  Double_t xVtx, Double_t yVtx, Double_t zVtx,
1181  Double_t errXVtx, Double_t errYVtx)
1182 {
1185  ExtrapToVertex(trackParam, xVtx, yVtx, zVtx, errXVtx, errYVtx, kTRUE, kTRUE);
1186 }
1187 
1188 //__________________________________________________________________________
1190  Double_t xVtx, Double_t yVtx, Double_t zVtx,
1191  Double_t errXVtx, Double_t errYVtx)
1192 {
1195  ExtrapToVertex(trackParam, xVtx, yVtx, zVtx, errXVtx, errYVtx, kTRUE, kFALSE);
1196 }
1197 
1198 //__________________________________________________________________________
1200 {
1203  ExtrapToVertex(trackParam, 0., 0., zVtx, 0., 0., kFALSE, kTRUE);
1204 }
1205 
1206 //__________________________________________________________________________
1208 {
1211  ExtrapToVertex(trackParam, 0., 0., zVtx, 0., 0., kFALSE, kFALSE);
1212 }
1213 
1214 //__________________________________________________________________________
1215 Double_t AliMFTTrackExtrap::TotalMomentumEnergyLoss(AliMFTTrackParam* trackParam, Double_t xVtx, Double_t yVtx, Double_t zVtx)
1216 {
1218 
1219  if (trackParam->GetZ() == zVtx) return 0.; // nothing to be done if already at vertex
1220 
1221  // Check whether the geometry is available
1222  if (!gGeoManager) {
1223  cout<<"E-AliMFTTrackExtrap::TotalMomentumEnergyLoss: no TGeo"<<endl;
1224  return 0.;
1225  }
1226 
1227  // Get encountered material correction parameters assuming linear propagation from vertex to the track position
1228  Double_t trackXYZOut[3];
1229  trackXYZOut[0] = trackParam->GetX();
1230  trackXYZOut[1] = trackParam->GetY();
1231  trackXYZOut[2] = trackParam->GetZ();
1232  Double_t trackXYZIn[3];
1233  trackXYZIn[0] = xVtx;
1234  trackXYZIn[1] = yVtx;
1235  trackXYZIn[2] = zVtx;
1236  Double_t pTot = trackParam->P();
1237  Double_t pathLength, f0, f1, f2, meanRho, totalELoss, sigmaELoss2;
1238  GetAbsorberCorrectionParam(trackXYZIn,trackXYZOut,pTot,pathLength,f0,f1,f2,meanRho,totalELoss,sigmaELoss2);
1239 
1240  // total momentum corrected for energy loss
1241  Double_t muMass = TDatabasePDG::Instance()->GetParticle("mu-")->Mass(); // GeV
1242  Double_t e = TMath::Sqrt(pTot*pTot + muMass*muMass);
1243  Double_t eCorr = e + totalELoss;
1244  Double_t pTotCorr = TMath::Sqrt(eCorr*eCorr - muMass*muMass);
1245 
1246  return pTotCorr - pTot;
1247 }
1248 
1249 //__________________________________________________________________________
1250 Double_t AliMFTTrackExtrap::BetheBloch(Double_t pTotal, Double_t pathLength, Double_t rho, Double_t atomicZ, Double_t atomicZoverA)
1251 {
1254  Double_t muMass = TDatabasePDG::Instance()->GetParticle("mu-")->Mass(); // GeV
1255 
1256  // mean exitation energy (GeV)
1257  Double_t i;
1258  if (atomicZ < 13) i = (12. * atomicZ + 7.) * 1.e-9;
1259  else i = (9.76 * atomicZ + 58.8 * TMath::Power(atomicZ,-0.19)) * 1.e-9;
1260 
1261  return pathLength * rho * AliExternalTrackParam::BetheBlochGeant(pTotal/muMass, rho, 0.20, 3.00, i, atomicZoverA);
1262 }
1263 
1264 //__________________________________________________________________________
1265 Double_t AliMFTTrackExtrap::EnergyLossFluctuation(Double_t pTotal, Double_t pathLength, Double_t rho, Double_t atomicZoverA)
1266 {
1269  Double_t muMass = TDatabasePDG::Instance()->GetParticle("mu-")->Mass(); // GeV
1270  //Double_t eMass = 0.510998918e-3; // GeV
1271  Double_t k = 0.307075e-3; // GeV.g^-1.cm^2
1272  Double_t p2=pTotal*pTotal;
1273  Double_t beta2=p2/(p2 + muMass*muMass);
1274 
1275  Double_t fwhm = 2. * k * rho * pathLength * atomicZoverA / beta2; // FWHM of the energy loss Landau distribution
1276  Double_t sigma = fwhm / TMath::Sqrt(8.*log(2.)); // gaussian: fwmh = 2 * srqt(2*ln(2)) * sigma (i.e. fwmh = 2.35 * sigma)
1277 
1278  //sigma2 = k * rho * pathLength * atomicZ / atomicA * eMass; // sigma2 of the energy loss gaussian distribution
1279 
1280  return sigma;
1281 }
1282 
1283 //__________________________________________________________________________
1285 {
1288 
1289  // charge * total momentum
1290  Double_t qPTot = TMath::Sqrt(1. + param(1,0)*param(1,0) + param(3,0)*param(3,0)) /
1291  TMath::Sqrt(1. + param(3,0)*param(3,0)) / param(4,0);
1292 
1293  // Jacobian of the opposite transformation
1294  TMatrixD jacob(5,5);
1295  jacob.UnitMatrix();
1296  jacob(4,1) = qPTot * param(1,0) / (1. + param(1,0)*param(1,0) + param(3,0)*param(3,0));
1297  jacob(4,3) = - qPTot * param(1,0) * param(1,0) * param(3,0) /
1298  (1. + param(3,0)*param(3,0)) / (1. + param(1,0)*param(1,0) + param(3,0)*param(3,0));
1299  jacob(4,4) = - qPTot / param(4,0);
1300 
1301  // compute covariances in new coordinate system
1302  TMatrixD tmp(cov,TMatrixD::kMultTranspose,jacob);
1303  cov.Mult(jacob,tmp);
1304 }
1305 
1306 //__________________________________________________________________________
1308 {
1311 
1312  // charge * total momentum
1313  Double_t qPTot = TMath::Sqrt(1. + param(1,0)*param(1,0) + param(3,0)*param(3,0)) /
1314  TMath::Sqrt(1. + param(3,0)*param(3,0)) / param(4,0);
1315 
1316  // Jacobian of the transformation
1317  TMatrixD jacob(5,5);
1318  jacob.UnitMatrix();
1319  jacob(4,1) = param(4,0) * param(1,0) / (1. + param(1,0)*param(1,0) + param(3,0)*param(3,0));
1320  jacob(4,3) = - param(4,0) * param(1,0) * param(1,0) * param(3,0) /
1321  (1. + param(3,0)*param(3,0)) / (1. + param(1,0)*param(1,0) + param(3,0)*param(3,0));
1322  jacob(4,4) = - param(4,0) / qPTot;
1323 
1324  // compute covariances in new coordinate system
1325  TMatrixD tmp(covP,TMatrixD::kMultTranspose,jacob);
1326  covP.Mult(jacob,tmp);
1327 }
1328 
1329  //__________________________________________________________________________
1330 void AliMFTTrackExtrap::ExtrapOneStepHelix(Double_t charge, Double_t step, const Double_t *vect, Double_t *vout)
1331 {
1353 
1354 // modif: everything in double precision
1355 
1356  Double_t xyz[3], h[4], hxp[3];
1357  Double_t h2xy, hp, rho, tet;
1358  Double_t sint, sintt, tsint, cos1t;
1359  Double_t f1, f2, f3, f4, f5, f6;
1360 
1361  const Int_t kix = 0;
1362  const Int_t kiy = 1;
1363  const Int_t kiz = 2;
1364  const Int_t kipx = 3;
1365  const Int_t kipy = 4;
1366  const Int_t kipz = 5;
1367  const Int_t kipp = 6;
1368 // cout<<"vin ="<< vect[kix]<<" "<< vect[kiy]<<" "<< vect[kiz]<<" pxyz/ptot "<< vect[kipx]<<" "<< vect[kipy]<<" "<< vect[kipz]<<endl;
1369 
1370  const Double_t kec = 2.9979251e-4;
1371  //
1372  // ------------------------------------------------------------------
1373  //
1374  // units are kgauss,centimeters,gev/c
1375  //
1376  vout[kipp] = vect[kipp];
1377  if (TMath::Abs(charge) < 0.00001) {
1378  for (Int_t i = 0; i < 3; i++) {
1379  vout[i] = vect[i] + step * vect[i+3];
1380  vout[i+3] = vect[i+3];
1381  }
1382  return;
1383  }
1384  xyz[0] = vect[kix] + 0.5 * step * vect[kipx];
1385  xyz[1] = vect[kiy] + 0.5 * step * vect[kipy];
1386  xyz[2] = vect[kiz] + 0.5 * step * vect[kipz];
1387 
1388  //cmodif: call gufld (xyz, h) changed into:
1389  TGeoGlobalMagField::Instance()->Field(xyz,h);
1390  h2xy = h[0]*h[0] + h[1]*h[1];
1391  h[3] = h[2]*h[2]+ h2xy;
1392 // cout<<"Field ="<< h[0]<<" "<< h[1]<<" "<< h[2]<<" "<< h[3]<<endl;
1393  if (h[3] < 1.e-12) {
1394  for (Int_t i = 0; i < 3; i++) {
1395  vout[i] = vect[i] + step * vect[i+3];
1396  vout[i+3] = vect[i+3];
1397  }
1398  return;
1399  }
1400  if (h2xy < 1.e-12*h[3]) {
1401  ExtrapOneStepHelix3(charge*h[2], step, vect, vout);
1402  return;
1403  }
1404  h[3] = TMath::Sqrt(h[3]);
1405  h[0] /= h[3];
1406  h[1] /= h[3];
1407  h[2] /= h[3];
1408  h[3] *= kec;
1409 
1410  hxp[0] = h[1]*vect[kipz] - h[2]*vect[kipy];
1411  hxp[1] = h[2]*vect[kipx] - h[0]*vect[kipz];
1412  hxp[2] = h[0]*vect[kipy] - h[1]*vect[kipx];
1413 
1414  hp = h[0]*vect[kipx] + h[1]*vect[kipy] + h[2]*vect[kipz];
1415 
1416  rho = -charge*h[3]/vect[kipp];
1417  tet = rho * step;
1418 
1419  if (TMath::Abs(tet) > 0.15) {
1420  sint = TMath::Sin(tet);
1421  sintt = (sint/tet);
1422  tsint = (tet-sint)/tet;
1423  cos1t = 2.*(TMath::Sin(0.5*tet))*(TMath::Sin(0.5*tet))/tet;
1424  } else {
1425  tsint = tet*tet/36.;
1426  sintt = (1. - tsint);
1427  sint = tet*sintt;
1428  cos1t = 0.5*tet;
1429  }
1430 
1431  f1 = step * sintt;
1432  f2 = step * cos1t;
1433  f3 = step * tsint * hp;
1434  f4 = -tet*cos1t;
1435  f5 = sint;
1436  f6 = tet * cos1t * hp;
1437 
1438  vout[kix] = vect[kix] + f1*vect[kipx] + f2*hxp[0] + f3*h[0];
1439  vout[kiy] = vect[kiy] + f1*vect[kipy] + f2*hxp[1] + f3*h[1];
1440  vout[kiz] = vect[kiz] + f1*vect[kipz] + f2*hxp[2] + f3*h[2];
1441 
1442  vout[kipx] = vect[kipx] + f4*vect[kipx] + f5*hxp[0] + f6*h[0];
1443  vout[kipy] = vect[kipy] + f4*vect[kipy] + f5*hxp[1] + f6*h[1];
1444  vout[kipz] = vect[kipz] + f4*vect[kipz] + f5*hxp[2] + f6*h[2];
1445 // cout<<"vout ="<< vout[kix]<<" "<< vout[kiy]<<" "<< vout[kiz]<<" pxyz/ptot "<< vout[kipx]<<" "<< vout[kipy]<<" "<< vout[kipz]<<endl;
1446 // cout<<"vlin ="<< vect[kix] + step * vect[kix+3]<<" "<< vect[kiy] + step * vect[kiy+3]<<" "<< vect[kiz] + step * vect[kiz+3]<<" pxyz/ptot "<< vect[kipx]<<" "<< vect[kipy]<<" "<< vect[kipz]<<endl;
1447 
1448 
1449  return;
1450 }
1451 
1452  //__________________________________________________________________________
1453 void AliMFTTrackExtrap::ExtrapOneStepHelix3(Double_t field, Double_t step, const Double_t *vect, Double_t *vout)
1454 {
1469 
1470  Double_t hxp[3];
1471  Double_t h4, hp, rho, tet;
1472  Double_t sint, sintt, tsint, cos1t;
1473  Double_t f1, f2, f3, f4, f5, f6;
1474 
1475  const Int_t kix = 0;
1476  const Int_t kiy = 1;
1477  const Int_t kiz = 2;
1478  const Int_t kipx = 3;
1479  const Int_t kipy = 4;
1480  const Int_t kipz = 5;
1481  const Int_t kipp = 6;
1482 
1483  const Double_t kec = 2.9979251e-4;
1484 
1485 //
1486 // ------------------------------------------------------------------
1487 //
1488 // units are kgauss,centimeters,gev/c
1489 //
1490  vout[kipp] = vect[kipp];
1491  h4 = field * kec;
1492 
1493  hxp[0] = - vect[kipy];
1494  hxp[1] = + vect[kipx];
1495 
1496  hp = vect[kipz];
1497 
1498  rho = -h4/vect[kipp];
1499  tet = rho * step;
1500  if (TMath::Abs(tet) > 0.15) {
1501  sint = TMath::Sin(tet);
1502  sintt = (sint/tet);
1503  tsint = (tet-sint)/tet;
1504  cos1t = 2.* TMath::Sin(0.5*tet) * TMath::Sin(0.5*tet)/tet;
1505  } else {
1506  tsint = tet*tet/36.;
1507  sintt = (1. - tsint);
1508  sint = tet*sintt;
1509  cos1t = 0.5*tet;
1510  }
1511 
1512  f1 = step * sintt;
1513  f2 = step * cos1t;
1514  f3 = step * tsint * hp;
1515  f4 = -tet*cos1t;
1516  f5 = sint;
1517  f6 = tet * cos1t * hp;
1518 
1519  vout[kix] = vect[kix] + f1*vect[kipx] + f2*hxp[0];
1520  vout[kiy] = vect[kiy] + f1*vect[kipy] + f2*hxp[1];
1521  vout[kiz] = vect[kiz] + f1*vect[kipz] + f3;
1522 
1523  vout[kipx] = vect[kipx] + f4*vect[kipx] + f5*hxp[0];
1524  vout[kipy] = vect[kipy] + f4*vect[kipy] + f5*hxp[1];
1525  vout[kipz] = vect[kipz] + f4*vect[kipz] + f6;
1526 
1527  return;
1528 }
1529 
1530  //__________________________________________________________________________
1531 Bool_t AliMFTTrackExtrap::ExtrapOneStepRungekutta(Double_t charge, Double_t step, const Double_t* vect, Double_t* vout)
1532 {
1556 
1557  Double_t h2, h4, f[4];
1558  Double_t xyzt[3] = {FLT_MAX, FLT_MAX, FLT_MAX};
1559  Double_t a, b, c, ph,ph2;
1560  Double_t secxs[4],secys[4],seczs[4],hxp[3];
1561  Double_t g1, g2, g3, g4, g5, g6, ang2, dxt, dyt, dzt;
1562  Double_t est, at, bt, ct, cba;
1563  Double_t f1, f2, f3, f4, rho, tet, hnorm, hp, rho1, sint, cost;
1564 
1565  Double_t x;
1566  Double_t y;
1567  Double_t z;
1568 
1569  Double_t xt;
1570  Double_t yt;
1571  Double_t zt;
1572 
1573  Double_t maxit = 1992;
1574  Double_t maxcut = 11;
1575 
1576  const Double_t kdlt = 1e-4;
1577  const Double_t kdlt32 = kdlt/32.;
1578  const Double_t kthird = 1./3.;
1579  const Double_t khalf = 0.5;
1580  const Double_t kec = 2.9979251e-4;
1581 
1582  const Double_t kpisqua = 9.86960440109;
1583  const Int_t kix = 0;
1584  const Int_t kiy = 1;
1585  const Int_t kiz = 2;
1586  const Int_t kipx = 3;
1587  const Int_t kipy = 4;
1588  const Int_t kipz = 5;
1589 
1590  // *.
1591  // *. ------------------------------------------------------------------
1592  // *.
1593  // * this constant is for units cm,gev/c and kgauss
1594  // *
1595  Int_t iter = 0;
1596  Int_t ncut = 0;
1597  for(Int_t j = 0; j < 7; j++)
1598  vout[j] = vect[j];
1599 
1600  Double_t pinv = kec * charge / vect[6];
1601  Double_t tl = 0.;
1602  Double_t h = step;
1603  Double_t rest;
1604 
1605 
1606  do {
1607  rest = step - tl;
1608  if (TMath::Abs(h) > TMath::Abs(rest)) h = rest;
1609  //cmodif: call gufld(vout,f) changed into:
1610  TGeoGlobalMagField::Instance()->Field(vout,f);
1611 
1612  // *
1613  // * start of integration
1614  // *
1615  x = vout[0];
1616  y = vout[1];
1617  z = vout[2];
1618  a = vout[3];
1619  b = vout[4];
1620  c = vout[5];
1621 
1622  h2 = khalf * h;
1623  h4 = khalf * h2;
1624  ph = pinv * h;
1625  ph2 = khalf * ph;
1626  secxs[0] = (b * f[2] - c * f[1]) * ph2;
1627  secys[0] = (c * f[0] - a * f[2]) * ph2;
1628  seczs[0] = (a * f[1] - b * f[0]) * ph2;
1629  ang2 = (secxs[0]*secxs[0] + secys[0]*secys[0] + seczs[0]*seczs[0]);
1630  if (ang2 > kpisqua) break;
1631 
1632  dxt = h2 * a + h4 * secxs[0];
1633  dyt = h2 * b + h4 * secys[0];
1634  dzt = h2 * c + h4 * seczs[0];
1635  xt = x + dxt;
1636  yt = y + dyt;
1637  zt = z + dzt;
1638  // *
1639  // * second intermediate point
1640  // *
1641 
1642  est = TMath::Abs(dxt) + TMath::Abs(dyt) + TMath::Abs(dzt);
1643  if (est > h) {
1644  if (ncut++ > maxcut) break;
1645  h *= khalf;
1646  continue;
1647  }
1648 
1649  xyzt[0] = xt;
1650  xyzt[1] = yt;
1651  xyzt[2] = zt;
1652 
1653  //cmodif: call gufld(xyzt,f) changed into:
1654  TGeoGlobalMagField::Instance()->Field(xyzt,f);
1655 
1656  at = a + secxs[0];
1657  bt = b + secys[0];
1658  ct = c + seczs[0];
1659 
1660  secxs[1] = (bt * f[2] - ct * f[1]) * ph2;
1661  secys[1] = (ct * f[0] - at * f[2]) * ph2;
1662  seczs[1] = (at * f[1] - bt * f[0]) * ph2;
1663  at = a + secxs[1];
1664  bt = b + secys[1];
1665  ct = c + seczs[1];
1666  secxs[2] = (bt * f[2] - ct * f[1]) * ph2;
1667  secys[2] = (ct * f[0] - at * f[2]) * ph2;
1668  seczs[2] = (at * f[1] - bt * f[0]) * ph2;
1669  dxt = h * (a + secxs[2]);
1670  dyt = h * (b + secys[2]);
1671  dzt = h * (c + seczs[2]);
1672  xt = x + dxt;
1673  yt = y + dyt;
1674  zt = z + dzt;
1675  at = a + 2.*secxs[2];
1676  bt = b + 2.*secys[2];
1677  ct = c + 2.*seczs[2];
1678 
1679  est = TMath::Abs(dxt)+TMath::Abs(dyt)+TMath::Abs(dzt);
1680  if (est > 2.*TMath::Abs(h)) {
1681  if (ncut++ > maxcut) break;
1682  h *= khalf;
1683  continue;
1684  }
1685 
1686  xyzt[0] = xt;
1687  xyzt[1] = yt;
1688  xyzt[2] = zt;
1689 
1690  //cmodif: call gufld(xyzt,f) changed into:
1691  TGeoGlobalMagField::Instance()->Field(xyzt,f);
1692 
1693  z = z + (c + (seczs[0] + seczs[1] + seczs[2]) * kthird) * h;
1694  y = y + (b + (secys[0] + secys[1] + secys[2]) * kthird) * h;
1695  x = x + (a + (secxs[0] + secxs[1] + secxs[2]) * kthird) * h;
1696 
1697  secxs[3] = (bt*f[2] - ct*f[1])* ph2;
1698  secys[3] = (ct*f[0] - at*f[2])* ph2;
1699  seczs[3] = (at*f[1] - bt*f[0])* ph2;
1700  a = a+(secxs[0]+secxs[3]+2. * (secxs[1]+secxs[2])) * kthird;
1701  b = b+(secys[0]+secys[3]+2. * (secys[1]+secys[2])) * kthird;
1702  c = c+(seczs[0]+seczs[3]+2. * (seczs[1]+seczs[2])) * kthird;
1703 
1704  est = TMath::Abs(secxs[0]+secxs[3] - (secxs[1]+secxs[2]))
1705  + TMath::Abs(secys[0]+secys[3] - (secys[1]+secys[2]))
1706  + TMath::Abs(seczs[0]+seczs[3] - (seczs[1]+seczs[2]));
1707 
1708  if (est > kdlt && TMath::Abs(h) > 1.e-4) {
1709  if (ncut++ > maxcut) break;
1710  h *= khalf;
1711  continue;
1712  }
1713 
1714  ncut = 0;
1715  // * if too many iterations, go to helix
1716  if (iter++ > maxit) break;
1717 
1718  tl += h;
1719  if (est < kdlt32)
1720  h *= 2.;
1721  cba = 1./ TMath::Sqrt(a*a + b*b + c*c);
1722  vout[0] = x;
1723  vout[1] = y;
1724  vout[2] = z;
1725  vout[3] = cba*a;
1726  vout[4] = cba*b;
1727  vout[5] = cba*c;
1728  rest = step - tl;
1729  if (step < 0.) rest = -rest;
1730  if (rest < 1.e-5*TMath::Abs(step)) return kTRUE;
1731 
1732  } while(1);
1733 
1734  // angle too big, use helix
1735  cout<<"W-AliMFTTrackExtrap::ExtrapOneStepRungekutta: Ruge-Kutta failed: switch to helix"<<endl;
1736 
1737  f1 = f[0];
1738  f2 = f[1];
1739  f3 = f[2];
1740  f4 = TMath::Sqrt(f1*f1+f2*f2+f3*f3);
1741  if (f4 < 1.e-10) {
1742  cout<<"E-AliMFTTrackExtrap::ExtrapOneStepRungekutta: magnetic field at (";
1743  cout<<xyzt[0]<<", "<<xyzt[1]<<", "<<xyzt[2]<<") = "<<f4<<": giving up"<<endl;
1744  return kFALSE;
1745  }
1746  rho = -f4*pinv;
1747  tet = rho * step;
1748 
1749  hnorm = 1./f4;
1750  f1 = f1*hnorm;
1751  f2 = f2*hnorm;
1752  f3 = f3*hnorm;
1753 
1754  hxp[0] = f2*vect[kipz] - f3*vect[kipy];
1755  hxp[1] = f3*vect[kipx] - f1*vect[kipz];
1756  hxp[2] = f1*vect[kipy] - f2*vect[kipx];
1757 
1758  hp = f1*vect[kipx] + f2*vect[kipy] + f3*vect[kipz];
1759 
1760  rho1 = 1./rho;
1761  sint = TMath::Sin(tet);
1762  cost = 2.*TMath::Sin(khalf*tet)*TMath::Sin(khalf*tet);
1763 
1764  g1 = sint*rho1;
1765  g2 = cost*rho1;
1766  g3 = (tet-sint) * hp*rho1;
1767  g4 = -cost;
1768  g5 = sint;
1769  g6 = cost * hp;
1770 
1771  vout[kix] = vect[kix] + g1*vect[kipx] + g2*hxp[0] + g3*f1;
1772  vout[kiy] = vect[kiy] + g1*vect[kipy] + g2*hxp[1] + g3*f2;
1773  vout[kiz] = vect[kiz] + g1*vect[kipz] + g2*hxp[2] + g3*f3;
1774 
1775  vout[kipx] = vect[kipx] + g4*vect[kipx] + g5*hxp[0] + g6*f1;
1776  vout[kipy] = vect[kipy] + g4*vect[kipy] + g5*hxp[1] + g6*f2;
1777  vout[kipz] = vect[kipz] + g4*vect[kipz] + g5*hxp[2] + g6*f3;
1778 
1779  return kTRUE;
1780 }
1781 
TBrowser b
Definition: RunAnaESD.C:12
Double_t GetSlopeY() const
return Y slope
const TMatrixD & GetParameters() const
return track parameters
Track parameters in ALICE dimuon spectrometer.
static Double_t Sagitta(Int_t nVal, Double_t *xVal, Double_t *yVal, Double_t &distL, Double_t &q2)
Bool_t CovariancesExist() const
return kTRUE if the covariance matrix exist, kFALSE if not
static Double_t fgSimpleBValue
! magnetic field value at the centre
void SetZ(Double_t z)
set Z coordinate (cm)
static Double_t GetBendingMomentumFromImpactParam(Double_t impactParam)
static void RecoverTrackParam(Double_t *v3, Double_t Charge, AliMFTTrackParam *trackParam)
static Double_t EnergyLossFluctuation(Double_t pTotal, Double_t pathLength, Double_t rho, Double_t atomicZoverA)
void dZ()
Definition: CalibAlign.C:517
static Bool_t ExtrapToZRungekutta(AliMFTTrackParam *trackParam, Double_t Z)
static Double_t BetheBlochGeant(Double_t bg, Double_t kp0=2.33, Double_t kp1=0.20, Double_t kp2=3.00, Double_t kp3=173e-9, Double_t kp4=0.49848)
static Double_t TotalMomentumEnergyLoss(AliMFTTrackParam *trackParam, Double_t xVtx, Double_t yVtx, Double_t zVtx)
const TMatrixD & GetCovariances() const
static void LinearExtrapToZ(AliMFTTrackParam *trackParam, Double_t zEnd)
static void Cov2CovP(const TMatrixD &param, TMatrixD &cov)
static void CovP2Cov(const TMatrixD &param, TMatrixD &cov)
static void LinearExtrapToZCov(AliMFTTrackParam *trackParam, Double_t zEnd, Bool_t updatePropagator=kFALSE)
void SetInverseTransverseMomentum(Double_t val)
set Inverse Momentum
static const Double_t fgkHelixStepLength
! Step lenght for track extrapolation (used in Helix)
void AddParameters(const TMatrixD &parameters)
add track parameters
void SetCovariances(const TMatrixD &covariances)
static void ConvertTrackParamForExtrap(AliMFTTrackParam *trackParam, Double_t forwardBackward, Double_t *v3)
static Double_t GetMCSAngle2(const AliMFTTrackParam &param, Double_t dZ, Double_t x0)
void SetX(Double_t val)
set X coordinate (cm)
Double_t GetInverseTransverseMomentum() const
return Inverse Momentum
static const Double_t fgkRungeKuttaMaxResidue
! Maximal distance (in Z) to destination to stop the track extrapolation (used in Runge-Kutta) ...
static void ExtrapToVertex(AliMFTTrackParam *trackParam, Double_t xVtx, Double_t yVtx, Double_t zVtx, Double_t errXVtx, Double_t errYVtx)
void sum()
Double_t chi2
Definition: AnalyzeLaser.C:7
static Double_t GetImpactParamFromBendingMomentum(Double_t bendingMomentum)
static Bool_t ExtrapToZHelix(AliMFTTrackParam *trackParam, Double_t Z)
static Bool_t fgFieldON
! kTRUE if the field is switched ON
TGeoManager * gGeoManager
Double_t GetSlopeX() const
return X slope
static void AddMCSEffect(AliMFTTrackParam *param, Double_t dZ, Double_t x0)
static const Int_t fgkMaxStepNumber
! Maximum number of steps for track extrapolation
static Bool_t ExtrapToZ(AliMFTTrackParam *trackParam, Double_t zEnd)
Double_t GetZ() const
return Z coordinate (cm)
void SetSlopeY(Double_t val)
set Y slope
static void ExtrapToVertexUncorrected(AliMFTTrackParam *trackParam, Double_t zVtx)
Double_t P() const
return total momentum
static void CorrectMCSEffectInAbsorber(AliMFTTrackParam *param, Double_t xVtx, Double_t yVtx, Double_t zVtx, Double_t errXVtx, Double_t errYVtx, Double_t absZBeg, Double_t pathLength, Double_t f0, Double_t f1, Double_t f2)
TF1 * f
Definition: interpolTest.C:21
void SetY(Double_t val)
set Y coordinate (cm)
void UpdatePropagator(const TMatrixD &propagator)
static Bool_t GetAbsorberCorrectionParam(Double_t trackXYZIn[3], Double_t trackXYZOut[3], Double_t pTotal, Double_t &pathLength, Double_t &f0, Double_t &f1, Double_t &f2, Double_t &meanRho, Double_t &totalELoss, Double_t &sigmaELoss2)
static void ExtrapToVertexWithoutBranson(AliMFTTrackParam *trackParam, Double_t zVtx)
Constants for the Muon Forward Tracker.
static Double_t CircleRegression(Int_t nVal, Double_t *xVal, Double_t *yVal)
static void CorrectELossEffectInAbsorber(AliMFTTrackParam *param, Double_t eLoss, Double_t sigmaELoss2)
static const Bool_t fgkUseHelix
! Tell whether to use Helix or not (default is Runge-Kutta)
Double_t GetX() const
return X coordinate (cm)
Class holding the parameter of a MFT Standalone Track.
Double_t GetY() const
return Y coordinate (cm)
static Bool_t ExtrapOneStepRungekutta(Double_t charge, Double_t step, const Double_t *vect, Double_t *vout)
static void AddMCSEffectInAbsorber(AliMFTTrackParam *trackParam, Double_t signedPathLength, Double_t f0, Double_t f1, Double_t f2)
static Bool_t ExtrapToZCov(AliMFTTrackParam *trackParam, Double_t zEnd, Bool_t updatePropagator=kFALSE)
void SetParameters(const TMatrixD &parameters)
set track parameters
static void ExtrapOneStepHelix(Double_t charge, Double_t step, const Double_t *vect, Double_t *vout)
static Double_t BetheBloch(Double_t pTotal, Double_t pathLength, Double_t rho, Double_t atomicZ, Double_t atomicZoverA)
static void ExtrapToVertexWithoutELoss(AliMFTTrackParam *trackParam, Double_t xVtx, Double_t yVtx, Double_t zVtx, Double_t errXVtx, Double_t errYVtx)
static Double_t LinearRegression(Int_t nVal, Double_t *xVal, Double_t *yVal, Double_t &p0, Double_t &p1)
static Double_t QuadraticRegression(Int_t nVal, Double_t *xVal, Double_t *yVal, Double_t &p0, Double_t &p1, Double_t &p2)
static void ExtrapOneStepHelix3(Double_t field, Double_t step, const Double_t *vect, Double_t *vout)
void SetSlopeX(Double_t val)
set X slope
class TMatrixT< Double_t > TMatrixD