Line data Source code
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;
42 : /// \cond CLASSIMP
43 12 : ClassImp(AliMFTTrackExtrap) // Class implementation in ROOT context
44 : /// \endcond
45 :
46 12 : 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 : {
58 : /// set field on/off flag;
59 : /// set field at the centre of the MFT
60 0 : const Double_t x[3] = {0.,0.,fgkSimpleBPosition};
61 0 : Double_t b[3] = {0.,0.,0.};
62 0 : TGeoGlobalMagField::Instance()->Field(x,b);
63 0 : cout<<Form("Field = %e %e %e",b[0],b[1],b[2])<<endl;
64 0 : fgSimpleBValue = b[2];
65 0 : fgFieldON = (TMath::Abs(fgSimpleBValue) > 1.e-10) ? kTRUE : kFALSE;
66 : // fgFieldON = kFALSE;
67 0 : }
68 :
69 : //__________________________________________________________________________
70 : Double_t AliMFTTrackExtrap::GetImpactParamFromBendingMomentum(Double_t bendingMomentum)
71 : {
72 : /// Returns impact parameter at vertex in bending plane (cm),
73 : /// from the signed bending momentum "BendingMomentum" in bending plane (GeV/c),
74 : /// using simple values for dipole magnetic field.
75 : /// The sign of "BendingMomentum" is the sign of the charge.
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 0 : }
83 :
84 : //__________________________________________________________________________
85 : Double_t
86 : AliMFTTrackExtrap::GetBendingMomentumFromImpactParam(Double_t impactParam)
87 : {
88 : /// Returns signed bending momentum in bending plane (GeV/c),
89 : /// the sign being the sign of the charge for particles moving forward in Z,
90 : /// from the impact parameter "ImpactParam" at vertex in bending plane (cm),
91 : /// using simple values for dipole magnetic field.
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 0 : }
106 : //__________________________________________________________________________
107 : Double_t AliMFTTrackExtrap::Sagitta(Int_t nVal, Double_t *xVal, Double_t *yVal, Double_t &distL2, Double_t &q2)
108 : {
109 : /// Calculate sagitta of the track
110 : /// Return sagitta
111 0 : cout<<"Sagitta"<<endl;
112 0 : Double_t q0,q1,chi2;
113 : // fit by a polynom of 2nd order
114 0 : chi2 = QuadraticRegression(nVal,xVal ,yVal, q0, q1,q2);
115 0 : cout<<Form("q param = %f %f %f chi2 = %e ",q0, q1,q2,chi2)<<endl;
116 0 : 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 0 : for (int i=0; i<nVal-1; i++) {
127 0 : 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 0 : for (int j = i+1; j<nVal; j++) {
130 0 : Double_t dist = (xVal[i]-xVal[j]) * (xVal[i]-xVal[j]) + (yVal[i]-yVal[j]) * (yVal[i]-yVal[j]);
131 0 : if(dist>maxDist){
132 : i1 = i;
133 : i2 = j;
134 : maxDist = dist;
135 0 : }
136 : }
137 : }
138 :
139 0 : cout<<Form("distMAx = %f distTot =%f i1 = %d i2 =%d",TMath::Sqrt(maxDist),distTot,i1,i2)<<endl;
140 :
141 0 : distL2 = 1.e-2*TMath::Sqrt((xVal[i1]-xVal[i2]) * (xVal[i1]-xVal[i2]) + (yVal[i1]-yVal[i2]) * (yVal[i1]-yVal[i2]) );
142 :
143 0 : Double_t p1 = (yVal[i1]-yVal[i2]) / (xVal[i1]-xVal[i2]);
144 0 : Double_t p0 = yVal[i2] - xVal[i2] * p1;
145 0 : cout<<Form("p param = %f %f ",p0, p1)<<endl;
146 0 : q2 = TMath::Sign(q2,q1*q2*(-(p0 + p1 * (xVal[i1]+xVal[i2])/2.))*(fgSimpleBValue));
147 :
148 :
149 0 : Double_t y12 = q0+q1*xVal[i1]+q2*xVal[i1]*xVal[i1];
150 0 : Double_t y22 = q0+q1*xVal[i2]+q2*xVal[i2]*xVal[i2];
151 0 : cout<<Form("x1, y1 %f %f ",xVal[i1], y12)<<endl;
152 0 : 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 0 : for (int i=0; i<20; i++) {
160 0 : Double_t step =TMath::Abs(xVal[i1]-xVal[i2])/20.;
161 0 : Double_t xmiddle = xVal[i1] + i*step;
162 :
163 0 : Double_t y1 = p0 + p1 * xmiddle;
164 0 : Double_t p1perp = -1./p1;
165 0 : Double_t p0perp = p0 + xmiddle *(p1-p1perp);
166 : //cout<<Form("p param perp = %f %f ",p0perp, p1perp)<<endl;
167 :
168 0 : Double_t xmax = (p1perp - q1 + TMath::Sqrt(p1perp*p1perp - 2*p1perp*q1 + q1*q1 + 4*p0perp*q2 - 4*q0*q2))/(2*q2);
169 0 : 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 0 : if (TMath::Abs(xmax2-xmiddle) < TMath::Abs(xmax-xmiddle)) xmax = xmax2;
173 :
174 :
175 0 : Double_t y2 = q0 + q1 * xmax + q2 * xmax * xmax;
176 :
177 0 : 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 0 : sagitta = TMath::Sign(sagitta,q1*q2*(-y1)*(fgSimpleBValue));
182 :
183 0 : if(TMath::Abs(sagitta)>TMath::Abs(maxSagitta)) maxSagitta = sagitta;
184 :
185 : }
186 0 : 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;
195 : //// cout<<Form("x middle = %f ",xmiddle)<<endl;
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);
204 : //// cout<<Form("xmax = %f xmax2 %f",xmax,xmax2)<<endl;
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 : //
213 : //// cout<<Form("sagitta = %f Bz = %f",sagitta,fgSimpleBValue)<<endl;
214 : //
215 : // sagitta = TMath::Sign(sagitta,q1*q2*(-y1)*(fgSimpleBValue));
216 : //
217 : //// cout<<Form("Signed sagitta = %e ",sagitta)<<endl;
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 0 : return maxSagitta;
224 0 : }
225 :
226 : //__________________________________________________________________________
227 : Double_t AliMFTTrackExtrap::LinearRegression(Int_t nVal, Double_t *xVal, Double_t *yVal, Double_t &p0, Double_t &p1)
228 : {
229 : /// Perform a Linear Regression
230 : /// Return Chi2
231 : Double_t meanx =0, meany=0.;
232 :
233 0 : for (Int_t i = 0; i< nVal; i++) {
234 0 : meanx =(meanx*i+xVal[i])/(i+1);
235 0 : meany =(meany*i+yVal[i])/(i+1);
236 : }
237 : Double_t cov_xy = 0., var_x=0., var_y=0.;
238 0 : for (Int_t i = 0; i< nVal; i++) {
239 0 : var_x += (xVal[i] -meanx) * (xVal[i] -meanx);
240 0 : var_y += (yVal[i] -meany) * (yVal[i] -meany);
241 0 : cov_xy += (xVal[i] -meanx) * (yVal[i] -meany);
242 : }
243 0 : if(var_x<1.e-6) return 0.;
244 0 : p1 = cov_xy/var_x;
245 0 : p0 = meany - p1 * meanx;
246 :
247 : Double_t chi2 = 0.;
248 : Double_t yest=0.;
249 0 : for (Int_t i = 0; i< nVal; i++) {
250 0 : yest = xVal[i]*p1+p0;
251 0 : chi2 += (yest-yVal[i]) * (yest-yVal[i]);
252 : }
253 0 : chi2 /= var_y;
254 :
255 :
256 : return chi2;
257 0 : }
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 : {
261 : /// Perform a Quadratic Regression
262 : /// Assume same error on all clusters = 1
263 : /// Return ~ Chi2
264 :
265 0 : TMatrixD y(nVal,1);
266 0 : TMatrixD x(nVal,3);
267 0 : TMatrixD xtrans(3,nVal);
268 :
269 0 : for (int i=0; i<nVal; i++) {
270 0 : y(i,0) = yVal[i];
271 0 : x(i,0) = 1.;
272 0 : x(i,1) = xVal[i];
273 0 : x(i,2) = xVal[i]*xVal[i];
274 0 : xtrans(0,i) = 1.;
275 0 : xtrans(1,i) = xVal[i];
276 0 : xtrans(2,i) = xVal[i]*xVal[i];
277 : }
278 0 : TMatrixD tmp(xtrans,TMatrixD::kMult,x);
279 0 : tmp.Invert();
280 :
281 0 : TMatrixD tmp2(xtrans,TMatrixD::kMult,y);
282 0 : TMatrixD b(tmp,TMatrixD::kMult,tmp2);
283 :
284 0 : p0 = b(0,0);
285 0 : p1 = b(1,0);
286 0 : p2 = b(2,0);
287 :
288 : // chi2 = (y-xb)^t . W . (y-xb)
289 0 : TMatrixD tmp3(x,TMatrixD::kMult,b);
290 0 : TMatrixD tmp4(y,TMatrixD::kMinus,tmp3);
291 0 : TMatrixD chi2(tmp4,TMatrixD::kTransposeMult,tmp4);
292 :
293 :
294 0 : return chi2(0,0);
295 0 : }
296 : //__________________________________________________________________________
297 : Double_t AliMFTTrackExtrap::CircleRegression(Int_t nVal, Double_t *xVal, Double_t *yVal)
298 : {
299 : /// Perform a Circular Regression
300 : /// Assume same error on all clusters = 1
301 : /// Return Radius evaluation
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 0 : for (int i=0; i<nVal; i++) {
306 0 : xi = xVal[i]/100.;
307 0 : yi = yVal[i]/100.;
308 0 : ri = TMath::Sqrt(xi*xi + yi*yi);
309 : // xi /= ri*ri;
310 : // yi /= ri*ri;
311 :
312 0 : sumxi += xi;
313 0 : sumyi += yi;
314 0 : sumxi2 += xi*xi;
315 0 : sumyi2 += yi*yi;
316 0 : sumxi3 += xi*xi*xi;
317 0 : sumyi3 += yi*yi*yi;
318 0 : sumxiyi += xi*yi;
319 0 : sumxi2yi += xi*xi*yi;
320 0 : sumxiyi2 += xi*yi*yi;
321 : }
322 :
323 0 : Double_t A = nVal * sumxi2 - sumxi*sumxi;
324 0 : Double_t B = nVal * sumxiyi - sumxi*sumyi;
325 0 : Double_t C = nVal * sumyi2 - sumyi*sumyi;
326 0 : Double_t D = 0.5*(nVal*sumxiyi2 -sumxi*sumyi2 +nVal*sumxi3 -sumxi*sumxi2);
327 0 : Double_t E = 0.5*(nVal*sumxi2yi -sumyi*sumxi2 +nVal*sumyi3 -sumyi*sumyi2);
328 :
329 0 : Double_t aM = (D*C-B*E) / (A*C-B*B);
330 0 : Double_t bM = (A*E-B*D) / (A*C-B*B);
331 :
332 : Double_t rM = 0.;
333 : Double_t rK = 0.;
334 :
335 0 : for (int i=0; i<nVal; i++) {
336 0 : xi = xVal[i]/100.;
337 0 : yi = yVal[i]/100.;
338 :
339 0 : rM += TMath::Sqrt( (xi-aM)*(xi-aM) + (yi-bM)*(yi-bM) );
340 0 : rK += ((xi-aM)*(xi-aM) + (yi-bM)*(yi-bM) );
341 : }
342 0 : rM /= nVal;
343 0 : rK = TMath::Sqrt(rK/nVal);
344 :
345 0 : 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 0 : return (rM+rK)/2.;
348 : }
349 :
350 : //__________________________________________________________________________
351 : void AliMFTTrackExtrap::LinearExtrapToZ(AliMFTTrackParam* trackParam, Double_t zEnd)
352 : {
353 : /// Track parameters linearly extrapolated to the plane at "zEnd".
354 : /// On return, results from the extrapolation are updated in trackParam.
355 :
356 0 : if (trackParam->GetZ() == zEnd) return; // nothing to be done if same z
357 :
358 : // Compute track parameters
359 0 : Double_t dZ = zEnd - trackParam->GetZ();
360 :
361 0 : trackParam->SetX(trackParam->GetX() + trackParam->GetSlopeX() * dZ);
362 0 : trackParam->SetY(trackParam->GetY() + trackParam->GetSlopeY() * dZ);
363 0 : trackParam->SetZ(zEnd);
364 0 : }
365 :
366 : //__________________________________________________________________________
367 : void AliMFTTrackExtrap::LinearExtrapToZCov(AliMFTTrackParam* trackParam, Double_t zEnd, Bool_t updatePropagator)
368 : {
369 : /// Track parameters and their covariances linearly extrapolated to the plane at "zEnd".
370 : /// On return, results from the extrapolation are updated in trackParam.
371 :
372 0 : 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 0 : if (!trackParam->CovariancesExist()) {
376 0 : cout<<"W-AliMUONTrackExtrap::LinearExtrapToZCov: Covariance matrix does not exist"<<endl;
377 : // Extrapolate linearly track parameters to "zEnd"
378 0 : LinearExtrapToZ(trackParam,zEnd);
379 0 : return;
380 : }
381 :
382 : // Compute track parameters
383 0 : Double_t dZ = zEnd - trackParam->GetZ();
384 0 : trackParam->SetX(trackParam->GetX() + trackParam->GetSlopeX() * dZ);
385 0 : trackParam->SetY(trackParam->GetY() + trackParam->GetSlopeY() * dZ);
386 0 : trackParam->SetZ(zEnd);
387 :
388 : // Calculate the jacobian related to the track parameters linear extrapolation to "zEnd"
389 0 : TMatrixD jacob(5,5);
390 0 : jacob.UnitMatrix();
391 0 : jacob(0,2) = dZ;
392 0 : jacob(1,3) = dZ;
393 :
394 : // Extrapolate track parameter covariances to "zEnd"
395 0 : TMatrixD tmp(trackParam->GetCovariances(),TMatrixD::kMultTranspose,jacob);
396 0 : TMatrixD tmp2(jacob,TMatrixD::kMult,tmp);
397 0 : trackParam->SetCovariances(tmp2);
398 :
399 : // Update the propagator if required
400 0 : if (updatePropagator) trackParam->UpdatePropagator(jacob);
401 0 : }
402 :
403 : //__________________________________________________________________________
404 : Bool_t AliMFTTrackExtrap::ExtrapToZ(AliMFTTrackParam* trackParam, Double_t zEnd)
405 : {
406 : /// Interface to track parameter extrapolation to the plane at "Z" using Helix or Rungekutta algorithm.
407 : /// On return, the track parameters resulting from the extrapolation are updated in trackParam.
408 : //return AliMFTTrackExtrap::ExtrapToZHelix(trackParam,zEnd);
409 0 : if (!fgFieldON) {
410 0 : AliMFTTrackExtrap::LinearExtrapToZ(trackParam,zEnd);
411 0 : return kTRUE;
412 : }
413 0 : else if (fgkUseHelix) return AliMFTTrackExtrap::ExtrapToZHelix(trackParam,zEnd);
414 : else return AliMFTTrackExtrap::ExtrapToZRungekutta(trackParam,zEnd);
415 0 : }
416 :
417 : //__________________________________________________________________________
418 : Bool_t AliMFTTrackExtrap::ExtrapToZHelix(AliMFTTrackParam* trackParam, Double_t zEnd)
419 : {
420 : // cout<<"I-AliMFTTrackExtrap::ExtrapToZHelix: Entering ------ "<<endl;
421 : /// Track parameter extrapolation to the plane at "Z" using Helix algorithm.
422 : /// On return, the track parameters resulting from the extrapolation are updated in trackParam.
423 0 : if (trackParam->GetZ() == zEnd) return kTRUE; // nothing to be done if same Z
424 : Double_t forwardBackward; // +1 if forward, -1 if backward
425 0 : if (zEnd < trackParam->GetZ()) forwardBackward = 1.0; // spectro. z<0
426 : else forwardBackward = -1.0;
427 0 : 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 0 : ConvertTrackParamForExtrap(trackParam, forwardBackward, v3);
432 : // sign of charge (sign of fInverseBendingMomentum if forward motion)
433 : // must be changed if backward extrapolation
434 0 : 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 0 : while (((-forwardBackward * (v3[2] - zEnd)) <= 0.0) && (stepNumber < fgkMaxStepNumber)) { // spectro. z<0
441 0 : stepNumber++;
442 0 : ExtrapOneStepHelix(chargeExtrap, fgkHelixStepLength, v3, v3New);
443 : // cout<<" (-forwardBackward * (v3New[2] - zEnd)) = "<<(-forwardBackward * (v3New[2] - zEnd))<<endl;
444 :
445 0 : if ((-forwardBackward * (v3New[2] - zEnd)) > 0.0) break; // one is beyond Z spectro. z<0
446 : // better use TArray ????
447 0 : for (i3 = 0; i3 < 7; i3++) {
448 : // cout<<" v3New["<<i3<< "] = "<<v3New[i3]<<endl;
449 0 : 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 0 : Double_t dZ12 = v3New[2] - v3[2]; // 1->2
456 0 : if (TMath::Abs(dZ12) > 0) {
457 0 : Double_t dZ1i = zEnd - v3[2]; // 1-i
458 0 : Double_t dZi2 = v3New[2] - zEnd; // i->2
459 0 : Double_t xPrime = (v3New[0] - v3[0]) / dZ12;
460 0 : Double_t xSecond = ((v3New[3] / v3New[5]) - (v3[3] / v3[5])) / dZ12;
461 0 : Double_t yPrime = (v3New[1] - v3[1]) / dZ12;
462 0 : Double_t ySecond = ((v3New[4] / v3New[5]) - (v3[4] / v3[5])) / dZ12;
463 0 : v3[0] = v3[0] + xPrime * dZ1i - 0.5 * xSecond * dZ1i * dZi2; // X
464 0 : v3[1] = v3[1] + yPrime * dZ1i - 0.5 * ySecond * dZ1i * dZi2; // Y
465 0 : v3[2] = zEnd; // Z
466 0 : Double_t xPrimeI = xPrime - 0.5 * xSecond * (dZi2 - dZ1i);
467 0 : Double_t yPrimeI = yPrime - 0.5 * ySecond * (dZi2 - dZ1i);
468 : // (PX, PY, PZ)/PTOT assuming forward motion
469 0 : v3[5] = 1.0 / TMath::Sqrt(1.0 + xPrimeI * xPrimeI + yPrimeI * yPrimeI); // PZ/PTOT
470 0 : v3[3] = xPrimeI * v3[5]; // PX/PTOT
471 0 : v3[4] = yPrimeI * v3[5]; // PY/PTOT
472 0 : } else {
473 0 : cout<<"W-AliMFTTrackExtrap::ExtrapToZHelix: Extrap. to Z not reached, Z = "<<zEnd<<endl;
474 : }
475 : // Recover track parameters (charge back for forward motion)
476 0 : RecoverTrackParam(v3, chargeExtrap * forwardBackward, trackParam);
477 : return kTRUE;
478 0 : }
479 :
480 : //__________________________________________________________________________
481 : Bool_t AliMFTTrackExtrap::ExtrapToZRungekutta(AliMFTTrackParam* trackParam, Double_t zEnd)
482 : {
483 : /// Track parameter extrapolation to the plane at "Z" using Rungekutta algorithm.
484 : /// On return, the track parameters resulting from the extrapolation are updated in trackParam.
485 0 : if (trackParam->GetZ() == zEnd) return kTRUE; // nothing to be done if same Z
486 : Double_t forwardBackward; // +1 if forward, -1 if backward
487 0 : 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 0 : Double_t chargeExtrap = forwardBackward * TMath::Sign(Double_t(1.0), trackParam->GetInverseTransverseMomentum());
492 0 : 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 0 : Double_t residue = zEnd - trackParam->GetZ();
498 : Bool_t uturn = kFALSE;
499 : Bool_t trackingFailed = kFALSE;
500 : Bool_t tooManyStep = kFALSE;
501 0 : while (TMath::Abs(residue) > fgkRungeKuttaMaxResidue && stepNumber <= fgkMaxStepNumber) {
502 :
503 0 : dZ = zEnd - trackParam->GetZ();
504 : // step lenght assuming linear trajectory
505 0 : step = dZ * TMath::Sqrt(1.0 + trackParam->GetSlopeY()*trackParam->GetSlopeY() +
506 0 : trackParam->GetSlopeX()*trackParam->GetSlopeX());
507 0 : ConvertTrackParamForExtrap(trackParam, forwardBackward, v3);
508 :
509 0 : do { // reduce step lenght while zEnd oversteped
510 0 : if (stepNumber > fgkMaxStepNumber) {
511 0 : cout<<"W-AliMFTTrackExtrap::ExtrapToZRungekutta: Too many trials: "<<stepNumber<<endl;
512 : tooManyStep = kTRUE;
513 0 : break;
514 : }
515 0 : stepNumber ++;
516 0 : step = TMath::Abs(step);
517 0 : if (!AliMFTTrackExtrap::ExtrapOneStepRungekutta(chargeExtrap,step,v3,v3New)) {
518 : trackingFailed = kTRUE;
519 0 : break;
520 : }
521 0 : residue = zEnd - v3New[2];
522 0 : step *= dZ/(v3New[2]-trackParam->GetZ());
523 0 : } while (residue*dZ < 0 && TMath::Abs(residue) > fgkRungeKuttaMaxResidue);
524 :
525 0 : if (trackingFailed) break;
526 0 : else if (v3New[5]*v3[5] < 0) { // the track turned around
527 0 : cout<<"W-AliMFTTrackExtrap::ExtrapToZRungekutta: The track turned around"<<endl;
528 : uturn = kTRUE;
529 0 : break;
530 0 : } else RecoverTrackParam(v3New, chargeExtrap * forwardBackward, trackParam);
531 :
532 : }
533 :
534 : // terminate the extropolation with a straight line up to the exact "zEnd" value
535 : /// todo : change that
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 0 : }
563 :
564 : //__________________________________________________________________________
565 : void AliMFTTrackExtrap::ConvertTrackParamForExtrap(AliMFTTrackParam* trackParam, Double_t forwardBackward, Double_t *v3)
566 : {
567 : /// Set vector of Geant3 parameters pointed to by "v3" from track parameters in trackParam.
568 : /// Since AliMFTTrackParam is only geometry, one uses "forwardBackward"
569 : /// to know whether the particle is going forward (+1) or backward (-1).
570 0 : v3[0] = trackParam->GetX(); // X
571 0 : v3[1] = trackParam->GetY(); // Y
572 0 : v3[2] = trackParam->GetZ(); // Z
573 0 : Double_t slopeX = trackParam->GetSlopeX() ;
574 0 : Double_t slopeY = trackParam->GetSlopeY() ;
575 0 : Double_t slope2 = TMath::Sqrt(1.+slopeX*slopeX +slopeY*slopeY);
576 0 : Double_t pt = TMath::Abs(1.0 / trackParam->GetInverseTransverseMomentum());
577 0 : v3[6] = pt*slope2/TMath::Sqrt(slopeX*slopeX +slopeY*slopeY); // PTOT
578 :
579 0 : v3[5] = -forwardBackward / slope2; // PZ/PTOT spectro. z<0
580 0 : v3[3] = slopeX / slope2; // PX/PTOT
581 0 : 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 0 : }
595 :
596 : //__________________________________________________________________________
597 : void AliMFTTrackExtrap::RecoverTrackParam(Double_t *v3, Double_t charge, AliMFTTrackParam* trackParam)
598 : {
599 : /// Set track parameters in trackParam from Geant3 parameters pointed to by "v3",
600 : /// assumed to be calculated for forward motion in Z.
601 : /// "InverseBendingMomentum" is signed with "charge".
602 0 : trackParam->SetX(v3[0]); // X
603 0 : trackParam->SetY(v3[1]); // Y
604 0 : trackParam->SetZ(v3[2]); // Z
605 0 : Double_t pt = v3[6]*TMath::Sqrt(1. - v3[5]*v3[5]);
606 0 : trackParam->SetInverseTransverseMomentum(charge/pt);
607 0 : trackParam->SetSlopeY(v3[4]/v3[5]);
608 0 : 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 0 : }
619 :
620 : //__________________________________________________________________________
621 : Bool_t AliMFTTrackExtrap::ExtrapToZCov(AliMFTTrackParam* trackParam, Double_t zEnd, Bool_t updatePropagator)
622 : {
623 : /// Track parameters and their covariances extrapolated to the plane at "zEnd".
624 : /// On return, results from the extrapolation are updated in trackParam.
625 0 : if (trackParam->GetZ() == zEnd) return kTRUE; // nothing to be done if same z
626 :
627 0 : if (!fgFieldON) { // linear extrapolation if no magnetic field
628 0 : AliMFTTrackExtrap::LinearExtrapToZCov(trackParam,zEnd,updatePropagator);
629 0 : return kTRUE;
630 : }
631 :
632 : // No need to propagate the covariance matrix if it does not exist
633 0 : if (!trackParam->CovariancesExist()) {
634 0 : cout<<"W-AliMFTTrackExtrap::ExtrapToZCov: Covariance matrix does not exist"<<endl;
635 : // Extrapolate track parameters to "zEnd"
636 0 : return ExtrapToZ(trackParam,zEnd);
637 : }
638 :
639 : // Save the actual track parameters
640 0 : AliMFTTrackParam trackParamSave(*trackParam);
641 0 : TMatrixD paramSave(trackParamSave.GetParameters());
642 0 : Double_t zBegin = trackParamSave.GetZ();
643 :
644 : // Get reference to the parameter covariance matrix
645 0 : const TMatrixD& kParamCov = trackParam->GetCovariances();
646 :
647 : // Extrapolate track parameters to "zEnd"
648 : // Do not update the covariance matrix if the extrapolation failed
649 0 : if (!ExtrapToZ(trackParam,zEnd)) return kFALSE;
650 :
651 : // Get reference to the extrapolated parameters
652 0 : const TMatrixD& extrapParam = trackParam->GetParameters();
653 :
654 : // Calculate the jacobian related to the track parameters extrapolation to "zEnd"
655 : Bool_t extrapStatus = kTRUE;
656 0 : TMatrixD jacob(5,5);
657 0 : jacob.Zero();
658 0 : TMatrixD dParam(5,1);
659 0 : Double_t direction[5] = {-1.,-1.,1.,1.,-1.};
660 0 : for (Int_t i=0; i<5; i++) {
661 : // Skip jacobian calculation for parameters with no associated error
662 0 : if (kParamCov(i,i) <= 0.) continue;
663 :
664 : // Small variation of parameter i only
665 0 : for (Int_t j=0; j<5; j++) {
666 0 : if (j==i) {
667 0 : dParam(j,0) = TMath::Sqrt(kParamCov(i,i));
668 0 : dParam(j,0) *= TMath::Sign(1.,direction[j]*paramSave(j,0)); // variation always in the same direction
669 0 : } else dParam(j,0) = 0.;
670 : }
671 : // Set new parameters
672 0 : trackParamSave.SetParameters(paramSave);
673 0 : trackParamSave.AddParameters(dParam);
674 0 : trackParamSave.SetZ(zBegin);
675 : // Extrapolate new track parameters to "zEnd"
676 0 : if (!ExtrapToZ(&trackParamSave,zEnd)) {
677 0 : cout<<"W-AliMFTTrackExtrap::ExtrapToZCov: Bad covariance matrix"<<endl;
678 : extrapStatus = kFALSE;
679 0 : }
680 :
681 : // Calculate the jacobian
682 0 : TMatrixD jacobji(trackParamSave.GetParameters(),TMatrixD::kMinus,extrapParam);
683 :
684 0 : jacobji *= 1. / dParam(i,0);
685 0 : jacob.SetSub(0,i,jacobji);
686 0 : }
687 0 : cout<<"jacob"<<endl;
688 0 : jacob.Print();
689 :
690 : // Extrapolate track parameter covariances to "zEnd"
691 0 : cout<<"Initial Cov MAtrix "<<endl;
692 0 : kParamCov.Print();
693 0 : TMatrixD tmp(kParamCov,TMatrixD::kMultTranspose,jacob);
694 0 : TMatrixD tmp2(jacob,TMatrixD::kMult,tmp);
695 0 : cout<<"Extrapolated Cov MAtrix "<<endl;
696 0 : tmp2.Print();
697 :
698 0 : trackParam->SetCovariances(tmp2);
699 : // Update the propagator if required
700 0 : if (updatePropagator) trackParam->UpdatePropagator(jacob);
701 :
702 : // return extrapStatus;
703 : return kTRUE;
704 0 : }
705 :
706 : ////__________________________________________________________________________
707 : void AliMFTTrackExtrap::AddMCSEffectInAbsorber(AliMFTTrackParam* param, Double_t signedPathLength, Double_t f0, Double_t f1, Double_t f2)
708 : {
709 : /// Add to the track parameter covariances the effects of multiple Coulomb scattering
710 : /// signedPathLength must have the sign of (zOut - zIn) where all other parameters are assumed to be given at zOut.
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 0 : }
749 :
750 : //__________________________________________________________________________
751 : void AliMFTTrackExtrap::CorrectMCSEffectInAbsorber(AliMFTTrackParam* param,
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 0 : }
820 :
821 : //__________________________________________________________________________
822 : void AliMFTTrackExtrap::CorrectELossEffectInAbsorber(AliMFTTrackParam* param, Double_t eLoss, Double_t sigmaELoss2)
823 : {
824 : /// Correct parameters for energy loss and add energy loss fluctuation effect to covariances
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 0 : }
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 : {
857 : /// Parameters used to correct for Multiple Coulomb Scattering and energy loss in absorber
858 : /// Calculated assuming a linear propagation from trackXYZIn to trackXYZOut (order is important)
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 0 : pathLength = 0.;
868 0 : f0 = 0.;
869 0 : f1 = 0.;
870 0 : f2 = 0.;
871 0 : meanRho = 0.;
872 0 : totalELoss = 0.;
873 0 : sigmaELoss2 = 0.;
874 :
875 : // Check whether the geometry is available
876 0 : if (!gGeoManager) {
877 0 : cout<<"E-AliMFTTrackExtrap::GetAbsorberCorrectionParam: no TGeo"<<endl;
878 0 : return kFALSE;
879 : }
880 :
881 : // Initialize starting point and direction
882 0 : pathLength = TMath::Sqrt((trackXYZOut[0] - trackXYZIn[0])*(trackXYZOut[0] - trackXYZIn[0])+
883 0 : (trackXYZOut[1] - trackXYZIn[1])*(trackXYZOut[1] - trackXYZIn[1])+
884 0 : (trackXYZOut[2] - trackXYZIn[2])*(trackXYZOut[2] - trackXYZIn[2]));
885 0 : if (pathLength < TGeoShape::Tolerance()) return kFALSE;
886 0 : Double_t b[3];
887 0 : b[0] = (trackXYZOut[0] - trackXYZIn[0]) / pathLength;
888 0 : b[1] = (trackXYZOut[1] - trackXYZIn[1]) / pathLength;
889 0 : b[2] = (trackXYZOut[2] - trackXYZIn[2]) / pathLength;
890 0 : TGeoNode *currentnode = gGeoManager->InitTrack(trackXYZIn, b);
891 0 : if (!currentnode) {
892 0 : cout<<"E-AliMFTTrackExtrap::GetAbsorberCorrectionParam: start point out of geometry"<<endl;
893 0 : 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 0 : Double_t remainingPathLength = pathLength;
904 : Double_t sigmaELoss = 0.;
905 0 : Double_t zB = trackXYZIn[2];
906 : Double_t zE, dzB, dzE;
907 0 : do {
908 : // Get material properties
909 0 : TGeoMaterial *material = currentnode->GetVolume()->GetMedium()->GetMaterial();
910 0 : rho = material->GetDensity();
911 0 : x0 = material->GetRadLen();
912 0 : atomicA = material->GetA();
913 0 : atomicZ = material->GetZ();
914 0 : if(material->IsMixture()){
915 0 : TGeoMixture * mixture = (TGeoMixture*)material;
916 : atomicZoverA = 0.;
917 : Double_t sum = 0.;
918 0 : for (Int_t iel=0;iel<mixture->GetNelements();iel++){
919 0 : sum += mixture->GetWmixt()[iel];
920 0 : atomicZoverA += mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel];
921 : }
922 0 : atomicZoverA/=sum;
923 0 : }
924 0 : else atomicZoverA = atomicZ/atomicA;
925 :
926 : // Get path length within this material
927 0 : gGeoManager->FindNextBoundary(remainingPathLength);
928 0 : 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 0 : if (localPathLength >= remainingPathLength) localPathLength = remainingPathLength;
931 : else {
932 0 : currentnode = gGeoManager->Step();
933 0 : if (!currentnode) {
934 0 : cout<<"E-AliMFTTrackExtrap::GetAbsorberCorrectionParam: navigation failed"<<endl;
935 0 : f0 = f1 = f2 = meanRho = totalELoss = sigmaELoss2 = 0.;
936 0 : return kFALSE;
937 : }
938 0 : if (!gGeoManager->IsEntering()) {
939 : // make another small step to try to enter in new absorber slice
940 0 : gGeoManager->SetStep(0.001);
941 0 : currentnode = gGeoManager->Step();
942 0 : if (!gGeoManager->IsEntering() || !currentnode) {
943 0 : cout<<"E-AliMFTTrackExtrap::GetAbsorberCorrectionParam: navigation failed"<<endl;
944 0 : f0 = f1 = f2 = meanRho = totalELoss = sigmaELoss2 = 0.;
945 0 : return kFALSE;
946 : }
947 0 : localPathLength += 0.001;
948 0 : }
949 : }
950 :
951 : // calculate absorber's parameters
952 0 : zE = b[2] * localPathLength + zB;
953 0 : dzB = zB - trackXYZIn[2];
954 0 : dzE = zE - trackXYZIn[2];
955 0 : f0 += localPathLength / x0;
956 0 : f1 += (dzE*dzE - dzB*dzB) / b[2] / b[2] / x0 / 2.;
957 0 : f2 += (dzE*dzE*dzE - dzB*dzB*dzB) / b[2] / b[2] / b[2] / x0 / 3.;
958 0 : meanRho += localPathLength * rho;
959 0 : totalELoss += BetheBloch(pTotal, localPathLength, rho, atomicZ, atomicZoverA);
960 0 : sigmaELoss += EnergyLossFluctuation(pTotal, localPathLength, rho, atomicZoverA);
961 :
962 : // prepare next step
963 : zB = zE;
964 0 : remainingPathLength -= localPathLength;
965 0 : } while (remainingPathLength > TGeoShape::Tolerance());
966 :
967 0 : meanRho /= pathLength;
968 0 : sigmaELoss2 = sigmaELoss*sigmaELoss;
969 :
970 0 : return kTRUE;
971 0 : }
972 :
973 : //__________________________________________________________________________
974 : Double_t AliMFTTrackExtrap::GetMCSAngle2(const AliMFTTrackParam& param, Double_t dZ, Double_t x0)
975 : {
976 : /// Return the angular dispersion square due to multiple Coulomb scattering
977 : /// through a material of thickness "dZ" and of radiation length "x0"
978 : /// assuming linear propagation and using the small angle approximation.
979 0 : 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 : {
1000 : /// Add to the track parameter covariances the effects of multiple Coulomb scattering
1001 : /// through a material of thickness "Abs(dZ)" and of radiation length "x0"
1002 : /// assuming linear propagation and using the small angle approximation.
1003 : /// dZ = zOut - zIn (sign is important) and "param" is assumed to be given zOut.
1004 : /// If x0 <= 0., assume dZ = pathLength/x0 and consider the material thickness as negligible.
1005 0 : Double_t slopeX = param->GetSlopeX();
1006 0 : Double_t slopeY = param->GetSlopeY();
1007 0 : Double_t slope2 = slopeX*slopeX+slopeY*slopeY;
1008 :
1009 0 : Double_t inversePt = TMath::Abs(param->GetInverseTransverseMomentum());
1010 :
1011 0 : Double_t inverseTotalMomentum2 = inversePt*inversePt / (1. + 1./slope2 );
1012 : // Path length in the material
1013 0 : Double_t signedPathLength = dZ * TMath::Sqrt(1.0 + slope2);
1014 0 : 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 0 : Double_t theta02 = 0.0136 / velo * (1 + 0.038 * TMath::Log(pathLengthOverX0));
1019 0 : theta02 *= theta02 * inverseTotalMomentum2 * pathLengthOverX0;
1020 :
1021 0 : Double_t varCoor = (x0 > 0.) ? signedPathLength * signedPathLength * theta02 / 3. : 0.;
1022 : Double_t varSlop = theta02;
1023 0 : Double_t covCorrSlope = (x0 > 0.) ? signedPathLength * theta02/ 2. : 0.;
1024 :
1025 0 : cout<<Form("theta02=%e inverseTotalMomentum2=%e signedPathLength=%e pathLengthOverX0=%e ",theta02,inverseTotalMomentum2,signedPathLength,pathLengthOverX0 )<<endl;
1026 :
1027 0 : cout<<Form("dz=%e x0=%e varCoor=%e varSlop=%e covCorrSlope=%e ",dZ,x0,varCoor,varSlop,covCorrSlope )<<endl;
1028 : // Set MCS covariance matrix
1029 0 : TMatrixD newParamCov(param->GetCovariances());
1030 0 : cout<<"Covariance avant MCS"<<endl;
1031 0 : newParamCov.Print();
1032 : // Non bending plane
1033 0 : newParamCov(0,0) += varCoor; newParamCov(0,2) += covCorrSlope;
1034 0 : newParamCov(2,0) += covCorrSlope; newParamCov(2,2) += varSlop;
1035 : // Bending plane
1036 0 : newParamCov(1,1) += varCoor; newParamCov(1,3) += covCorrSlope;
1037 0 : newParamCov(3,1) += covCorrSlope; newParamCov(3,3) += varSlop;
1038 :
1039 0 : cout<<"Covariance apres MCS"<<endl;
1040 0 : 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 0 : param->SetCovariances(newParamCov);
1060 0 : }
1061 :
1062 : //__________________________________________________________________________
1063 : void AliMFTTrackExtrap::ExtrapToVertex(AliMFTTrackParam* trackParam,
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 : {
1068 : /// Main method for extrapolation to the vertex:
1069 : /// Returns the track parameters and covariances resulting from the extrapolation of the current trackParam
1070 : /// Changes parameters and covariances according to multiple scattering and energy loss corrections:
1071 : /// if correctForMCS=kTRUE: compute parameters using Branson correction and add correction resolution to covariances
1072 : /// if correctForMCS=kFALSE: add parameter dispersion due to MCS in parameter covariances
1073 : /// if correctForEnergyLoss=kTRUE: correct parameters for energy loss and add energy loss fluctuation to covariances
1074 : /// if correctForEnergyLoss=kFALSE: do nothing about energy loss
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 0 : }
1177 :
1178 : //__________________________________________________________________________
1179 : void AliMFTTrackExtrap::ExtrapToVertex(AliMFTTrackParam* trackParam,
1180 : Double_t xVtx, Double_t yVtx, Double_t zVtx,
1181 : Double_t errXVtx, Double_t errYVtx)
1182 : {
1183 : /// Extrapolate track parameters to vertex, corrected for multiple scattering and energy loss effects
1184 : /// Add branson correction resolution and energy loss fluctuation to parameter covariances
1185 0 : ExtrapToVertex(trackParam, xVtx, yVtx, zVtx, errXVtx, errYVtx, kTRUE, kTRUE);
1186 0 : }
1187 :
1188 : //__________________________________________________________________________
1189 : void AliMFTTrackExtrap::ExtrapToVertexWithoutELoss(AliMFTTrackParam* trackParam,
1190 : Double_t xVtx, Double_t yVtx, Double_t zVtx,
1191 : Double_t errXVtx, Double_t errYVtx)
1192 : {
1193 : /// Extrapolate track parameters to vertex, corrected for multiple scattering effects only
1194 : /// Add branson correction resolution to parameter covariances
1195 0 : ExtrapToVertex(trackParam, xVtx, yVtx, zVtx, errXVtx, errYVtx, kTRUE, kFALSE);
1196 0 : }
1197 :
1198 : //__________________________________________________________________________
1199 : void AliMFTTrackExtrap::ExtrapToVertexWithoutBranson(AliMFTTrackParam* trackParam, Double_t zVtx)
1200 : {
1201 : /// Extrapolate track parameters to vertex, corrected for energy loss effects only
1202 : /// Add dispersion due to multiple scattering and energy loss fluctuation to parameter covariances
1203 0 : ExtrapToVertex(trackParam, 0., 0., zVtx, 0., 0., kFALSE, kTRUE);
1204 0 : }
1205 :
1206 : //__________________________________________________________________________
1207 : void AliMFTTrackExtrap::ExtrapToVertexUncorrected(AliMFTTrackParam* trackParam, Double_t zVtx)
1208 : {
1209 : /// Extrapolate track parameters to vertex without multiple scattering and energy loss corrections
1210 : /// Add dispersion due to multiple scattering to parameter covariances
1211 0 : ExtrapToVertex(trackParam, 0., 0., zVtx, 0., 0., kFALSE, kFALSE);
1212 0 : }
1213 :
1214 : //__________________________________________________________________________
1215 : Double_t AliMFTTrackExtrap::TotalMomentumEnergyLoss(AliMFTTrackParam* trackParam, Double_t xVtx, Double_t yVtx, Double_t zVtx)
1216 : {
1217 : /// Calculate the total momentum energy loss in-between the track position and the vertex assuming a linear propagation
1218 :
1219 0 : if (trackParam->GetZ() == zVtx) return 0.; // nothing to be done if already at vertex
1220 :
1221 : // Check whether the geometry is available
1222 0 : if (!gGeoManager) {
1223 0 : cout<<"E-AliMFTTrackExtrap::TotalMomentumEnergyLoss: no TGeo"<<endl;
1224 0 : return 0.;
1225 : }
1226 :
1227 : // Get encountered material correction parameters assuming linear propagation from vertex to the track position
1228 0 : Double_t trackXYZOut[3];
1229 0 : trackXYZOut[0] = trackParam->GetX();
1230 0 : trackXYZOut[1] = trackParam->GetY();
1231 0 : trackXYZOut[2] = trackParam->GetZ();
1232 0 : Double_t trackXYZIn[3];
1233 0 : trackXYZIn[0] = xVtx;
1234 0 : trackXYZIn[1] = yVtx;
1235 0 : trackXYZIn[2] = zVtx;
1236 0 : Double_t pTot = trackParam->P();
1237 0 : Double_t pathLength, f0, f1, f2, meanRho, totalELoss, sigmaELoss2;
1238 0 : GetAbsorberCorrectionParam(trackXYZIn,trackXYZOut,pTot,pathLength,f0,f1,f2,meanRho,totalELoss,sigmaELoss2);
1239 :
1240 : // total momentum corrected for energy loss
1241 0 : Double_t muMass = TDatabasePDG::Instance()->GetParticle("mu-")->Mass(); // GeV
1242 0 : Double_t e = TMath::Sqrt(pTot*pTot + muMass*muMass);
1243 0 : Double_t eCorr = e + totalELoss;
1244 0 : Double_t pTotCorr = TMath::Sqrt(eCorr*eCorr - muMass*muMass);
1245 :
1246 0 : return pTotCorr - pTot;
1247 0 : }
1248 :
1249 : //__________________________________________________________________________
1250 : Double_t AliMFTTrackExtrap::BetheBloch(Double_t pTotal, Double_t pathLength, Double_t rho, Double_t atomicZ, Double_t atomicZoverA)
1251 : {
1252 : /// Returns the mean total momentum energy loss of muon with total momentum='pTotal'
1253 : /// in the absorber layer of lenght='pathLength', density='rho', A='atomicA' and Z='atomicZ'
1254 0 : Double_t muMass = TDatabasePDG::Instance()->GetParticle("mu-")->Mass(); // GeV
1255 :
1256 : // mean exitation energy (GeV)
1257 : Double_t i;
1258 0 : if (atomicZ < 13) i = (12. * atomicZ + 7.) * 1.e-9;
1259 0 : else i = (9.76 * atomicZ + 58.8 * TMath::Power(atomicZ,-0.19)) * 1.e-9;
1260 :
1261 0 : 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 : {
1267 : /// Returns the total momentum energy loss fluctuation of muon with total momentum='pTotal'
1268 : /// in the absorber layer of lenght='pathLength', density='rho', A='atomicA' and Z='atomicZ'
1269 0 : 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 0 : Double_t p2=pTotal*pTotal;
1273 0 : Double_t beta2=p2/(p2 + muMass*muMass);
1274 :
1275 0 : Double_t fwhm = 2. * k * rho * pathLength * atomicZoverA / beta2; // FWHM of the energy loss Landau distribution
1276 0 : 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 0 : return sigma;
1281 : }
1282 :
1283 : //__________________________________________________________________________
1284 : void AliMFTTrackExtrap::Cov2CovP(const TMatrixD ¶m, TMatrixD &cov)
1285 : {
1286 : /// change coordinate system: (X, SlopeX, Y, SlopeY, q/Pyz) -> (X, SlopeX, Y, SlopeY, q*PTot)
1287 : /// parameters (param) are given in the (X, SlopeX, Y, SlopeY, q/Pyz) coordinate system
1288 :
1289 : // charge * total momentum
1290 0 : Double_t qPTot = TMath::Sqrt(1. + param(1,0)*param(1,0) + param(3,0)*param(3,0)) /
1291 0 : TMath::Sqrt(1. + param(3,0)*param(3,0)) / param(4,0);
1292 :
1293 : // Jacobian of the opposite transformation
1294 0 : TMatrixD jacob(5,5);
1295 0 : jacob.UnitMatrix();
1296 0 : jacob(4,1) = qPTot * param(1,0) / (1. + param(1,0)*param(1,0) + param(3,0)*param(3,0));
1297 0 : jacob(4,3) = - qPTot * param(1,0) * param(1,0) * param(3,0) /
1298 0 : (1. + param(3,0)*param(3,0)) / (1. + param(1,0)*param(1,0) + param(3,0)*param(3,0));
1299 0 : jacob(4,4) = - qPTot / param(4,0);
1300 :
1301 : // compute covariances in new coordinate system
1302 0 : TMatrixD tmp(cov,TMatrixD::kMultTranspose,jacob);
1303 0 : cov.Mult(jacob,tmp);
1304 0 : }
1305 :
1306 : //__________________________________________________________________________
1307 : void AliMFTTrackExtrap::CovP2Cov(const TMatrixD ¶m, TMatrixD &covP)
1308 : {
1309 : /// change coordinate system: (X, SlopeX, Y, SlopeY, q*PTot) -> (X, SlopeX, Y, SlopeY, q/Pyz)
1310 : /// parameters (param) are given in the (X, SlopeX, Y, SlopeY, q/Pyz) coordinate system
1311 :
1312 : // charge * total momentum
1313 0 : Double_t qPTot = TMath::Sqrt(1. + param(1,0)*param(1,0) + param(3,0)*param(3,0)) /
1314 0 : TMath::Sqrt(1. + param(3,0)*param(3,0)) / param(4,0);
1315 :
1316 : // Jacobian of the transformation
1317 0 : TMatrixD jacob(5,5);
1318 0 : jacob.UnitMatrix();
1319 0 : jacob(4,1) = param(4,0) * param(1,0) / (1. + param(1,0)*param(1,0) + param(3,0)*param(3,0));
1320 0 : jacob(4,3) = - param(4,0) * param(1,0) * param(1,0) * param(3,0) /
1321 0 : (1. + param(3,0)*param(3,0)) / (1. + param(1,0)*param(1,0) + param(3,0)*param(3,0));
1322 0 : jacob(4,4) = - param(4,0) / qPTot;
1323 :
1324 : // compute covariances in new coordinate system
1325 0 : TMatrixD tmp(covP,TMatrixD::kMultTranspose,jacob);
1326 0 : covP.Mult(jacob,tmp);
1327 0 : }
1328 :
1329 : //__________________________________________________________________________
1330 : void AliMFTTrackExtrap::ExtrapOneStepHelix(Double_t charge, Double_t step, const Double_t *vect, Double_t *vout)
1331 : {
1332 : /// <pre>
1333 : /// ******************************************************************
1334 : /// * *
1335 : /// * Performs the tracking of one step in a magnetic field *
1336 : /// * The trajectory is assumed to be a helix in a constant field *
1337 : /// * taken at the mid point of the step. *
1338 : /// * Parameters: *
1339 : /// * input *
1340 : /// * STEP =arc length of the step asked *
1341 : /// * VECT =input vector (position,direction cos and momentum) *
1342 : /// * CHARGE= electric charge of the particle *
1343 : /// * output *
1344 : /// * VOUT = same as VECT after completion of the step *
1345 : /// * *
1346 : /// * ==>Called by : USER, GUSWIM *
1347 : /// * Author m.hansroul ********* *
1348 : /// * modified s.egli, s.v.levonian *
1349 : /// * modified v.perevoztchikov
1350 : /// * *
1351 : /// ******************************************************************
1352 : /// </pre>
1353 :
1354 : // modif: everything in double precision
1355 :
1356 0 : 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 0 : vout[kipp] = vect[kipp];
1377 0 : if (TMath::Abs(charge) < 0.00001) {
1378 0 : for (Int_t i = 0; i < 3; i++) {
1379 0 : vout[i] = vect[i] + step * vect[i+3];
1380 0 : vout[i+3] = vect[i+3];
1381 : }
1382 0 : return;
1383 : }
1384 0 : xyz[0] = vect[kix] + 0.5 * step * vect[kipx];
1385 0 : xyz[1] = vect[kiy] + 0.5 * step * vect[kipy];
1386 0 : xyz[2] = vect[kiz] + 0.5 * step * vect[kipz];
1387 :
1388 : //cmodif: call gufld (xyz, h) changed into:
1389 0 : TGeoGlobalMagField::Instance()->Field(xyz,h);
1390 0 : h2xy = h[0]*h[0] + h[1]*h[1];
1391 0 : h[3] = h[2]*h[2]+ h2xy;
1392 : // cout<<"Field ="<< h[0]<<" "<< h[1]<<" "<< h[2]<<" "<< h[3]<<endl;
1393 0 : if (h[3] < 1.e-12) {
1394 0 : for (Int_t i = 0; i < 3; i++) {
1395 0 : vout[i] = vect[i] + step * vect[i+3];
1396 0 : vout[i+3] = vect[i+3];
1397 : }
1398 0 : return;
1399 : }
1400 0 : if (h2xy < 1.e-12*h[3]) {
1401 0 : ExtrapOneStepHelix3(charge*h[2], step, vect, vout);
1402 0 : return;
1403 : }
1404 0 : h[3] = TMath::Sqrt(h[3]);
1405 0 : h[0] /= h[3];
1406 0 : h[1] /= h[3];
1407 0 : h[2] /= h[3];
1408 0 : h[3] *= kec;
1409 :
1410 0 : hxp[0] = h[1]*vect[kipz] - h[2]*vect[kipy];
1411 0 : hxp[1] = h[2]*vect[kipx] - h[0]*vect[kipz];
1412 0 : hxp[2] = h[0]*vect[kipy] - h[1]*vect[kipx];
1413 :
1414 0 : hp = h[0]*vect[kipx] + h[1]*vect[kipy] + h[2]*vect[kipz];
1415 :
1416 0 : rho = -charge*h[3]/vect[kipp];
1417 0 : tet = rho * step;
1418 :
1419 0 : if (TMath::Abs(tet) > 0.15) {
1420 0 : sint = TMath::Sin(tet);
1421 0 : sintt = (sint/tet);
1422 0 : tsint = (tet-sint)/tet;
1423 0 : cos1t = 2.*(TMath::Sin(0.5*tet))*(TMath::Sin(0.5*tet))/tet;
1424 0 : } else {
1425 0 : tsint = tet*tet/36.;
1426 0 : sintt = (1. - tsint);
1427 0 : sint = tet*sintt;
1428 0 : cos1t = 0.5*tet;
1429 : }
1430 :
1431 0 : f1 = step * sintt;
1432 0 : f2 = step * cos1t;
1433 0 : f3 = step * tsint * hp;
1434 0 : f4 = -tet*cos1t;
1435 : f5 = sint;
1436 0 : f6 = tet * cos1t * hp;
1437 :
1438 0 : vout[kix] = vect[kix] + f1*vect[kipx] + f2*hxp[0] + f3*h[0];
1439 0 : vout[kiy] = vect[kiy] + f1*vect[kipy] + f2*hxp[1] + f3*h[1];
1440 0 : vout[kiz] = vect[kiz] + f1*vect[kipz] + f2*hxp[2] + f3*h[2];
1441 :
1442 0 : vout[kipx] = vect[kipx] + f4*vect[kipx] + f5*hxp[0] + f6*h[0];
1443 0 : vout[kipy] = vect[kipy] + f4*vect[kipy] + f5*hxp[1] + f6*h[1];
1444 0 : 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 0 : return;
1450 0 : }
1451 :
1452 : //__________________________________________________________________________
1453 : void AliMFTTrackExtrap::ExtrapOneStepHelix3(Double_t field, Double_t step, const Double_t *vect, Double_t *vout)
1454 : {
1455 : /// <pre>
1456 : /// ******************************************************************
1457 : /// * *
1458 : /// * Tracking routine in a constant field oriented *
1459 : /// * along axis 3 *
1460 : /// * Tracking is performed with a conventional *
1461 : /// * helix step method *
1462 : /// * *
1463 : /// * ==>Called by : USER, GUSWIM *
1464 : /// * Authors R.Brun, M.Hansroul ********* *
1465 : /// * Rewritten V.Perevoztchikov
1466 : /// * *
1467 : /// ******************************************************************
1468 : /// </pre>
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 0 : vout[kipp] = vect[kipp];
1491 0 : h4 = field * kec;
1492 :
1493 0 : hxp[0] = - vect[kipy];
1494 0 : hxp[1] = + vect[kipx];
1495 :
1496 0 : hp = vect[kipz];
1497 :
1498 0 : rho = -h4/vect[kipp];
1499 0 : tet = rho * step;
1500 0 : if (TMath::Abs(tet) > 0.15) {
1501 0 : sint = TMath::Sin(tet);
1502 0 : sintt = (sint/tet);
1503 0 : tsint = (tet-sint)/tet;
1504 0 : cos1t = 2.* TMath::Sin(0.5*tet) * TMath::Sin(0.5*tet)/tet;
1505 0 : } else {
1506 0 : tsint = tet*tet/36.;
1507 0 : sintt = (1. - tsint);
1508 0 : sint = tet*sintt;
1509 0 : cos1t = 0.5*tet;
1510 : }
1511 :
1512 0 : f1 = step * sintt;
1513 0 : f2 = step * cos1t;
1514 0 : f3 = step * tsint * hp;
1515 0 : f4 = -tet*cos1t;
1516 : f5 = sint;
1517 0 : f6 = tet * cos1t * hp;
1518 :
1519 0 : vout[kix] = vect[kix] + f1*vect[kipx] + f2*hxp[0];
1520 0 : vout[kiy] = vect[kiy] + f1*vect[kipy] + f2*hxp[1];
1521 0 : vout[kiz] = vect[kiz] + f1*vect[kipz] + f3;
1522 :
1523 0 : vout[kipx] = vect[kipx] + f4*vect[kipx] + f5*hxp[0];
1524 0 : vout[kipy] = vect[kipy] + f4*vect[kipy] + f5*hxp[1];
1525 0 : vout[kipz] = vect[kipz] + f4*vect[kipz] + f6;
1526 :
1527 : return;
1528 0 : }
1529 :
1530 : //__________________________________________________________________________
1531 : Bool_t AliMFTTrackExtrap::ExtrapOneStepRungekutta(Double_t charge, Double_t step, const Double_t* vect, Double_t* vout)
1532 : {
1533 : /// <pre>
1534 : /// ******************************************************************
1535 : /// * *
1536 : /// * Runge-Kutta method for tracking a particle through a magnetic *
1537 : /// * field. Uses Nystroem algorithm (See Handbook Nat. Bur. of *
1538 : /// * Standards, procedure 25.5.20) *
1539 : /// * *
1540 : /// * Input parameters *
1541 : /// * CHARGE Particle charge *
1542 : /// * STEP Step size *
1543 : /// * VECT Initial co-ords,direction cosines,momentum *
1544 : /// * Output parameters *
1545 : /// * VOUT Output co-ords,direction cosines,momentum *
1546 : /// * User routine called *
1547 : /// * CALL GUFLD(X,F) *
1548 : /// * *
1549 : /// * ==>Called by : USER, GUSWIM *
1550 : /// * Authors R.Brun, M.Hansroul ********* *
1551 : /// * V.Perevoztchikov (CUT STEP implementation) *
1552 : /// * *
1553 : /// * *
1554 : /// ******************************************************************
1555 : /// </pre>
1556 :
1557 0 : Double_t h2, h4, f[4];
1558 0 : 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 0 : for(Int_t j = 0; j < 7; j++)
1598 0 : vout[j] = vect[j];
1599 :
1600 0 : Double_t pinv = kec * charge / vect[6];
1601 : Double_t tl = 0.;
1602 : Double_t h = step;
1603 : Double_t rest;
1604 :
1605 :
1606 0 : do {
1607 0 : rest = step - tl;
1608 0 : if (TMath::Abs(h) > TMath::Abs(rest)) h = rest;
1609 : //cmodif: call gufld(vout,f) changed into:
1610 0 : TGeoGlobalMagField::Instance()->Field(vout,f);
1611 :
1612 : // *
1613 : // * start of integration
1614 : // *
1615 0 : x = vout[0];
1616 0 : y = vout[1];
1617 0 : z = vout[2];
1618 0 : a = vout[3];
1619 0 : b = vout[4];
1620 0 : c = vout[5];
1621 :
1622 0 : h2 = khalf * h;
1623 0 : h4 = khalf * h2;
1624 0 : ph = pinv * h;
1625 0 : ph2 = khalf * ph;
1626 0 : secxs[0] = (b * f[2] - c * f[1]) * ph2;
1627 0 : secys[0] = (c * f[0] - a * f[2]) * ph2;
1628 0 : seczs[0] = (a * f[1] - b * f[0]) * ph2;
1629 0 : ang2 = (secxs[0]*secxs[0] + secys[0]*secys[0] + seczs[0]*seczs[0]);
1630 0 : if (ang2 > kpisqua) break;
1631 :
1632 0 : dxt = h2 * a + h4 * secxs[0];
1633 0 : dyt = h2 * b + h4 * secys[0];
1634 0 : dzt = h2 * c + h4 * seczs[0];
1635 0 : xt = x + dxt;
1636 0 : yt = y + dyt;
1637 0 : zt = z + dzt;
1638 : // *
1639 : // * second intermediate point
1640 : // *
1641 :
1642 0 : est = TMath::Abs(dxt) + TMath::Abs(dyt) + TMath::Abs(dzt);
1643 0 : if (est > h) {
1644 0 : if (ncut++ > maxcut) break;
1645 : h *= khalf;
1646 0 : continue;
1647 : }
1648 :
1649 0 : xyzt[0] = xt;
1650 0 : xyzt[1] = yt;
1651 0 : xyzt[2] = zt;
1652 :
1653 : //cmodif: call gufld(xyzt,f) changed into:
1654 0 : TGeoGlobalMagField::Instance()->Field(xyzt,f);
1655 :
1656 0 : at = a + secxs[0];
1657 0 : bt = b + secys[0];
1658 0 : ct = c + seczs[0];
1659 :
1660 0 : secxs[1] = (bt * f[2] - ct * f[1]) * ph2;
1661 0 : secys[1] = (ct * f[0] - at * f[2]) * ph2;
1662 0 : seczs[1] = (at * f[1] - bt * f[0]) * ph2;
1663 0 : at = a + secxs[1];
1664 0 : bt = b + secys[1];
1665 0 : ct = c + seczs[1];
1666 0 : secxs[2] = (bt * f[2] - ct * f[1]) * ph2;
1667 0 : secys[2] = (ct * f[0] - at * f[2]) * ph2;
1668 0 : seczs[2] = (at * f[1] - bt * f[0]) * ph2;
1669 0 : dxt = h * (a + secxs[2]);
1670 0 : dyt = h * (b + secys[2]);
1671 0 : dzt = h * (c + seczs[2]);
1672 0 : xt = x + dxt;
1673 0 : yt = y + dyt;
1674 0 : zt = z + dzt;
1675 0 : at = a + 2.*secxs[2];
1676 0 : bt = b + 2.*secys[2];
1677 0 : ct = c + 2.*seczs[2];
1678 :
1679 0 : est = TMath::Abs(dxt)+TMath::Abs(dyt)+TMath::Abs(dzt);
1680 0 : if (est > 2.*TMath::Abs(h)) {
1681 0 : if (ncut++ > maxcut) break;
1682 : h *= khalf;
1683 0 : continue;
1684 : }
1685 :
1686 0 : xyzt[0] = xt;
1687 0 : xyzt[1] = yt;
1688 0 : xyzt[2] = zt;
1689 :
1690 : //cmodif: call gufld(xyzt,f) changed into:
1691 0 : TGeoGlobalMagField::Instance()->Field(xyzt,f);
1692 :
1693 0 : z = z + (c + (seczs[0] + seczs[1] + seczs[2]) * kthird) * h;
1694 0 : y = y + (b + (secys[0] + secys[1] + secys[2]) * kthird) * h;
1695 0 : x = x + (a + (secxs[0] + secxs[1] + secxs[2]) * kthird) * h;
1696 :
1697 0 : secxs[3] = (bt*f[2] - ct*f[1])* ph2;
1698 0 : secys[3] = (ct*f[0] - at*f[2])* ph2;
1699 0 : seczs[3] = (at*f[1] - bt*f[0])* ph2;
1700 0 : a = a+(secxs[0]+secxs[3]+2. * (secxs[1]+secxs[2])) * kthird;
1701 0 : b = b+(secys[0]+secys[3]+2. * (secys[1]+secys[2])) * kthird;
1702 0 : c = c+(seczs[0]+seczs[3]+2. * (seczs[1]+seczs[2])) * kthird;
1703 :
1704 0 : est = TMath::Abs(secxs[0]+secxs[3] - (secxs[1]+secxs[2]))
1705 0 : + TMath::Abs(secys[0]+secys[3] - (secys[1]+secys[2]))
1706 0 : + TMath::Abs(seczs[0]+seczs[3] - (seczs[1]+seczs[2]));
1707 :
1708 0 : if (est > kdlt && TMath::Abs(h) > 1.e-4) {
1709 0 : if (ncut++ > maxcut) break;
1710 : h *= khalf;
1711 0 : continue;
1712 : }
1713 :
1714 : ncut = 0;
1715 : // * if too many iterations, go to helix
1716 0 : if (iter++ > maxit) break;
1717 :
1718 0 : tl += h;
1719 0 : if (est < kdlt32)
1720 0 : h *= 2.;
1721 0 : cba = 1./ TMath::Sqrt(a*a + b*b + c*c);
1722 0 : vout[0] = x;
1723 0 : vout[1] = y;
1724 0 : vout[2] = z;
1725 0 : vout[3] = cba*a;
1726 0 : vout[4] = cba*b;
1727 0 : vout[5] = cba*c;
1728 0 : rest = step - tl;
1729 0 : if (step < 0.) rest = -rest;
1730 0 : if (rest < 1.e-5*TMath::Abs(step)) return kTRUE;
1731 :
1732 : } while(1);
1733 :
1734 : // angle too big, use helix
1735 0 : cout<<"W-AliMFTTrackExtrap::ExtrapOneStepRungekutta: Ruge-Kutta failed: switch to helix"<<endl;
1736 :
1737 0 : f1 = f[0];
1738 0 : f2 = f[1];
1739 0 : f3 = f[2];
1740 0 : f4 = TMath::Sqrt(f1*f1+f2*f2+f3*f3);
1741 0 : if (f4 < 1.e-10) {
1742 0 : cout<<"E-AliMFTTrackExtrap::ExtrapOneStepRungekutta: magnetic field at (";
1743 0 : cout<<xyzt[0]<<", "<<xyzt[1]<<", "<<xyzt[2]<<") = "<<f4<<": giving up"<<endl;
1744 0 : return kFALSE;
1745 : }
1746 0 : rho = -f4*pinv;
1747 0 : tet = rho * step;
1748 :
1749 0 : hnorm = 1./f4;
1750 0 : f1 = f1*hnorm;
1751 0 : f2 = f2*hnorm;
1752 0 : f3 = f3*hnorm;
1753 :
1754 0 : hxp[0] = f2*vect[kipz] - f3*vect[kipy];
1755 0 : hxp[1] = f3*vect[kipx] - f1*vect[kipz];
1756 0 : hxp[2] = f1*vect[kipy] - f2*vect[kipx];
1757 :
1758 0 : hp = f1*vect[kipx] + f2*vect[kipy] + f3*vect[kipz];
1759 :
1760 0 : rho1 = 1./rho;
1761 0 : sint = TMath::Sin(tet);
1762 0 : cost = 2.*TMath::Sin(khalf*tet)*TMath::Sin(khalf*tet);
1763 :
1764 0 : g1 = sint*rho1;
1765 0 : g2 = cost*rho1;
1766 0 : g3 = (tet-sint) * hp*rho1;
1767 0 : g4 = -cost;
1768 : g5 = sint;
1769 0 : g6 = cost * hp;
1770 :
1771 0 : vout[kix] = vect[kix] + g1*vect[kipx] + g2*hxp[0] + g3*f1;
1772 0 : vout[kiy] = vect[kiy] + g1*vect[kipy] + g2*hxp[1] + g3*f2;
1773 0 : vout[kiz] = vect[kiz] + g1*vect[kipz] + g2*hxp[2] + g3*f3;
1774 :
1775 0 : vout[kipx] = vect[kipx] + g4*vect[kipx] + g5*hxp[0] + g6*f1;
1776 0 : vout[kipy] = vect[kipy] + g4*vect[kipy] + g5*hxp[1] + g6*f2;
1777 0 : vout[kipz] = vect[kipz] + g4*vect[kipz] + g5*hxp[2] + g6*f3;
1778 :
1779 0 : return kTRUE;
1780 0 : }
1781 :
|
