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 : // //
20 : // Implementation of the external track parameterisation class. //
21 : // //
22 : // This parameterisation is used to exchange tracks between the detectors. //
23 : // A set of functions returning the position and the momentum of tracks //
24 : // in the global coordinate system as well as the track impact parameters //
25 : // are implemented.
26 : // Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch //
27 : ///////////////////////////////////////////////////////////////////////////////
28 : #include <cassert>
29 :
30 : #include <TVectorD.h>
31 : #include <TMatrixDSym.h>
32 : #include <TPolyMarker3D.h>
33 : #include <TVector3.h>
34 : #include <TMatrixD.h>
35 :
36 : #include "AliExternalTrackParam.h"
37 : #include "AliVVertex.h"
38 : #include "AliLog.h"
39 :
40 176 : ClassImp(AliExternalTrackParam)
41 :
42 : Double32_t AliExternalTrackParam::fgMostProbablePt=kMostProbablePt;
43 : Bool_t AliExternalTrackParam::fgUseLogTermMS = kFALSE;;
44 : //_____________________________________________________________________________
45 : AliExternalTrackParam::AliExternalTrackParam() :
46 40576 : AliVTrack(),
47 40576 : fX(0),
48 40576 : fAlpha(0)
49 160388 : {
50 : //
51 : // default constructor
52 : //
53 486912 : for (Int_t i = 0; i < 5; i++) fP[i] = 0;
54 1298432 : for (Int_t i = 0; i < 15; i++) fC[i] = 0;
55 59906 : }
56 :
57 : //_____________________________________________________________________________
58 : AliExternalTrackParam::AliExternalTrackParam(const AliExternalTrackParam &track):
59 24232 : AliVTrack(track),
60 24232 : fX(track.fX),
61 24232 : fAlpha(track.fAlpha)
62 87370 : {
63 : //
64 : // copy constructor
65 : //
66 290784 : for (Int_t i = 0; i < 5; i++) fP[i] = track.fP[i];
67 775424 : for (Int_t i = 0; i < 15; i++) fC[i] = track.fC[i];
68 24232 : CheckCovariance();
69 31569 : }
70 :
71 : //_____________________________________________________________________________
72 : AliExternalTrackParam& AliExternalTrackParam::operator=(const AliExternalTrackParam &trkPar)
73 : {
74 : //
75 : // assignment operator
76 : //
77 :
78 31792 : if (this!=&trkPar) {
79 15896 : AliVTrack::operator=(trkPar);
80 15896 : fX = trkPar.fX;
81 15896 : fAlpha = trkPar.fAlpha;
82 :
83 190752 : for (Int_t i = 0; i < 5; i++) fP[i] = trkPar.fP[i];
84 508672 : for (Int_t i = 0; i < 15; i++) fC[i] = trkPar.fC[i];
85 15896 : CheckCovariance();
86 15896 : }
87 :
88 15896 : return *this;
89 : }
90 :
91 : //_____________________________________________________________________________
92 : AliExternalTrackParam::AliExternalTrackParam(Double_t x, Double_t alpha,
93 : const Double_t param[5],
94 : const Double_t covar[15]) :
95 96 : AliVTrack(),
96 96 : fX(x),
97 96 : fAlpha(alpha)
98 480 : {
99 : //
100 : // create external track parameters from given arguments
101 : //
102 1152 : for (Int_t i = 0; i < 5; i++) fP[i] = param[i];
103 3072 : for (Int_t i = 0; i < 15; i++) fC[i] = covar[i];
104 96 : CheckCovariance();
105 192 : }
106 :
107 : //_____________________________________________________________________________
108 : void AliExternalTrackParam::CopyFromVTrack(const AliVTrack *vTrack)
109 : {
110 : //
111 : // Recreate TrackParams from VTrack
112 : // This is not a copy contructor !
113 : //
114 0 : if (!vTrack) {
115 0 : AliError("Source VTrack is NULL");
116 0 : return;
117 : }
118 0 : if (this==vTrack) {
119 0 : AliError("Copy of itself is requested");
120 0 : return;
121 : }
122 : //
123 0 : if (vTrack->InheritsFrom(AliExternalTrackParam::Class())) {
124 0 : AliDebug(1,"Source itself is AliExternalTrackParam, using assignment operator");
125 0 : *this = *(AliExternalTrackParam*)vTrack;
126 0 : return;
127 : }
128 : //
129 0 : AliVTrack::operator=( *vTrack );
130 : //
131 0 : Double_t xyz[3],pxpypz[3],cv[21];
132 0 : vTrack->GetXYZ(xyz);
133 0 : pxpypz[0]=vTrack->Px();
134 0 : pxpypz[1]=vTrack->Py();
135 0 : pxpypz[2]=vTrack->Pz();
136 0 : vTrack->GetCovarianceXYZPxPyPz(cv);
137 0 : Short_t sign = (Short_t)vTrack->Charge();
138 0 : Set(xyz,pxpypz,cv,sign);
139 0 : }
140 :
141 : //_____________________________________________________________________________
142 : AliExternalTrackParam::AliExternalTrackParam(const AliVTrack *vTrack) :
143 0 : AliVTrack(),
144 0 : fX(0.),
145 0 : fAlpha(0.)
146 0 : {
147 : //
148 : // Constructor from virtual track,
149 : // This is not a copy contructor !
150 : //
151 :
152 0 : if (vTrack->InheritsFrom("AliExternalTrackParam")) {
153 0 : AliError("This is not a copy constructor. Use AliExternalTrackParam(const AliExternalTrackParam &) !");
154 0 : AliWarning("Calling the default constructor...");
155 0 : AliExternalTrackParam();
156 0 : return;
157 : }
158 :
159 0 : Double_t xyz[3],pxpypz[3],cv[21];
160 0 : vTrack->GetXYZ(xyz);
161 0 : pxpypz[0]=vTrack->Px();
162 0 : pxpypz[1]=vTrack->Py();
163 0 : pxpypz[2]=vTrack->Pz();
164 0 : vTrack->GetCovarianceXYZPxPyPz(cv);
165 0 : Short_t sign = (Short_t)vTrack->Charge();
166 :
167 0 : Set(xyz,pxpypz,cv,sign);
168 0 : }
169 :
170 : //_____________________________________________________________________________
171 : AliExternalTrackParam::AliExternalTrackParam(Double_t xyz[3],Double_t pxpypz[3],
172 : Double_t cv[21],Short_t sign) :
173 0 : AliVTrack(),
174 0 : fX(0.),
175 0 : fAlpha(0.)
176 0 : {
177 : //
178 : // constructor from the global parameters
179 : //
180 :
181 0 : Set(xyz,pxpypz,cv,sign);
182 0 : }
183 :
184 : /*
185 : //_____________________________________________________________________________
186 : void AliExternalTrackParam::Set(Double_t xyz[3],Double_t pxpypz[3],
187 : Double_t cv[21],Short_t sign)
188 : {
189 : //
190 : // create external track parameters from the global parameters
191 : // x,y,z,px,py,pz and their 6x6 covariance matrix
192 : // A.Dainese 10.10.08
193 :
194 : // Calculate alpha: the rotation angle of the corresponding local system.
195 : //
196 : // For global radial position inside the beam pipe, alpha is the
197 : // azimuthal angle of the momentum projected on (x,y).
198 : //
199 : // For global radial position outside the ITS, alpha is the
200 : // azimuthal angle of the centre of the TPC sector in which the point
201 : // xyz lies
202 : //
203 : const double kSafe = 1e-5;
204 : Double_t radPos2 = xyz[0]*xyz[0]+xyz[1]*xyz[1];
205 : Double_t radMax = 45.; // approximately ITS outer radius
206 : if (radPos2 < radMax*radMax) { // inside the ITS
207 : fAlpha = TMath::ATan2(pxpypz[1],pxpypz[0]);
208 : } else { // outside the ITS
209 : Float_t phiPos = TMath::Pi()+TMath::ATan2(-xyz[1], -xyz[0]);
210 : fAlpha =
211 : TMath::DegToRad()*(20*((((Int_t)(phiPos*TMath::RadToDeg()))/20))+10);
212 : }
213 : //
214 : Double_t cs=TMath::Cos(fAlpha), sn=TMath::Sin(fAlpha);
215 : // protection: avoid alpha being too close to 0 or +-pi/2
216 : if (TMath::Abs(sn)<2*kSafe) {
217 : if (fAlpha>0) fAlpha += fAlpha< TMath::Pi()/2. ? 2*kSafe : -2*kSafe;
218 : else fAlpha += fAlpha>-TMath::Pi()/2. ? -2*kSafe : 2*kSafe;
219 : cs=TMath::Cos(fAlpha);
220 : sn=TMath::Sin(fAlpha);
221 : }
222 : else if (TMath::Abs(cs)<2*kSafe) {
223 : if (fAlpha>0) fAlpha += fAlpha> TMath::Pi()/2. ? 2*kSafe : -2*kSafe;
224 : else fAlpha += fAlpha>-TMath::Pi()/2. ? 2*kSafe : -2*kSafe;
225 : cs=TMath::Cos(fAlpha);
226 : sn=TMath::Sin(fAlpha);
227 : }
228 : // Get the vertex of origin and the momentum
229 : TVector3 ver(xyz[0],xyz[1],xyz[2]);
230 : TVector3 mom(pxpypz[0],pxpypz[1],pxpypz[2]);
231 : //
232 : // avoid momenta along axis
233 : if (TMath::Abs(mom[0])<kSafe) mom[0] = TMath::Sign(kSafe*TMath::Abs(mom[1]), mom[0]);
234 : if (TMath::Abs(mom[1])<kSafe) mom[1] = TMath::Sign(kSafe*TMath::Abs(mom[0]), mom[1]);
235 :
236 : // Rotate to the local coordinate system
237 : ver.RotateZ(-fAlpha);
238 : mom.RotateZ(-fAlpha);
239 :
240 : //
241 : // x of the reference plane
242 : fX = ver.X();
243 :
244 : Double_t charge = (Double_t)sign;
245 :
246 : fP[0] = ver.Y();
247 : fP[1] = ver.Z();
248 : fP[2] = TMath::Sin(mom.Phi());
249 : fP[3] = mom.Pz()/mom.Pt();
250 : fP[4] = TMath::Sign(1/mom.Pt(),charge);
251 : //
252 : if (TMath::Abs( 1-fP[2]) < 3*kSafe) fP[2] = 1.- 3*kSafe; //Protection
253 : else if (TMath::Abs(-1-fP[2]) < 3*kSafe) fP[2] =-1.+ 3*kSafe; //Protection
254 : //
255 : // Covariance matrix (formulas to be simplified)
256 : Double_t pt=1./TMath::Abs(fP[4]);
257 : // avoid alpha+phi being to close to +-pi/2 in the cov.matrix evaluation
258 : double fp2 = fP[2];
259 : Double_t r=TMath::Sqrt((1.-fp2)*(1.+fp2));
260 : //
261 : Double_t m00=-sn;// m10=cs;
262 : Double_t m23=-pt*(sn + fp2*cs/r), m43=-pt*pt*(r*cs - fp2*sn);
263 : Double_t m24= pt*(cs - fp2*sn/r), m44=-pt*pt*(r*sn + fp2*cs);
264 : Double_t m35=pt, m45=-pt*pt*fP[3];
265 :
266 : m43*=GetSign();
267 : m44*=GetSign();
268 : m45*=GetSign();
269 :
270 : Double_t cv34 = TMath::Sqrt(cv[3 ]*cv[3 ]+cv[4 ]*cv[4 ]);
271 : Double_t a1=cv[13]-cv[9]*(m23*m44+m43*m24)/m23/m43;
272 : Double_t a2=m23*m24-m23*(m23*m44+m43*m24)/m43;
273 : Double_t a3=m43*m44-m43*(m23*m44+m43*m24)/m23;
274 : Double_t a4=cv[14]+2.*cv[9]; //cv[14]-2.*cv[9]*m24*m44/m23/m43;
275 : Double_t a5=m24*m24-2.*m24*m44*m23/m43;
276 : Double_t a6=m44*m44-2.*m24*m44*m43/m23;
277 :
278 : fC[0 ] = cv[0 ]+cv[2 ];
279 : fC[1 ] = TMath::Sign(cv34,cv[3 ]/m00);
280 : fC[2 ] = cv[5 ];
281 : fC[3 ] = (cv[10]*m43-cv[6]*m44)/(m24*m43-m23*m44)/m00;
282 : fC[10] = (cv[6]/m00-fC[3 ]*m23)/m43;
283 : fC[6 ] = (cv[15]/m00-fC[10]*m45)/m35;
284 : fC[4 ] = (cv[12]*m43-cv[8]*m44)/(m24*m43-m23*m44);
285 : fC[11] = (cv[8]-fC[4]*m23)/m43;
286 : fC[7 ] = cv[17]/m35-fC[11]*m45/m35;
287 : fC[5 ] = TMath::Abs((a4*a3-a6*a1)/(a5*a3-a6*a2));
288 : fC[14] = TMath::Abs((a1-a2*fC[5])/a3);
289 : fC[12] = (cv[9]-fC[5]*m23*m23-fC[14]*m43*m43)/m23/m43;
290 : Double_t b1=cv[18]-fC[12]*m23*m45-fC[14]*m43*m45;
291 : Double_t b2=m23*m35;
292 : Double_t b3=m43*m35;
293 : Double_t b4=cv[19]-fC[12]*m24*m45-fC[14]*m44*m45;
294 : Double_t b5=m24*m35;
295 : Double_t b6=m44*m35;
296 : fC[8 ] = (b4-b6*b1/b3)/(b5-b6*b2/b3);
297 : fC[13] = b1/b3-b2*fC[8]/b3;
298 : fC[9 ] = TMath::Abs((cv[20]-fC[14]*(m45*m45)-fC[13]*2.*m35*m45)/(m35*m35));
299 :
300 : CheckCovariance();
301 :
302 : return;
303 : }
304 : */
305 :
306 : //_____________________________________________________________________________
307 : void AliExternalTrackParam::Set(Double_t xyz[3],Double_t pxpypz[3],
308 : Double_t cv[21],Short_t sign)
309 : {
310 : //
311 : // create external track parameters from the global parameters
312 : // x,y,z,px,py,pz and their 6x6 covariance matrix
313 : // A.Dainese 10.10.08
314 :
315 : // Calculate alpha: the rotation angle of the corresponding local system.
316 : //
317 : // For global radial position inside the beam pipe, alpha is the
318 : // azimuthal angle of the momentum projected on (x,y).
319 : //
320 : // For global radial position outside the ITS, alpha is the
321 : // azimuthal angle of the centre of the TPC sector in which the point
322 : // xyz lies
323 : //
324 : const double kSafe = 1e-5;
325 8 : Double_t radPos2 = xyz[0]*xyz[0]+xyz[1]*xyz[1];
326 : Double_t radMax = 45.; // approximately ITS outer radius
327 4 : if (radPos2 < radMax*radMax) { // inside the ITS
328 0 : fAlpha = TMath::ATan2(pxpypz[1],pxpypz[0]);
329 0 : } else { // outside the ITS
330 4 : Float_t phiPos = TMath::Pi()+TMath::ATan2(-xyz[1], -xyz[0]);
331 4 : fAlpha =
332 4 : TMath::DegToRad()*(20*((((Int_t)(phiPos*TMath::RadToDeg()))/20))+10);
333 : }
334 : //
335 4 : Double_t cs=TMath::Cos(fAlpha), sn=TMath::Sin(fAlpha);
336 : // protection: avoid alpha being too close to 0 or +-pi/2
337 4 : if (TMath::Abs(sn)<2*kSafe) {
338 0 : if (fAlpha>0) fAlpha += fAlpha< TMath::Pi()/2. ? 2*kSafe : -2*kSafe;
339 0 : else fAlpha += fAlpha>-TMath::Pi()/2. ? -2*kSafe : 2*kSafe;
340 0 : cs=TMath::Cos(fAlpha);
341 0 : sn=TMath::Sin(fAlpha);
342 0 : }
343 4 : else if (TMath::Abs(cs)<2*kSafe) {
344 8 : if (fAlpha>0) fAlpha += fAlpha> TMath::Pi()/2. ? 2*kSafe : -2*kSafe;
345 0 : else fAlpha += fAlpha>-TMath::Pi()/2. ? 2*kSafe : -2*kSafe;
346 4 : cs=TMath::Cos(fAlpha);
347 4 : sn=TMath::Sin(fAlpha);
348 4 : }
349 : // Get the vertex of origin and the momentum
350 4 : TVector3 ver(xyz[0],xyz[1],xyz[2]);
351 4 : TVector3 mom(pxpypz[0],pxpypz[1],pxpypz[2]);
352 : //
353 : // Rotate to the local coordinate system
354 4 : ver.RotateZ(-fAlpha);
355 4 : mom.RotateZ(-fAlpha);
356 :
357 : //
358 : // x of the reference plane
359 4 : fX = ver.X();
360 :
361 4 : Double_t charge = (Double_t)sign;
362 :
363 4 : fP[0] = ver.Y();
364 4 : fP[1] = ver.Z();
365 8 : fP[2] = TMath::Sin(mom.Phi());
366 8 : fP[3] = mom.Pz()/mom.Pt();
367 8 : fP[4] = TMath::Sign(1/mom.Pt(),charge);
368 : //
369 4 : if (TMath::Abs( 1-fP[2]) < kSafe) fP[2] = 1.- kSafe; //Protection
370 4 : else if (TMath::Abs(-1-fP[2]) < kSafe) fP[2] =-1.+ kSafe; //Protection
371 : //
372 : // Covariance matrix (formulas to be simplified)
373 4 : Double_t pt=1./TMath::Abs(fP[4]);
374 4 : Double_t r=TMath::Sqrt((1.-fP[2])*(1.+fP[2]));
375 : //
376 4 : Double_t cv34 = TMath::Sqrt(cv[3 ]*cv[3 ]+cv[4 ]*cv[4 ]);
377 : //
378 : Int_t special = 0;
379 4 : double sgcheck = r*sn + fP[2]*cs;
380 4 : if (TMath::Abs(sgcheck)>=1-kSafe) { // special case: lab phi is +-pi/2
381 : special = 1;
382 0 : sgcheck = TMath::Sign(1.0,sgcheck);
383 0 : }
384 4 : else if (TMath::Abs(sgcheck)<kSafe) {
385 0 : sgcheck = TMath::Sign(1.0,cs);
386 : special = 2; // special case: lab phi is 0
387 0 : }
388 : //
389 4 : fC[0 ] = cv[0 ]+cv[2 ];
390 4 : fC[1 ] = TMath::Sign(cv34,-cv[3 ]*sn);
391 4 : fC[2 ] = cv[5 ];
392 : //
393 4 : if (special==1) {
394 0 : double pti = 1./pt;
395 0 : double pti2 = pti*pti;
396 0 : int q = GetSign();
397 0 : fC[3 ] = cv[6]*pti;
398 0 : fC[4 ] = -sgcheck*cv[8]*r*pti;
399 0 : fC[5 ] = TMath::Abs(cv[9]*r*r*pti2);
400 0 : fC[6 ] = (cv[10]*fP[3]-sgcheck*cv[15])*pti/r;
401 0 : fC[7 ] = (cv[17]-sgcheck*cv[12]*fP[3])*pti;
402 0 : fC[8 ] = (-sgcheck*cv[18]+cv[13]*fP[3])*r*pti2;
403 0 : fC[9 ] = TMath::Abs( cv[20]-2*sgcheck*cv[19]*fP[3]+cv[14]*fP[3]*fP[3])*pti2;
404 0 : fC[10] = cv[10]*pti2/r*q;
405 0 : fC[11] = -sgcheck*cv[12]*pti2*q;
406 0 : fC[12] = cv[13]*r*pti*pti2*q;
407 0 : fC[13] = (-sgcheck*cv[19]+cv[14]*fP[3])*r*pti2*pti;
408 0 : fC[14] = TMath::Abs(cv[14]*pti2*pti2);
409 4 : } else if (special==2) {
410 0 : double pti = 1./pt;
411 0 : double pti2 = pti*pti;
412 0 : int q = GetSign();
413 0 : fC[3 ] = -cv[10]*pti*cs/sn;
414 0 : fC[4 ] = cv[12]*cs*pti;
415 0 : fC[5 ] = TMath::Abs(cv[14]*cs*cs*pti2);
416 0 : fC[6 ] = (sgcheck*cv[6]*fP[3]-cv[15])*pti/sn;
417 0 : fC[7 ] = (cv[17]-sgcheck*cv[8]*fP[3])*pti;
418 0 : fC[8 ] = (cv[19]-sgcheck*cv[13]*fP[3])*cs*pti2;
419 0 : fC[9 ] = TMath::Abs( cv[20]-2*sgcheck*cv[18]*fP[3]+cv[9]*fP[3]*fP[3])*pti2;
420 0 : fC[10] = sgcheck*cv[6]*pti2/sn*q;
421 0 : fC[11] = -sgcheck*cv[8]*pti2*q;
422 0 : fC[12] = -sgcheck*cv[13]*cs*pti*pti2*q;
423 0 : fC[13] = (-sgcheck*cv[18]+cv[9]*fP[3])*pti2*pti*q;
424 0 : fC[14] = TMath::Abs(cv[9]*pti2*pti2);
425 0 : }
426 : else {
427 4 : Double_t m00=-sn;// m10=cs;
428 4 : Double_t m23=-pt*(sn + fP[2]*cs/r), m43=-pt*pt*(r*cs - fP[2]*sn);
429 4 : Double_t m24= pt*(cs - fP[2]*sn/r), m44=-pt*pt*(r*sn + fP[2]*cs);
430 4 : Double_t m35=pt, m45=-pt*pt*fP[3];
431 : //
432 8 : m43*=GetSign();
433 8 : m44*=GetSign();
434 8 : m45*=GetSign();
435 : //
436 4 : Double_t a1=cv[13]-cv[9]*(m23*m44+m43*m24)/m23/m43;
437 4 : Double_t a2=m23*m24-m23*(m23*m44+m43*m24)/m43;
438 4 : Double_t a3=m43*m44-m43*(m23*m44+m43*m24)/m23;
439 4 : Double_t a4=cv[14]+2.*cv[9]; //cv[14]-2.*cv[9]*m24*m44/m23/m43;
440 4 : Double_t a5=m24*m24-2.*m24*m44*m23/m43;
441 4 : Double_t a6=m44*m44-2.*m24*m44*m43/m23;
442 : //
443 4 : fC[3 ] = (cv[10]*m43-cv[6]*m44)/(m24*m43-m23*m44)/m00;
444 4 : fC[10] = (cv[6]/m00-fC[3 ]*m23)/m43;
445 4 : fC[6 ] = (cv[15]/m00-fC[10]*m45)/m35;
446 4 : fC[4 ] = (cv[12]*m43-cv[8]*m44)/(m24*m43-m23*m44);
447 4 : fC[11] = (cv[8]-fC[4]*m23)/m43;
448 4 : fC[7 ] = cv[17]/m35-fC[11]*m45/m35;
449 4 : fC[5 ] = TMath::Abs((a4*a3-a6*a1)/(a5*a3-a6*a2));
450 4 : fC[14] = TMath::Abs((a1-a2*fC[5])/a3);
451 4 : fC[12] = (cv[9]-fC[5]*m23*m23-fC[14]*m43*m43)/m23/m43;
452 4 : Double_t b1=cv[18]-fC[12]*m23*m45-fC[14]*m43*m45;
453 4 : Double_t b2=m23*m35;
454 4 : Double_t b3=m43*m35;
455 4 : Double_t b4=cv[19]-fC[12]*m24*m45-fC[14]*m44*m45;
456 4 : Double_t b5=m24*m35;
457 4 : Double_t b6=m44*m35;
458 4 : fC[8 ] = (b4-b6*b1/b3)/(b5-b6*b2/b3);
459 4 : fC[13] = b1/b3-b2*fC[8]/b3;
460 4 : fC[9 ] = TMath::Abs((cv[20]-fC[14]*(m45*m45)-fC[13]*2.*m35*m45)/(m35*m35));
461 : }
462 4 : CheckCovariance();
463 :
464 : return;
465 4 : }
466 :
467 : //_____________________________________________________________________________
468 : void AliExternalTrackParam::Reset() {
469 : //
470 : // Resets all the parameters to 0
471 : //
472 0 : fX=fAlpha=0.;
473 0 : for (Int_t i = 0; i < 5; i++) fP[i] = 0;
474 0 : for (Int_t i = 0; i < 15; i++) fC[i] = 0;
475 0 : }
476 :
477 : //_____________________________________________________________________________
478 : void AliExternalTrackParam::AddCovariance(const Double_t c[15]) {
479 : //
480 : // Add "something" to the track covarince matrix.
481 : // May be needed to account for unknown mis-calibration/mis-alignment
482 : //
483 24 : fC[0] +=c[0];
484 12 : fC[1] +=c[1]; fC[2] +=c[2];
485 12 : fC[3] +=c[3]; fC[4] +=c[4]; fC[5] +=c[5];
486 12 : fC[6] +=c[6]; fC[7] +=c[7]; fC[8] +=c[8]; fC[9] +=c[9];
487 12 : fC[10]+=c[10]; fC[11]+=c[11]; fC[12]+=c[12]; fC[13]+=c[13]; fC[14]+=c[14];
488 12 : CheckCovariance();
489 12 : }
490 :
491 :
492 : Double_t AliExternalTrackParam::GetP() const {
493 : //---------------------------------------------------------------------
494 : // This function returns the track momentum
495 : // Results for (nearly) straight tracks are meaningless !
496 : //---------------------------------------------------------------------
497 1437636 : if (TMath::Abs(fP[4])<=kAlmost0) return kVeryBig;
498 718818 : return TMath::Sqrt(1.+ fP[3]*fP[3])/TMath::Abs(fP[4]);
499 718818 : }
500 :
501 : Double_t AliExternalTrackParam::Get1P() const {
502 : //---------------------------------------------------------------------
503 : // This function returns the 1/(track momentum)
504 : //---------------------------------------------------------------------
505 0 : return TMath::Abs(fP[4])/TMath::Sqrt(1.+ fP[3]*fP[3]);
506 : }
507 :
508 : //_______________________________________________________________________
509 : Double_t AliExternalTrackParam::GetD(Double_t x,Double_t y,Double_t b) const {
510 : //------------------------------------------------------------------
511 : // This function calculates the transverse impact parameter
512 : // with respect to a point with global coordinates (x,y)
513 : // in the magnetic field "b" (kG)
514 : //------------------------------------------------------------------
515 101878 : if (TMath::Abs(b) < kAlmost0Field) return GetLinearD(x,y);
516 50939 : Double_t rp4=GetC(b);
517 :
518 50939 : Double_t xt=fX, yt=fP[0];
519 :
520 50939 : Double_t sn=TMath::Sin(fAlpha), cs=TMath::Cos(fAlpha);
521 50939 : Double_t a = x*cs + y*sn;
522 50939 : y = -x*sn + y*cs; x=a;
523 50939 : xt-=x; yt-=y;
524 :
525 50939 : sn=rp4*xt - fP[2]; cs=rp4*yt + TMath::Sqrt((1.- fP[2])*(1.+fP[2]));
526 50939 : a=2*(xt*fP[2] - yt*TMath::Sqrt((1.-fP[2])*(1.+fP[2])))-rp4*(xt*xt + yt*yt);
527 50939 : return -a/(1 + TMath::Sqrt(sn*sn + cs*cs));
528 50939 : }
529 :
530 : //_______________________________________________________________________
531 : void AliExternalTrackParam::
532 : GetDZ(Double_t x, Double_t y, Double_t z, Double_t b, Float_t dz[2]) const {
533 : //------------------------------------------------------------------
534 : // This function calculates the transverse and longitudinal impact parameters
535 : // with respect to a point with global coordinates (x,y)
536 : // in the magnetic field "b" (kG)
537 : //------------------------------------------------------------------
538 6356 : Double_t f1 = fP[2], r1 = TMath::Sqrt((1.-f1)*(1.+f1));
539 3178 : Double_t xt=fX, yt=fP[0];
540 3178 : Double_t sn=TMath::Sin(fAlpha), cs=TMath::Cos(fAlpha);
541 3178 : Double_t a = x*cs + y*sn;
542 3178 : y = -x*sn + y*cs; x=a;
543 3178 : xt-=x; yt-=y;
544 :
545 3178 : Double_t rp4=GetC(b);
546 6356 : if ((TMath::Abs(b) < kAlmost0Field) || (TMath::Abs(rp4) < kAlmost0)) {
547 0 : dz[0] = -(xt*f1 - yt*r1);
548 0 : dz[1] = fP[1] + (dz[0]*f1 - xt)/r1*fP[3] - z;
549 0 : return;
550 : }
551 :
552 3178 : sn=rp4*xt - f1; cs=rp4*yt + r1;
553 3178 : a=2*(xt*f1 - yt*r1)-rp4*(xt*xt + yt*yt);
554 3178 : Double_t rr=TMath::Sqrt(sn*sn + cs*cs);
555 3178 : dz[0] = -a/(1 + rr);
556 3178 : Double_t f2 = -sn/rr, r2 = TMath::Sqrt((1.-f2)*(1.+f2));
557 3178 : dz[1] = fP[1] + fP[3]/rp4*TMath::ASin(f2*r1 - f1*r2) - z;
558 6356 : }
559 :
560 : //_______________________________________________________________________
561 : Double_t AliExternalTrackParam::GetLinearD(Double_t xv,Double_t yv) const {
562 : //------------------------------------------------------------------
563 : // This function calculates the transverse impact parameter
564 : // with respect to a point with global coordinates (xv,yv)
565 : // neglecting the track curvature.
566 : //------------------------------------------------------------------
567 0 : Double_t sn=TMath::Sin(fAlpha), cs=TMath::Cos(fAlpha);
568 0 : Double_t x= xv*cs + yv*sn;
569 0 : Double_t y=-xv*sn + yv*cs;
570 :
571 0 : Double_t d = (fX-x)*fP[2] - (fP[0]-y)*TMath::Sqrt((1.-fP[2])*(1.+fP[2]));
572 :
573 0 : return -d;
574 : }
575 :
576 : Bool_t AliExternalTrackParam::CorrectForMeanMaterialdEdx
577 : (Double_t xOverX0, Double_t xTimesRho, Double_t mass,
578 : Double_t dEdx,
579 : Bool_t anglecorr) {
580 : //------------------------------------------------------------------
581 : // This function corrects the track parameters for the crossed material.
582 : // "xOverX0" - X/X0, the thickness in units of the radiation length.
583 : // "xTimesRho" - is the product length*density (g/cm^2).
584 : // It should be passed as negative when propagating tracks
585 : // from the intreaction point to the outside of the central barrel.
586 : // "mass" - the mass of this particle (GeV/c^2). Negative mass means charge=2 particle
587 : // "dEdx" - mean enery loss (GeV/(g/cm^2)
588 : // "anglecorr" - switch for the angular correction
589 : //------------------------------------------------------------------
590 202868 : Double_t &fP2=fP[2];
591 202868 : Double_t &fP3=fP[3];
592 202868 : Double_t &fP4=fP[4];
593 :
594 202868 : Double_t &fC22=fC[5];
595 202868 : Double_t &fC33=fC[9];
596 202868 : Double_t &fC43=fC[13];
597 202868 : Double_t &fC44=fC[14];
598 :
599 : //Apply angle correction, if requested
600 202868 : if(anglecorr) {
601 19360 : Double_t angle=TMath::Sqrt((1.+ fP3*fP3)/((1-fP2)*(1.+fP2)));
602 19360 : xOverX0 *=angle;
603 19360 : xTimesRho *=angle;
604 19360 : }
605 :
606 202868 : Double_t p=GetP();
607 202868 : if (mass<0) p += p; // q=2 particle
608 202868 : Double_t p2=p*p;
609 202868 : Double_t beta2=p2/(p2 + mass*mass);
610 :
611 : //Calculating the multiple scattering corrections******************
612 : Double_t cC22 = 0.;
613 : Double_t cC33 = 0.;
614 : Double_t cC43 = 0.;
615 : Double_t cC44 = 0.;
616 202868 : if (xOverX0 != 0) {
617 : //Double_t theta2=1.0259e-6*14*14/28/(beta2*p2)*TMath::Abs(d)*9.36*2.33;
618 200892 : Double_t theta2=0.0136*0.0136/(beta2*p2)*TMath::Abs(xOverX0);
619 200892 : if (GetUseLogTermMS()) {
620 0 : double lt = 1+0.038*TMath::Log(TMath::Abs(xOverX0));
621 0 : if (lt>0) theta2 *= lt*lt;
622 0 : }
623 200892 : if (mass<0) theta2 *= 4; // q=2 particle
624 200893 : if(theta2>TMath::Pi()*TMath::Pi()) return kFALSE;
625 200891 : cC22 = theta2*((1.-fP2)*(1.+fP2))*(1. + fP3*fP3);
626 200891 : cC33 = theta2*(1. + fP3*fP3)*(1. + fP3*fP3);
627 200891 : cC43 = theta2*fP3*fP4*(1. + fP3*fP3);
628 200891 : cC44 = theta2*fP3*fP4*fP3*fP4;
629 200891 : }
630 :
631 : //Calculating the energy loss corrections************************
632 : Double_t cP4=1.;
633 202867 : if ((xTimesRho != 0.) && (beta2 < 1.)) {
634 200891 : Double_t dE=dEdx*xTimesRho;
635 200891 : Double_t e=TMath::Sqrt(p2 + mass*mass);
636 200893 : if ( TMath::Abs(dE) > 0.3*e ) return kFALSE; //30% energy loss is too much!
637 200895 : if ( (1.+ dE/p2*(dE + 2*e)) < 0. ) return kFALSE;
638 200883 : cP4 = 1./TMath::Sqrt(1.+ dE/p2*(dE + 2*e)); //A precise formula by Ruben !
639 200883 : if (TMath::Abs(fP4*cP4)>100.) return kFALSE; //Do not track below 10 MeV/c
640 :
641 :
642 : // Approximate energy loss fluctuation (M.Ivanov)
643 : const Double_t knst=0.07; // To be tuned.
644 200883 : Double_t sigmadE=knst*TMath::Sqrt(TMath::Abs(dE));
645 200883 : cC44 += ((sigmadE*e/p2*fP4)*(sigmadE*e/p2*fP4));
646 :
647 200883 : }
648 :
649 : //Applying the corrections*****************************
650 202859 : fC22 += cC22;
651 202859 : fC33 += cC33;
652 202859 : fC43 += cC43;
653 202859 : fC44 += cC44;
654 202859 : fP4 *= cP4;
655 :
656 202859 : CheckCovariance();
657 :
658 202859 : return kTRUE;
659 202868 : }
660 :
661 : Bool_t AliExternalTrackParam::CorrectForMeanMaterial
662 : (Double_t xOverX0, Double_t xTimesRho, Double_t mass,
663 : Bool_t anglecorr,
664 : Double_t (*Bethe)(Double_t)) {
665 : //------------------------------------------------------------------
666 : // This function corrects the track parameters for the crossed material.
667 : // "xOverX0" - X/X0, the thickness in units of the radiation length.
668 : // "xTimesRho" - is the product length*density (g/cm^2).
669 : // It should be passed as negative when propagating tracks
670 : // from the intreaction point to the outside of the central barrel.
671 : // "mass" - the mass of this particle (GeV/c^2). mass<0 means charge=2
672 : // "anglecorr" - switch for the angular correction
673 : // "Bethe" - function calculating the energy loss (GeV/(g/cm^2))
674 : //------------------------------------------------------------------
675 :
676 202868 : Double_t bg=GetP()/mass;
677 202868 : if (mass<0) {
678 0 : if (mass<-990) {
679 0 : AliDebug(2,Form("Mass %f corresponds to unknown PID particle",mass));
680 0 : return kFALSE;
681 : }
682 0 : bg = -2*bg;
683 0 : }
684 202868 : Double_t dEdx=Bethe(bg);
685 202868 : if (mass<0) dEdx *= 4;
686 :
687 202868 : return CorrectForMeanMaterialdEdx(xOverX0,xTimesRho,mass,dEdx,anglecorr);
688 202868 : }
689 :
690 : Bool_t AliExternalTrackParam::CorrectForMeanMaterialZA
691 : (Double_t xOverX0, Double_t xTimesRho, Double_t mass,
692 : Double_t zOverA,
693 : Double_t density,
694 : Double_t exEnergy,
695 : Double_t jp1,
696 : Double_t jp2,
697 : Bool_t anglecorr) {
698 : //------------------------------------------------------------------
699 : // This function corrects the track parameters for the crossed material
700 : // using the full Geant-like Bethe-Bloch formula parameterization
701 : // "xOverX0" - X/X0, the thickness in units of the radiation length.
702 : // "xTimesRho" - is the product length*density (g/cm^2).
703 : // It should be passed as negative when propagating tracks
704 : // from the intreaction point to the outside of the central barrel.
705 : // "mass" - the mass of this particle (GeV/c^2). mass<0 means charge=2 particle
706 : // "density" - mean density (g/cm^3)
707 : // "zOverA" - mean Z/A
708 : // "exEnergy" - mean exitation energy (GeV)
709 : // "jp1" - density effect first junction point
710 : // "jp2" - density effect second junction point
711 : // "anglecorr" - switch for the angular correction
712 : //
713 : // The default values of the parameters are for silicon
714 : //
715 : //------------------------------------------------------------------
716 :
717 0 : Double_t bg=GetP()/mass;
718 0 : if (mass<0) {
719 0 : if (mass<-990) {
720 0 : AliDebug(2,Form("Mass %f corresponds to unknown PID particle",mass));
721 0 : return kFALSE;
722 : }
723 0 : bg = -2*bg;
724 0 : }
725 0 : Double_t dEdx=BetheBlochGeant(bg,density,jp1,jp2,exEnergy,zOverA);
726 :
727 0 : if (mass<0) dEdx *= 4;
728 0 : return CorrectForMeanMaterialdEdx(xOverX0,xTimesRho,mass,dEdx,anglecorr);
729 0 : }
730 :
731 :
732 :
733 : Bool_t AliExternalTrackParam::CorrectForMaterial
734 : (Double_t d, Double_t x0, Double_t mass, Double_t (*Bethe)(Double_t)) {
735 : //------------------------------------------------------------------
736 : // Deprecated function !
737 : // Better use CorrectForMeanMaterial instead of it.
738 : //
739 : // This function corrects the track parameters for the crossed material
740 : // "d" - the thickness (fraction of the radiation length)
741 : // It should be passed as negative when propagating tracks
742 : // from the intreaction point to the outside of the central barrel.
743 : // "x0" - the radiation length (g/cm^2)
744 : // "mass" - the mass of this particle (GeV/c^2)
745 : //------------------------------------------------------------------
746 :
747 0 : return CorrectForMeanMaterial(d,x0*d,mass,kTRUE,Bethe);
748 :
749 : }
750 :
751 : Double_t AliExternalTrackParam::BetheBlochAleph(Double_t bg,
752 : Double_t kp1,
753 : Double_t kp2,
754 : Double_t kp3,
755 : Double_t kp4,
756 : Double_t kp5) {
757 : //
758 : // This is the empirical ALEPH parameterization of the Bethe-Bloch formula.
759 : // It is normalized to 1 at the minimum.
760 : //
761 : // bg - beta*gamma
762 : //
763 : // The default values for the kp* parameters are for ALICE TPC.
764 : // The returned value is in MIP units
765 : //
766 :
767 1654632 : Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
768 :
769 827316 : Double_t aa = TMath::Power(beta,kp4);
770 827316 : Double_t bb = TMath::Power(1./bg,kp5);
771 :
772 827316 : bb=TMath::Log(kp3+bb);
773 :
774 827316 : return (kp2-aa-bb)*kp1/aa;
775 : }
776 :
777 : Double_t AliExternalTrackParam::BetheBlochGeant(Double_t bg,
778 : Double_t kp0,
779 : Double_t kp1,
780 : Double_t kp2,
781 : Double_t kp3,
782 : Double_t kp4) {
783 : //
784 : // This is the parameterization of the Bethe-Bloch formula inspired by Geant.
785 : //
786 : // bg - beta*gamma
787 : // kp0 - density [g/cm^3]
788 : // kp1 - density effect first junction point
789 : // kp2 - density effect second junction point
790 : // kp3 - mean excitation energy [GeV]
791 : // kp4 - mean Z/A
792 : //
793 : // The default values for the kp* parameters are for silicon.
794 : // The returned value is in [GeV/(g/cm^2)].
795 : //
796 :
797 : const Double_t mK = 0.307075e-3; // [GeV*cm^2/g]
798 : const Double_t me = 0.511e-3; // [GeV/c^2]
799 : const Double_t rho = kp0;
800 406304 : const Double_t x0 = kp1*2.303;
801 203152 : const Double_t x1 = kp2*2.303;
802 : const Double_t mI = kp3;
803 : const Double_t mZA = kp4;
804 203152 : const Double_t bg2 = bg*bg;
805 203152 : const Double_t maxT= 2*me*bg2; // neglecting the electron mass
806 :
807 : //*** Density effect
808 : Double_t d2=0.;
809 203152 : const Double_t x=TMath::Log(bg);
810 203152 : const Double_t lhwI=TMath::Log(28.816*1e-9*TMath::Sqrt(rho*mZA)/mI);
811 203152 : if (x > x1) {
812 829 : d2 = lhwI + x - 0.5;
813 203152 : } else if (x > x0) {
814 91217 : const Double_t r=(x1-x)/(x1-x0);
815 91217 : d2 = lhwI + x - 0.5 + (0.5 - lhwI - x0)*r*r*r;
816 91217 : }
817 :
818 406304 : return mK*mZA*(1+bg2)/bg2*
819 203152 : (0.5*TMath::Log(2*me*bg2*maxT/(mI*mI)) - bg2/(1+bg2) - d2);
820 : }
821 :
822 : Double_t AliExternalTrackParam::BetheBlochSolid(Double_t bg) {
823 : //------------------------------------------------------------------
824 : // This is an approximation of the Bethe-Bloch formula,
825 : // reasonable for solid materials.
826 : // All the parameters are, in fact, for Si.
827 : // The returned value is in [GeV/(g/cm^2)]
828 : //------------------------------------------------------------------
829 :
830 171368 : return BetheBlochGeant(bg);
831 : }
832 :
833 : Double_t AliExternalTrackParam::BetheBlochGas(Double_t bg) {
834 : //------------------------------------------------------------------
835 : // This is an approximation of the Bethe-Bloch formula,
836 : // reasonable for gas materials.
837 : // All the parameters are, in fact, for Ne.
838 : // The returned value is in [GeV/(g/cm^2)]
839 : //------------------------------------------------------------------
840 :
841 : const Double_t rho = 0.9e-3;
842 : const Double_t x0 = 2.;
843 : const Double_t x1 = 4.;
844 : const Double_t mI = 140.e-9;
845 : const Double_t mZA = 0.49555;
846 :
847 234368 : return BetheBlochGeant(bg,rho,x0,x1,mI,mZA);
848 : }
849 :
850 : Bool_t AliExternalTrackParam::Rotate(Double_t alpha) {
851 : //------------------------------------------------------------------
852 : // Transform this track to the local coord. system rotated
853 : // by angle "alpha" (rad) with respect to the global coord. system.
854 : //------------------------------------------------------------------
855 212690 : if (TMath::Abs(fP[2]) >= kAlmost1) {
856 0 : AliError(Form("Precondition is not satisfied: |sin(phi)|>1 ! %f",fP[2]));
857 0 : return kFALSE;
858 : }
859 :
860 106425 : if (alpha < -TMath::Pi()) alpha += 2*TMath::Pi();
861 114165 : else if (alpha >= TMath::Pi()) alpha -= 2*TMath::Pi();
862 :
863 106345 : Double_t &fP0=fP[0];
864 : Double_t &fP2=fP[2];
865 106345 : Double_t &fC00=fC[0];
866 106345 : Double_t &fC10=fC[1];
867 106345 : Double_t &fC20=fC[3];
868 106345 : Double_t &fC21=fC[4];
869 106345 : Double_t &fC22=fC[5];
870 106345 : Double_t &fC30=fC[6];
871 106345 : Double_t &fC32=fC[8];
872 106345 : Double_t &fC40=fC[10];
873 106345 : Double_t &fC42=fC[12];
874 :
875 106345 : Double_t x=fX;
876 106345 : Double_t ca=TMath::Cos(alpha-fAlpha), sa=TMath::Sin(alpha-fAlpha);
877 106345 : Double_t sf=fP2, cf=TMath::Sqrt((1.- fP2)*(1.+fP2)); // Improve precision
878 : // RS: check if rotation does no invalidate track model (cos(local_phi)>=0, i.e. particle
879 : // direction in local frame is along the X axis
880 106345 : if ((cf*ca+sf*sa)<0) {
881 480 : AliDebug(1,Form("Rotation failed: local cos(phi) would become %.2f",cf*ca+sf*sa));
882 160 : return kFALSE;
883 : }
884 : //
885 106185 : Double_t tmp=sf*ca - cf*sa;
886 :
887 106185 : if (TMath::Abs(tmp) >= kAlmost1) {
888 0 : if (TMath::Abs(tmp) > 1.+ Double_t(FLT_EPSILON))
889 0 : AliWarning(Form("Rotation failed ! %.10e",tmp));
890 0 : return kFALSE;
891 : }
892 106185 : fAlpha = alpha;
893 106185 : fX = x*ca + fP0*sa;
894 106185 : fP0= -x*sa + fP0*ca;
895 106185 : fP2= tmp;
896 :
897 106185 : if (TMath::Abs(cf)<kAlmost0) {
898 0 : AliError(Form("Too small cosine value %f",cf));
899 : cf = kAlmost0;
900 0 : }
901 :
902 106185 : Double_t rr=(ca+sf/cf*sa);
903 :
904 106185 : fC00 *= (ca*ca);
905 106185 : fC10 *= ca;
906 106185 : fC20 *= ca*rr;
907 106185 : fC21 *= rr;
908 106185 : fC22 *= rr*rr;
909 106185 : fC30 *= ca;
910 106185 : fC32 *= rr;
911 106185 : fC40 *= ca;
912 106185 : fC42 *= rr;
913 :
914 106185 : CheckCovariance();
915 :
916 : return kTRUE;
917 106345 : }
918 :
919 : //______________________________________________________
920 : Bool_t AliExternalTrackParam::RotateParamOnly(Double_t alpha)
921 : {
922 : // rotate to new frame, ignore covariance
923 168 : if (TMath::Abs(fP[2]) >= kAlmost1) {
924 0 : AliError(Form("Precondition is not satisfied: |sin(phi)|>1 ! %f",fP[2]));
925 0 : return kFALSE;
926 : }
927 : //
928 84 : if (alpha < -TMath::Pi()) alpha += 2*TMath::Pi();
929 168 : else if (alpha >= TMath::Pi()) alpha -= 2*TMath::Pi();
930 : //
931 84 : Double_t &fP0=fP[0];
932 : Double_t &fP2=fP[2];
933 : //
934 84 : Double_t x=fX;
935 84 : Double_t ca=TMath::Cos(alpha-fAlpha), sa=TMath::Sin(alpha-fAlpha);
936 84 : Double_t sf=fP2, cf=TMath::Sqrt((1.- fP2)*(1.+fP2)); // Improve precision
937 : // RS: check if rotation does no invalidate track model (cos(local_phi)>=0, i.e. particle
938 : // direction in local frame is along the X axis
939 84 : if ((cf*ca+sf*sa)<0) {
940 180 : AliDebug(1,Form("Rotation failed: local cos(phi) would become %.2f",cf*ca+sf*sa));
941 60 : return kFALSE;
942 : }
943 : //
944 24 : Double_t tmp=sf*ca - cf*sa;
945 :
946 24 : if (TMath::Abs(tmp) >= kAlmost1) {
947 0 : if (TMath::Abs(tmp) > 1.+ Double_t(FLT_EPSILON))
948 0 : AliWarning(Form("Rotation failed ! %.10e",tmp));
949 0 : return kFALSE;
950 : }
951 24 : fAlpha = alpha;
952 24 : fX = x*ca + fP0*sa;
953 24 : fP0= -x*sa + fP0*ca;
954 24 : fP2= tmp;
955 24 : return kTRUE;
956 84 : }
957 :
958 : Bool_t AliExternalTrackParam::Invert() {
959 : //------------------------------------------------------------------
960 : // Transform this track to the local coord. system rotated by 180 deg.
961 : //------------------------------------------------------------------
962 0 : fX = -fX;
963 0 : fAlpha += TMath::Pi();
964 0 : while (fAlpha < -TMath::Pi()) fAlpha += 2*TMath::Pi();
965 0 : while (fAlpha >= TMath::Pi()) fAlpha -= 2*TMath::Pi();
966 : //
967 0 : fP[0] = -fP[0];
968 : //fP[2] = -fP[2];
969 0 : fP[3] = -fP[3];
970 0 : fP[4] = -fP[4];
971 : //
972 0 : fC[1] = -fC[1]; // since the fP1 and fP2 are not inverted, their covariances with others change sign
973 0 : fC[3] = -fC[3];
974 0 : fC[7] = -fC[7];
975 0 : fC[8] = -fC[8];
976 0 : fC[11] = -fC[11];
977 0 : fC[12] = -fC[12];
978 : //
979 0 : return kTRUE;
980 : }
981 :
982 : Bool_t AliExternalTrackParam::PropagateTo(Double_t xk, Double_t b) {
983 : //----------------------------------------------------------------
984 : // Propagate this track to the plane X=xk (cm) in the field "b" (kG)
985 : //----------------------------------------------------------------
986 41482 : Double_t dx=xk-fX;
987 20867 : if (TMath::Abs(dx)<=kAlmost0) return kTRUE;
988 :
989 20615 : Double_t crv=GetC(b);
990 20615 : if (TMath::Abs(b) < kAlmost0Field) crv=0.;
991 :
992 20615 : Double_t x2r = crv*dx;
993 20615 : Double_t f1=fP[2], f2=f1 + x2r;
994 20615 : if (TMath::Abs(f1) >= kAlmost1) return kFALSE;
995 20616 : if (TMath::Abs(f2) >= kAlmost1) return kFALSE;
996 20614 : if (TMath::Abs(fP[4])< kAlmost0) return kFALSE;
997 :
998 20614 : Double_t &fP0=fP[0], &fP1=fP[1], &fP2=fP[2], &fP3=fP[3], &fP4=fP[4];
999 : Double_t
1000 20614 : &fC00=fC[0],
1001 20614 : &fC10=fC[1], &fC11=fC[2],
1002 20614 : &fC20=fC[3], &fC21=fC[4], &fC22=fC[5],
1003 20614 : &fC30=fC[6], &fC31=fC[7], &fC32=fC[8], &fC33=fC[9],
1004 20614 : &fC40=fC[10], &fC41=fC[11], &fC42=fC[12], &fC43=fC[13], &fC44=fC[14];
1005 :
1006 20614 : Double_t r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2));
1007 20614 : if (TMath::Abs(r1)<kAlmost0) return kFALSE;
1008 20614 : if (TMath::Abs(r2)<kAlmost0) return kFALSE;
1009 :
1010 20614 : fX=xk;
1011 20614 : double dy2dx = (f1+f2)/(r1+r2);
1012 20614 : fP0 += dx*dy2dx;
1013 20614 : fP2 += x2r;
1014 41134 : if (TMath::Abs(x2r)<0.05) fP1 += dx*(r2 + f2*dy2dx)*fP3; // Many thanks to P.Hristov !
1015 : else {
1016 : // for small dx/R the linear apporximation of the arc by the segment is OK,
1017 : // but at large dx/R the error is very large and leads to incorrect Z propagation
1018 : // angle traversed delta = 2*asin(dist_start_end / R / 2), hence the arc is: R*deltaPhi
1019 : // The dist_start_end is obtained from sqrt(dx^2+dy^2) = x/(r1+r2)*sqrt(2+f1*f2+r1*r2)
1020 : // double chord = dx*TMath::Sqrt(1+dy2dx*dy2dx); // distance from old position to new one
1021 : // double rot = 2*TMath::ASin(0.5*chord*crv); // angular difference seen from the circle center
1022 : // fP1 += rot/crv*fP3;
1023 : //
1024 94 : double rot = TMath::ASin(r1*f2 - r2*f1); // more economic version from Yura.
1025 98 : if (f1*f1+f2*f2>1 && f1*f2<0) { // special cases of large rotations or large abs angles
1026 0 : if (f2>0) rot = TMath::Pi() - rot; //
1027 0 : else rot = -TMath::Pi() - rot;
1028 : }
1029 94 : fP1 += fP3/crv*rot;
1030 : }
1031 :
1032 : //f = F - 1
1033 : /*
1034 : Double_t f02= dx/(r1*r1*r1); Double_t cc=crv/fP4;
1035 : Double_t f04=0.5*dx*dx/(r1*r1*r1); f04*=cc;
1036 : Double_t f12= dx*fP3*f1/(r1*r1*r1);
1037 : Double_t f14=0.5*dx*dx*fP3*f1/(r1*r1*r1); f14*=cc;
1038 : Double_t f13= dx/r1;
1039 : Double_t f24= dx; f24*=cc;
1040 : */
1041 20614 : Double_t rinv = 1./r1;
1042 20614 : Double_t r3inv = rinv*rinv*rinv;
1043 20614 : Double_t f24= x2r/fP4;
1044 20614 : Double_t f02= dx*r3inv;
1045 20614 : Double_t f04=0.5*f24*f02;
1046 20614 : Double_t f12= f02*fP3*f1;
1047 20614 : Double_t f14=0.5*f24*f02*fP3*f1;
1048 20614 : Double_t f13= dx*rinv;
1049 :
1050 : //b = C*ft
1051 20614 : Double_t b00=f02*fC20 + f04*fC40, b01=f12*fC20 + f14*fC40 + f13*fC30;
1052 20614 : Double_t b02=f24*fC40;
1053 20614 : Double_t b10=f02*fC21 + f04*fC41, b11=f12*fC21 + f14*fC41 + f13*fC31;
1054 20614 : Double_t b12=f24*fC41;
1055 20614 : Double_t b20=f02*fC22 + f04*fC42, b21=f12*fC22 + f14*fC42 + f13*fC32;
1056 20614 : Double_t b22=f24*fC42;
1057 20614 : Double_t b40=f02*fC42 + f04*fC44, b41=f12*fC42 + f14*fC44 + f13*fC43;
1058 20614 : Double_t b42=f24*fC44;
1059 20614 : Double_t b30=f02*fC32 + f04*fC43, b31=f12*fC32 + f14*fC43 + f13*fC33;
1060 20614 : Double_t b32=f24*fC43;
1061 :
1062 : //a = f*b = f*C*ft
1063 20614 : Double_t a00=f02*b20+f04*b40,a01=f02*b21+f04*b41,a02=f02*b22+f04*b42;
1064 20614 : Double_t a11=f12*b21+f14*b41+f13*b31,a12=f12*b22+f14*b42+f13*b32;
1065 20614 : Double_t a22=f24*b42;
1066 :
1067 : //F*C*Ft = C + (b + bt + a)
1068 20614 : fC00 += b00 + b00 + a00;
1069 20614 : fC10 += b10 + b01 + a01;
1070 20614 : fC20 += b20 + b02 + a02;
1071 20614 : fC30 += b30;
1072 20614 : fC40 += b40;
1073 20614 : fC11 += b11 + b11 + a11;
1074 20614 : fC21 += b21 + b12 + a12;
1075 20614 : fC31 += b31;
1076 20614 : fC41 += b41;
1077 20614 : fC22 += b22 + b22 + a22;
1078 20614 : fC32 += b32;
1079 20614 : fC42 += b42;
1080 :
1081 20614 : CheckCovariance();
1082 :
1083 : return kTRUE;
1084 20741 : }
1085 :
1086 : Bool_t AliExternalTrackParam::PropagateParamOnlyTo(Double_t xk, Double_t b) {
1087 : //----------------------------------------------------------------
1088 : // Propagate this track to the plane X=xk (cm) in the field "b" (kG)
1089 : // Only parameters are propagated, not the matrix. To be used for small
1090 : // distances only (<mm, i.e. misalignment)
1091 : //----------------------------------------------------------------
1092 0 : Double_t dx=xk-fX;
1093 0 : if (TMath::Abs(dx)<=kAlmost0) return kTRUE;
1094 :
1095 0 : Double_t crv=GetC(b);
1096 0 : if (TMath::Abs(b) < kAlmost0Field) crv=0.;
1097 :
1098 0 : Double_t x2r = crv*dx;
1099 0 : Double_t f1=fP[2], f2=f1 + x2r;
1100 0 : if (TMath::Abs(f1) >= kAlmost1) return kFALSE;
1101 0 : if (TMath::Abs(f2) >= kAlmost1) return kFALSE;
1102 0 : if (TMath::Abs(fP[4])< kAlmost0) return kFALSE;
1103 :
1104 0 : Double_t r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2));
1105 0 : if (TMath::Abs(r1)<kAlmost0) return kFALSE;
1106 0 : if (TMath::Abs(r2)<kAlmost0) return kFALSE;
1107 :
1108 0 : fX=xk;
1109 0 : double dy2dx = (f1+f2)/(r1+r2);
1110 0 : fP[0] += dx*dy2dx;
1111 0 : fP[2] += x2r;
1112 0 : if (TMath::Abs(x2r)<0.05) fP[1] += dx*(r2 + f2*dy2dx)*fP[3]; // Many thanks to P.Hristov !
1113 : else {
1114 : // for small dx/R the linear apporximation of the arc by the segment is OK,
1115 : // but at large dx/R the error is very large and leads to incorrect Z propagation
1116 : // angle traversed delta = 2*asin(dist_start_end / R / 2), hence the arc is: R*deltaPhi
1117 : // The dist_start_end is obtained from sqrt(dx^2+dy^2) = x/(r1+r2)*sqrt(2+f1*f2+r1*r2)
1118 : // double chord = dx*TMath::Sqrt(1+dy2dx*dy2dx); // distance from old position to new one
1119 : // double rot = 2*TMath::ASin(0.5*chord*crv); // angular difference seen from the circle center
1120 : // fP1 += rot/crv*fP3;
1121 : //
1122 0 : double rot = TMath::ASin(r1*f2 - r2*f1); // more economic version from Yura.
1123 0 : if (f1*f1+f2*f2>1 && f1*f2<0) { // special cases of large rotations or large abs angles
1124 0 : if (f2>0) rot = TMath::Pi() - rot; //
1125 0 : else rot = -TMath::Pi() - rot;
1126 : }
1127 0 : fP[1] += fP[3]/crv*rot;
1128 : }
1129 : return kTRUE;
1130 0 : }
1131 :
1132 : Bool_t
1133 : AliExternalTrackParam::Propagate(Double_t alpha, Double_t x, Double_t b) {
1134 : //------------------------------------------------------------------
1135 : // Transform this track to the local coord. system rotated
1136 : // by angle "alpha" (rad) with respect to the global coord. system,
1137 : // and propagate this track to the plane X=xk (cm) in the field "b" (kG)
1138 : //------------------------------------------------------------------
1139 :
1140 : //Save the parameters
1141 880 : Double_t as=fAlpha;
1142 440 : Double_t xs=fX;
1143 440 : Double_t ps[5], cs[15];
1144 5280 : for (Int_t i=0; i<5; i++) ps[i]=fP[i];
1145 14080 : for (Int_t i=0; i<15; i++) cs[i]=fC[i];
1146 :
1147 440 : if (Rotate(alpha))
1148 880 : if (PropagateTo(x,b)) return kTRUE;
1149 :
1150 : //Restore the parameters, if the operation failed
1151 0 : fAlpha=as;
1152 0 : fX=xs;
1153 0 : for (Int_t i=0; i<5; i++) fP[i]=ps[i];
1154 0 : for (Int_t i=0; i<15; i++) fC[i]=cs[i];
1155 0 : return kFALSE;
1156 440 : }
1157 :
1158 : Bool_t AliExternalTrackParam::PropagateBxByBz
1159 : (Double_t alpha, Double_t x, Double_t b[3]) {
1160 : //------------------------------------------------------------------
1161 : // Transform this track to the local coord. system rotated
1162 : // by angle "alpha" (rad) with respect to the global coord. system,
1163 : // and propagate this track to the plane X=xk (cm),
1164 : // taking into account all three components of the B field, "b[3]" (kG)
1165 : //------------------------------------------------------------------
1166 :
1167 : //Save the parameters
1168 193322 : Double_t as=fAlpha;
1169 96661 : Double_t xs=fX;
1170 96661 : Double_t ps[5], cs[15];
1171 1159932 : for (Int_t i=0; i<5; i++) ps[i]=fP[i];
1172 3093152 : for (Int_t i=0; i<15; i++) cs[i]=fC[i];
1173 :
1174 96661 : if (Rotate(alpha))
1175 193312 : if (PropagateToBxByBz(x,b)) return kTRUE;
1176 :
1177 : //Restore the parameters, if the operation failed
1178 6 : fAlpha=as;
1179 6 : fX=xs;
1180 72 : for (Int_t i=0; i<5; i++) fP[i]=ps[i];
1181 192 : for (Int_t i=0; i<15; i++) fC[i]=cs[i];
1182 6 : return kFALSE;
1183 96661 : }
1184 :
1185 :
1186 : void AliExternalTrackParam::Propagate(Double_t len, Double_t x[3],
1187 : Double_t p[3], Double_t bz) const {
1188 : //+++++++++++++++++++++++++++++++++++++++++
1189 : // Origin: K. Shileev (Kirill.Shileev@cern.ch)
1190 : // Extrapolate track along simple helix in magnetic field
1191 : // Arguments: len -distance alogn helix, [cm]
1192 : // bz - mag field, [kGaus]
1193 : // Returns: x and p contain extrapolated positon and momentum
1194 : // The momentum returned for straight-line tracks is meaningless !
1195 : //+++++++++++++++++++++++++++++++++++++++++
1196 0 : GetXYZ(x);
1197 :
1198 0 : if (OneOverPt() < kAlmost0 || TMath::Abs(bz) < kAlmost0Field || GetC(bz) < kAlmost0){ //straight-line tracks
1199 0 : Double_t unit[3]; GetDirection(unit);
1200 0 : x[0]+=unit[0]*len;
1201 0 : x[1]+=unit[1]*len;
1202 0 : x[2]+=unit[2]*len;
1203 :
1204 0 : p[0]=unit[0]/kAlmost0;
1205 0 : p[1]=unit[1]/kAlmost0;
1206 0 : p[2]=unit[2]/kAlmost0;
1207 0 : } else {
1208 0 : GetPxPyPz(p);
1209 0 : Double_t pp=GetP();
1210 0 : Double_t a = -kB2C*bz*GetSign();
1211 0 : Double_t rho = a/pp;
1212 0 : x[0] += p[0]*TMath::Sin(rho*len)/a - p[1]*(1-TMath::Cos(rho*len))/a;
1213 0 : x[1] += p[1]*TMath::Sin(rho*len)/a + p[0]*(1-TMath::Cos(rho*len))/a;
1214 0 : x[2] += p[2]*len/pp;
1215 :
1216 0 : Double_t p0=p[0];
1217 0 : p[0] = p0 *TMath::Cos(rho*len) - p[1]*TMath::Sin(rho*len);
1218 0 : p[1] = p[1]*TMath::Cos(rho*len) + p0 *TMath::Sin(rho*len);
1219 : }
1220 0 : }
1221 :
1222 : Bool_t AliExternalTrackParam::Intersect(Double_t pnt[3], Double_t norm[3],
1223 : Double_t bz) const {
1224 : //+++++++++++++++++++++++++++++++++++++++++
1225 : // Origin: K. Shileev (Kirill.Shileev@cern.ch)
1226 : // Finds point of intersection (if exists) of the helix with the plane.
1227 : // Stores result in fX and fP.
1228 : // Arguments: planePoint,planeNorm - the plane defined by any plane's point
1229 : // and vector, normal to the plane
1230 : // Returns: kTrue if helix intersects the plane, kFALSE otherwise.
1231 : //+++++++++++++++++++++++++++++++++++++++++
1232 0 : Double_t x0[3]; GetXYZ(x0); //get track position in MARS
1233 :
1234 : //estimates initial helix length up to plane
1235 : Double_t s=
1236 0 : (pnt[0]-x0[0])*norm[0] + (pnt[1]-x0[1])*norm[1] + (pnt[2]-x0[2])*norm[2];
1237 : Double_t dist=99999,distPrev=dist;
1238 0 : Double_t x[3],p[3];
1239 0 : while(TMath::Abs(dist)>0.00001){
1240 : //calculates helix at the distance s from x0 ALONG the helix
1241 0 : Propagate(s,x,p,bz);
1242 :
1243 : //distance between current helix position and plane
1244 0 : dist=(x[0]-pnt[0])*norm[0]+(x[1]-pnt[1])*norm[1]+(x[2]-pnt[2])*norm[2];
1245 :
1246 0 : if(TMath::Abs(dist) >= TMath::Abs(distPrev)) {return kFALSE;}
1247 : distPrev=dist;
1248 0 : s-=dist;
1249 : }
1250 : //on exit pnt is intersection point,norm is track vector at that point,
1251 : //all in MARS
1252 0 : for (Int_t i=0; i<3; i++) {pnt[i]=x[i]; norm[i]=p[i];}
1253 0 : return kTRUE;
1254 0 : }
1255 :
1256 : Double_t
1257 : AliExternalTrackParam::GetPredictedChi2(const Double_t p[2],const Double_t cov[3]) const {
1258 : //----------------------------------------------------------------
1259 : // Estimate the chi2 of the space point "p" with the cov. matrix "cov"
1260 : //----------------------------------------------------------------
1261 229880 : Double_t sdd = fC[0] + cov[0];
1262 114940 : Double_t sdz = fC[1] + cov[1];
1263 114940 : Double_t szz = fC[2] + cov[2];
1264 114940 : Double_t det = sdd*szz - sdz*sdz;
1265 :
1266 114940 : if (TMath::Abs(det) < kAlmost0) return kVeryBig;
1267 :
1268 114940 : Double_t d = fP[0] - p[0];
1269 114940 : Double_t z = fP[1] - p[1];
1270 :
1271 114940 : return (d*szz*d - 2*d*sdz*z + z*sdd*z)/det;
1272 114940 : }
1273 :
1274 : Double_t AliExternalTrackParam::
1275 : GetPredictedChi2(const Double_t p[3],const Double_t covyz[3],const Double_t covxyz[3]) const {
1276 : //----------------------------------------------------------------
1277 : // Estimate the chi2 of the 3D space point "p" and
1278 : // the full covariance matrix "covyz" and "covxyz"
1279 : //
1280 : // Cov(x,x) ... : covxyz[0]
1281 : // Cov(y,x) ... : covxyz[1] covyz[0]
1282 : // Cov(z,x) ... : covxyz[2] covyz[1] covyz[2]
1283 : //----------------------------------------------------------------
1284 :
1285 0 : Double_t res[3] = {
1286 0 : GetX() - p[0],
1287 0 : GetY() - p[1],
1288 0 : GetZ() - p[2]
1289 : };
1290 :
1291 0 : Double_t f=GetSnp();
1292 0 : if (TMath::Abs(f) >= kAlmost1) return kVeryBig;
1293 0 : Double_t r=TMath::Sqrt((1.-f)*(1.+f));
1294 0 : Double_t a=f/r, b=GetTgl()/r;
1295 :
1296 : Double_t s2=333.*333.; //something reasonably big (cm^2)
1297 :
1298 0 : TMatrixDSym v(3);
1299 0 : v(0,0)= s2; v(0,1)= a*s2; v(0,2)= b*s2;;
1300 0 : v(1,0)=a*s2; v(1,1)=a*a*s2 + GetSigmaY2(); v(1,2)=a*b*s2 + GetSigmaZY();
1301 0 : v(2,0)=b*s2; v(2,1)=a*b*s2 + GetSigmaZY(); v(2,2)=b*b*s2 + GetSigmaZ2();
1302 :
1303 0 : v(0,0)+=covxyz[0]; v(0,1)+=covxyz[1]; v(0,2)+=covxyz[2];
1304 0 : v(1,0)+=covxyz[1]; v(1,1)+=covyz[0]; v(1,2)+=covyz[1];
1305 0 : v(2,0)+=covxyz[2]; v(2,1)+=covyz[1]; v(2,2)+=covyz[2];
1306 :
1307 0 : v.Invert();
1308 0 : if (!v.IsValid()) return kVeryBig;
1309 :
1310 : Double_t chi2=0.;
1311 0 : for (Int_t i = 0; i < 3; i++)
1312 0 : for (Int_t j = 0; j < 3; j++) chi2 += res[i]*res[j]*v(i,j);
1313 :
1314 : return chi2;
1315 0 : }
1316 :
1317 : Double_t AliExternalTrackParam::
1318 : GetPredictedChi2(const AliExternalTrackParam *t) const {
1319 : //----------------------------------------------------------------
1320 : // Estimate the chi2 (5 dof) of this track with respect to the track
1321 : // given by the argument.
1322 : // The two tracks must be in the same reference system
1323 : // and estimated at the same reference plane.
1324 : //----------------------------------------------------------------
1325 :
1326 544 : if (TMath::Abs(t->GetAlpha()-GetAlpha()) > FLT_EPSILON) {
1327 0 : AliError("The reference systems of the tracks differ !");
1328 0 : return kVeryBig;
1329 : }
1330 272 : if (TMath::Abs(t->GetX()-GetX()) > FLT_EPSILON) {
1331 0 : AliError("The reference of the tracks planes differ !");
1332 0 : return kVeryBig;
1333 : }
1334 :
1335 272 : TMatrixDSym c(5);
1336 544 : c(0,0)=GetSigmaY2();
1337 816 : c(1,0)=GetSigmaZY(); c(1,1)=GetSigmaZ2();
1338 1360 : c(2,0)=GetSigmaSnpY(); c(2,1)=GetSigmaSnpZ(); c(2,2)=GetSigmaSnp2();
1339 1360 : c(3,0)=GetSigmaTglY(); c(3,1)=GetSigmaTglZ(); c(3,2)=GetSigmaTglSnp(); c(3,3)=GetSigmaTgl2();
1340 1632 : c(4,0)=GetSigma1PtY(); c(4,1)=GetSigma1PtZ(); c(4,2)=GetSigma1PtSnp(); c(4,3)=GetSigma1PtTgl(); c(4,4)=GetSigma1Pt2();
1341 :
1342 544 : c(0,0)+=t->GetSigmaY2();
1343 816 : c(1,0)+=t->GetSigmaZY(); c(1,1)+=t->GetSigmaZ2();
1344 1360 : c(2,0)+=t->GetSigmaSnpY();c(2,1)+=t->GetSigmaSnpZ();c(2,2)+=t->GetSigmaSnp2();
1345 1360 : c(3,0)+=t->GetSigmaTglY();c(3,1)+=t->GetSigmaTglZ();c(3,2)+=t->GetSigmaTglSnp();c(3,3)+=t->GetSigmaTgl2();
1346 1632 : c(4,0)+=t->GetSigma1PtY();c(4,1)+=t->GetSigma1PtZ();c(4,2)+=t->GetSigma1PtSnp();c(4,3)+=t->GetSigma1PtTgl();c(4,4)+=t->GetSigma1Pt2();
1347 816 : c(0,1)=c(1,0);
1348 1360 : c(0,2)=c(2,0); c(1,2)=c(2,1);
1349 1904 : c(0,3)=c(3,0); c(1,3)=c(3,1); c(2,3)=c(3,2);
1350 2448 : c(0,4)=c(4,0); c(1,4)=c(4,1); c(2,4)=c(4,2); c(3,4)=c(4,3);
1351 :
1352 272 : c.Invert();
1353 544 : if (!c.IsValid()) return kVeryBig;
1354 :
1355 :
1356 1632 : Double_t res[5] = {
1357 816 : GetY() - t->GetY(),
1358 816 : GetZ() - t->GetZ(),
1359 816 : GetSnp() - t->GetSnp(),
1360 816 : GetTgl() - t->GetTgl(),
1361 816 : GetSigned1Pt() - t->GetSigned1Pt()
1362 : };
1363 :
1364 : Double_t chi2=0.;
1365 3264 : for (Int_t i = 0; i < 5; i++)
1366 23120 : for (Int_t j = 0; j < 5; j++) chi2 += res[i]*res[j]*c(i,j);
1367 :
1368 : return chi2;
1369 816 : }
1370 :
1371 : Bool_t AliExternalTrackParam::
1372 : PropagateTo(Double_t p[3],Double_t covyz[3],Double_t covxyz[3],Double_t bz) {
1373 : //----------------------------------------------------------------
1374 : // Propagate this track to the plane
1375 : // the 3D space point "p" (with the covariance matrix "covyz" and "covxyz")
1376 : // belongs to.
1377 : // The magnetic field is "bz" (kG)
1378 : //
1379 : // The track curvature and the change of the covariance matrix
1380 : // of the track parameters are negleted !
1381 : // (So the "step" should be small compared with 1/curvature)
1382 : //----------------------------------------------------------------
1383 :
1384 0 : Double_t f=GetSnp();
1385 0 : if (TMath::Abs(f) >= kAlmost1) return kFALSE;
1386 0 : Double_t r=TMath::Sqrt((1.-f)*(1.+f));
1387 0 : Double_t a=f/r, b=GetTgl()/r;
1388 :
1389 : Double_t s2=333.*333.; //something reasonably big (cm^2)
1390 :
1391 0 : TMatrixDSym tV(3);
1392 0 : tV(0,0)= s2; tV(0,1)= a*s2; tV(0,2)= b*s2;
1393 0 : tV(1,0)=a*s2; tV(1,1)=a*a*s2; tV(1,2)=a*b*s2;
1394 0 : tV(2,0)=b*s2; tV(2,1)=a*b*s2; tV(2,2)=b*b*s2;
1395 :
1396 0 : TMatrixDSym pV(3);
1397 0 : pV(0,0)=covxyz[0]; pV(0,1)=covxyz[1]; pV(0,2)=covxyz[2];
1398 0 : pV(1,0)=covxyz[1]; pV(1,1)=covyz[0]; pV(1,2)=covyz[1];
1399 0 : pV(2,0)=covxyz[2]; pV(2,1)=covyz[1]; pV(2,2)=covyz[2];
1400 :
1401 0 : TMatrixDSym tpV(tV);
1402 0 : tpV+=pV;
1403 0 : tpV.Invert();
1404 0 : if (!tpV.IsValid()) return kFALSE;
1405 :
1406 0 : TMatrixDSym pW(3),tW(3);
1407 0 : for (Int_t i=0; i<3; i++)
1408 0 : for (Int_t j=0; j<3; j++) {
1409 0 : pW(i,j)=tW(i,j)=0.;
1410 0 : for (Int_t k=0; k<3; k++) {
1411 0 : pW(i,j) += tV(i,k)*tpV(k,j);
1412 0 : tW(i,j) += pV(i,k)*tpV(k,j);
1413 : }
1414 : }
1415 :
1416 0 : Double_t t[3] = {GetX(), GetY(), GetZ()};
1417 :
1418 : Double_t x=0.;
1419 0 : for (Int_t i=0; i<3; i++) x += (tW(0,i)*t[i] + pW(0,i)*p[i]);
1420 0 : Double_t crv=GetC(bz);
1421 0 : if (TMath::Abs(b) < kAlmost0Field) crv=0.;
1422 0 : f += crv*(x-fX);
1423 0 : if (TMath::Abs(f) >= kAlmost1) return kFALSE;
1424 0 : fX=x;
1425 :
1426 0 : fP[0]=0.;
1427 0 : for (Int_t i=0; i<3; i++) fP[0] += (tW(1,i)*t[i] + pW(1,i)*p[i]);
1428 0 : fP[1]=0.;
1429 0 : for (Int_t i=0; i<3; i++) fP[1] += (tW(2,i)*t[i] + pW(2,i)*p[i]);
1430 :
1431 0 : return kTRUE;
1432 0 : }
1433 :
1434 : Double_t *AliExternalTrackParam::GetResiduals(
1435 : Double_t *p,Double_t *cov,Bool_t updated) const {
1436 : //------------------------------------------------------------------
1437 : // Returns the track residuals with the space point "p" having
1438 : // the covariance matrix "cov".
1439 : // If "updated" is kTRUE, the track parameters expected to be updated,
1440 : // otherwise they must be predicted.
1441 : //------------------------------------------------------------------
1442 : static Double_t res[2];
1443 :
1444 0 : Double_t r00=cov[0], r01=cov[1], r11=cov[2];
1445 0 : if (updated) {
1446 0 : r00-=fC[0]; r01-=fC[1]; r11-=fC[2];
1447 0 : } else {
1448 0 : r00+=fC[0]; r01+=fC[1]; r11+=fC[2];
1449 : }
1450 0 : Double_t det=r00*r11 - r01*r01;
1451 :
1452 0 : if (TMath::Abs(det) < kAlmost0) return 0;
1453 :
1454 0 : Double_t tmp=r00; r00=r11/det; r11=tmp/det;
1455 :
1456 0 : if (r00 < 0.) return 0;
1457 0 : if (r11 < 0.) return 0;
1458 :
1459 0 : Double_t dy = fP[0] - p[0];
1460 0 : Double_t dz = fP[1] - p[1];
1461 :
1462 0 : res[0]=dy*TMath::Sqrt(r00);
1463 0 : res[1]=dz*TMath::Sqrt(r11);
1464 :
1465 : return res;
1466 0 : }
1467 :
1468 : Bool_t AliExternalTrackParam::Update(const Double_t p[2], const Double_t cov[3]) {
1469 : //------------------------------------------------------------------
1470 : // Update the track parameters with the space point "p" having
1471 : // the covariance matrix "cov"
1472 : //------------------------------------------------------------------
1473 229998 : Double_t &fP0=fP[0], &fP1=fP[1], &fP2=fP[2], &fP3=fP[3], &fP4=fP[4];
1474 : Double_t
1475 114999 : &fC00=fC[0],
1476 114999 : &fC10=fC[1], &fC11=fC[2],
1477 114999 : &fC20=fC[3], &fC21=fC[4], &fC22=fC[5],
1478 114999 : &fC30=fC[6], &fC31=fC[7], &fC32=fC[8], &fC33=fC[9],
1479 114999 : &fC40=fC[10], &fC41=fC[11], &fC42=fC[12], &fC43=fC[13], &fC44=fC[14];
1480 :
1481 114999 : Double_t r00=cov[0], r01=cov[1], r11=cov[2];
1482 114999 : r00+=fC00; r01+=fC10; r11+=fC11;
1483 114999 : Double_t det=r00*r11 - r01*r01;
1484 :
1485 114999 : if (TMath::Abs(det) < kAlmost0) return kFALSE;
1486 :
1487 :
1488 114999 : Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
1489 :
1490 114999 : Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11;
1491 114999 : Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11;
1492 114999 : Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11;
1493 114999 : Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11;
1494 114999 : Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11;
1495 :
1496 114999 : Double_t dy=p[0] - fP0, dz=p[1] - fP1;
1497 114999 : Double_t sf=fP2 + k20*dy + k21*dz;
1498 115011 : if (TMath::Abs(sf) > kAlmost1) return kFALSE;
1499 :
1500 114987 : fP0 += k00*dy + k01*dz;
1501 114987 : fP1 += k10*dy + k11*dz;
1502 114987 : fP2 = sf;
1503 114987 : fP3 += k30*dy + k31*dz;
1504 114987 : fP4 += k40*dy + k41*dz;
1505 :
1506 114987 : Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
1507 114987 : Double_t c12=fC21, c13=fC31, c14=fC41;
1508 :
1509 114987 : fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
1510 114987 : fC20-=k00*c02+k01*c12; fC30-=k00*c03+k01*c13;
1511 114987 : fC40-=k00*c04+k01*c14;
1512 :
1513 114987 : fC11-=k10*c01+k11*fC11;
1514 114987 : fC21-=k10*c02+k11*c12; fC31-=k10*c03+k11*c13;
1515 114987 : fC41-=k10*c04+k11*c14;
1516 :
1517 114987 : fC22-=k20*c02+k21*c12; fC32-=k20*c03+k21*c13;
1518 114987 : fC42-=k20*c04+k21*c14;
1519 :
1520 114987 : fC33-=k30*c03+k31*c13;
1521 114987 : fC43-=k30*c04+k31*c14;
1522 :
1523 114987 : fC44-=k40*c04+k41*c14;
1524 :
1525 114987 : CheckCovariance();
1526 :
1527 : return kTRUE;
1528 114999 : }
1529 :
1530 : void
1531 : AliExternalTrackParam::GetHelixParameters(Double_t hlx[6], Double_t b) const {
1532 : //--------------------------------------------------------------------
1533 : // External track parameters -> helix parameters
1534 : // "b" - magnetic field (kG)
1535 : //--------------------------------------------------------------------
1536 656 : Double_t cs=TMath::Cos(fAlpha), sn=TMath::Sin(fAlpha);
1537 :
1538 328 : hlx[0]=fP[0]; hlx[1]=fP[1]; hlx[2]=fP[2]; hlx[3]=fP[3];
1539 :
1540 328 : hlx[5]=fX*cs - hlx[0]*sn; // x0
1541 328 : hlx[0]=fX*sn + hlx[0]*cs; // y0
1542 : //hlx[1]= // z0
1543 328 : hlx[2]=TMath::ASin(hlx[2]) + fAlpha; // phi0
1544 : //hlx[3]= // tgl
1545 328 : hlx[4]=GetC(b); // C
1546 328 : }
1547 :
1548 :
1549 : static void Evaluate(const Double_t *h, Double_t t,
1550 : Double_t r[3], //radius vector
1551 : Double_t g[3], //first defivatives
1552 : Double_t gg[3]) //second derivatives
1553 : {
1554 : //--------------------------------------------------------------------
1555 : // Calculate position of a point on a track and some derivatives
1556 : //--------------------------------------------------------------------
1557 3312 : Double_t phase=h[4]*t+h[2];
1558 1656 : Double_t sn=TMath::Sin(phase), cs=TMath::Cos(phase);
1559 :
1560 1656 : r[0] = h[5];
1561 1656 : r[1] = h[0];
1562 1656 : if (TMath::Abs(h[4])>kAlmost0) {
1563 1656 : r[0] += (sn - h[6])/h[4];
1564 1656 : r[1] -= (cs - h[7])/h[4];
1565 1656 : } else {
1566 0 : r[0] += t*cs;
1567 0 : r[1] -= -t*sn;
1568 : }
1569 1656 : r[2] = h[1] + h[3]*t;
1570 :
1571 1656 : g[0] = cs; g[1]=sn; g[2]=h[3];
1572 :
1573 1656 : gg[0]=-h[4]*sn; gg[1]=h[4]*cs; gg[2]=0.;
1574 1656 : }
1575 :
1576 : Double_t AliExternalTrackParam::GetDCA(const AliExternalTrackParam *p,
1577 : Double_t b, Double_t &xthis, Double_t &xp) const {
1578 : //------------------------------------------------------------
1579 : // Returns the (weighed !) distance of closest approach between
1580 : // this track and the track "p".
1581 : // Other returned values:
1582 : // xthis, xt - coordinates of tracks' reference planes at the DCA
1583 : //-----------------------------------------------------------
1584 328 : Double_t dy2=GetSigmaY2() + p->GetSigmaY2();
1585 164 : Double_t dz2=GetSigmaZ2() + p->GetSigmaZ2();
1586 : Double_t dx2=dy2;
1587 :
1588 164 : Double_t p1[8]; GetHelixParameters(p1,b);
1589 164 : p1[6]=TMath::Sin(p1[2]); p1[7]=TMath::Cos(p1[2]);
1590 164 : Double_t p2[8]; p->GetHelixParameters(p2,b);
1591 164 : p2[6]=TMath::Sin(p2[2]); p2[7]=TMath::Cos(p2[2]);
1592 :
1593 :
1594 164 : Double_t r1[3],g1[3],gg1[3]; Double_t t1=0.;
1595 164 : Evaluate(p1,t1,r1,g1,gg1);
1596 164 : Double_t r2[3],g2[3],gg2[3]; Double_t t2=0.;
1597 164 : Evaluate(p2,t2,r2,g2,gg2);
1598 :
1599 164 : Double_t dx=r2[0]-r1[0], dy=r2[1]-r1[1], dz=r2[2]-r1[2];
1600 164 : Double_t dm=dx*dx/dx2 + dy*dy/dy2 + dz*dz/dz2;
1601 :
1602 : Int_t max=27;
1603 954 : while (max--) {
1604 790 : Double_t gt1=-(dx*g1[0]/dx2 + dy*g1[1]/dy2 + dz*g1[2]/dz2);
1605 790 : Double_t gt2=+(dx*g2[0]/dx2 + dy*g2[1]/dy2 + dz*g2[2]/dz2);
1606 1580 : Double_t h11=(g1[0]*g1[0] - dx*gg1[0])/dx2 +
1607 1580 : (g1[1]*g1[1] - dy*gg1[1])/dy2 +
1608 790 : (g1[2]*g1[2] - dz*gg1[2])/dz2;
1609 1580 : Double_t h22=(g2[0]*g2[0] + dx*gg2[0])/dx2 +
1610 1580 : (g2[1]*g2[1] + dy*gg2[1])/dy2 +
1611 790 : (g2[2]*g2[2] + dz*gg2[2])/dz2;
1612 790 : Double_t h12=-(g1[0]*g2[0]/dx2 + g1[1]*g2[1]/dy2 + g1[2]*g2[2]/dz2);
1613 :
1614 790 : Double_t det=h11*h22-h12*h12;
1615 :
1616 : Double_t dt1,dt2;
1617 790 : if (TMath::Abs(det)<1.e-33) {
1618 : //(quasi)singular Hessian
1619 0 : dt1=-gt1; dt2=-gt2;
1620 0 : } else {
1621 790 : dt1=-(gt1*h22 - gt2*h12)/det;
1622 790 : dt2=-(h11*gt2 - h12*gt1)/det;
1623 : }
1624 :
1625 819 : if ((dt1*gt1+dt2*gt2)>0) {dt1=-dt1; dt2=-dt2;}
1626 :
1627 : //check delta(phase1) ?
1628 : //check delta(phase2) ?
1629 :
1630 790 : if (TMath::Abs(dt1)/(TMath::Abs(t1)+1.e-3) < 1.e-4)
1631 166 : if (TMath::Abs(dt2)/(TMath::Abs(t2)+1.e-3) < 1.e-4) {
1632 164 : if ((gt1*gt1+gt2*gt2) > 1.e-4/dy2/dy2)
1633 0 : AliDebug(1," stopped at not a stationary point !");
1634 164 : Double_t lmb=h11+h22; lmb=lmb-TMath::Sqrt(lmb*lmb-4*det);
1635 164 : if (lmb < 0.)
1636 0 : AliDebug(1," stopped at not a minimum !");
1637 : break;
1638 : }
1639 :
1640 : Double_t dd=dm;
1641 664 : for (Int_t div=1 ; ; div*=2) {
1642 664 : Evaluate(p1,t1+dt1,r1,g1,gg1);
1643 664 : Evaluate(p2,t2+dt2,r2,g2,gg2);
1644 664 : dx=r2[0]-r1[0]; dy=r2[1]-r1[1]; dz=r2[2]-r1[2];
1645 664 : dd=dx*dx/dx2 + dy*dy/dy2 + dz*dz/dz2;
1646 1290 : if (dd<dm) break;
1647 38 : dt1*=0.5; dt2*=0.5;
1648 38 : if (div>512) {
1649 0 : AliDebug(1," overshoot !"); break;
1650 : }
1651 : }
1652 : dm=dd;
1653 :
1654 626 : t1+=dt1;
1655 626 : t2+=dt2;
1656 :
1657 626 : }
1658 :
1659 164 : if (max<=0) AliDebug(1," too many iterations !");
1660 :
1661 164 : Double_t cs=TMath::Cos(GetAlpha());
1662 164 : Double_t sn=TMath::Sin(GetAlpha());
1663 164 : xthis=r1[0]*cs + r1[1]*sn;
1664 :
1665 164 : cs=TMath::Cos(p->GetAlpha());
1666 164 : sn=TMath::Sin(p->GetAlpha());
1667 164 : xp=r2[0]*cs + r2[1]*sn;
1668 :
1669 328 : return TMath::Sqrt(dm*TMath::Sqrt(dy2*dz2));
1670 164 : }
1671 :
1672 : Double_t AliExternalTrackParam::
1673 : PropagateToDCA(AliExternalTrackParam *p, Double_t b) {
1674 : //--------------------------------------------------------------
1675 : // Propagates this track and the argument track to the position of the
1676 : // distance of closest approach.
1677 : // Returns the (weighed !) distance of closest approach.
1678 : //--------------------------------------------------------------
1679 0 : Double_t xthis,xp;
1680 0 : Double_t dca=GetDCA(p,b,xthis,xp);
1681 :
1682 0 : if (!PropagateTo(xthis,b)) {
1683 : //AliWarning(" propagation failed !");
1684 0 : return 1e+33;
1685 : }
1686 :
1687 0 : if (!p->PropagateTo(xp,b)) {
1688 : //AliWarning(" propagation failed !";
1689 0 : return 1e+33;
1690 : }
1691 :
1692 0 : return dca;
1693 0 : }
1694 :
1695 :
1696 : Bool_t AliExternalTrackParam::PropagateToDCA(const AliVVertex *vtx,
1697 : Double_t b, Double_t maxd, Double_t dz[2], Double_t covar[3]) {
1698 : //
1699 : // Propagate this track to the DCA to vertex "vtx",
1700 : // if the (rough) transverse impact parameter is not bigger then "maxd".
1701 : // Magnetic field is "b" (kG).
1702 : //
1703 : // a) The track gets extapolated to the DCA to the vertex.
1704 : // b) The impact parameters and their covariance matrix are calculated.
1705 : //
1706 : // In the case of success, the returned value is kTRUE
1707 : // (otherwise, it's kFALSE)
1708 : //
1709 880 : Double_t alpha=GetAlpha();
1710 440 : Double_t sn=TMath::Sin(alpha), cs=TMath::Cos(alpha);
1711 440 : Double_t x=GetX(), y=GetParameter()[0], snp=GetParameter()[2];
1712 440 : Double_t xv= vtx->GetX()*cs + vtx->GetY()*sn;
1713 440 : Double_t yv=-vtx->GetX()*sn + vtx->GetY()*cs, zv=vtx->GetZ();
1714 440 : x-=xv; y-=yv;
1715 :
1716 : //Estimate the impact parameter neglecting the track curvature
1717 440 : Double_t d=TMath::Abs(x*snp - y*TMath::Sqrt((1.-snp)*(1.+snp)));
1718 440 : if (d > maxd) return kFALSE;
1719 :
1720 : //Propagate to the DCA
1721 440 : Double_t crv=GetC(b);
1722 440 : if (TMath::Abs(b) < kAlmost0Field) crv=0.;
1723 :
1724 440 : Double_t tgfv=-(crv*x - snp)/(crv*y + TMath::Sqrt((1.-snp)*(1.+snp)));
1725 440 : sn=tgfv/TMath::Sqrt(1.+ tgfv*tgfv); cs=TMath::Sqrt((1.-sn)*(1.+sn));
1726 880 : if (TMath::Abs(tgfv)>0.) cs = sn/tgfv;
1727 : else cs=1.;
1728 :
1729 440 : x = xv*cs + yv*sn;
1730 440 : yv=-xv*sn + yv*cs; xv=x;
1731 :
1732 440 : if (!Propagate(alpha+TMath::ASin(sn),xv,b)) return kFALSE;
1733 :
1734 440 : if (dz==0) return kTRUE;
1735 440 : dz[0] = GetParameter()[0] - yv;
1736 440 : dz[1] = GetParameter()[1] - zv;
1737 :
1738 440 : if (covar==0) return kTRUE;
1739 440 : Double_t cov[6]; vtx->GetCovarianceMatrix(cov);
1740 :
1741 : //***** Improvements by A.Dainese
1742 440 : alpha=GetAlpha(); sn=TMath::Sin(alpha); cs=TMath::Cos(alpha);
1743 440 : Double_t s2ylocvtx = cov[0]*sn*sn + cov[2]*cs*cs - 2.*cov[1]*cs*sn;
1744 440 : covar[0] = GetCovariance()[0] + s2ylocvtx; // neglecting correlations
1745 440 : covar[1] = GetCovariance()[1]; // between (x,y) and z
1746 440 : covar[2] = GetCovariance()[2] + cov[5]; // in vertex's covariance matrix
1747 : //*****
1748 :
1749 : return kTRUE;
1750 880 : }
1751 :
1752 : Bool_t AliExternalTrackParam::PropagateToDCABxByBz(const AliVVertex *vtx,
1753 : Double_t b[3], Double_t maxd, Double_t dz[2], Double_t covar[3]) {
1754 : //
1755 : // Propagate this track to the DCA to vertex "vtx",
1756 : // if the (rough) transverse impact parameter is not bigger then "maxd".
1757 : //
1758 : // This function takes into account all three components of the magnetic
1759 : // field given by the b[3] arument (kG)
1760 : //
1761 : // a) The track gets extapolated to the DCA to the vertex.
1762 : // b) The impact parameters and their covariance matrix are calculated.
1763 : //
1764 : // In the case of success, the returned value is kTRUE
1765 : // (otherwise, it's kFALSE)
1766 : //
1767 696 : Double_t alpha=GetAlpha();
1768 348 : Double_t sn=TMath::Sin(alpha), cs=TMath::Cos(alpha);
1769 348 : Double_t x=GetX(), y=GetParameter()[0], snp=GetParameter()[2];
1770 348 : Double_t xv= vtx->GetX()*cs + vtx->GetY()*sn;
1771 348 : Double_t yv=-vtx->GetX()*sn + vtx->GetY()*cs, zv=vtx->GetZ();
1772 348 : x-=xv; y-=yv;
1773 :
1774 : //Estimate the impact parameter neglecting the track curvature
1775 348 : Double_t d=TMath::Abs(x*snp - y*TMath::Sqrt((1.-snp)*(1.+snp)));
1776 348 : if (d > maxd) return kFALSE;
1777 :
1778 : //Propagate to the DCA
1779 348 : Double_t crv=GetC(b[2]);
1780 348 : if (TMath::Abs(b[2]) < kAlmost0Field) crv=0.;
1781 :
1782 348 : Double_t tgfv=-(crv*x - snp)/(crv*y + TMath::Sqrt((1.-snp)*(1.+snp)));
1783 348 : sn=tgfv/TMath::Sqrt(1.+ tgfv*tgfv); cs=TMath::Sqrt((1.-sn)*(1.+sn));
1784 696 : if (TMath::Abs(tgfv)>0.) cs = sn/tgfv;
1785 : else cs=1.;
1786 :
1787 348 : x = xv*cs + yv*sn;
1788 348 : yv=-xv*sn + yv*cs; xv=x;
1789 :
1790 348 : if (!PropagateBxByBz(alpha+TMath::ASin(sn),xv,b)) return kFALSE;
1791 :
1792 348 : if (dz==0) return kTRUE;
1793 348 : dz[0] = GetParameter()[0] - yv;
1794 348 : dz[1] = GetParameter()[1] - zv;
1795 :
1796 348 : if (covar==0) return kTRUE;
1797 348 : Double_t cov[6]; vtx->GetCovarianceMatrix(cov);
1798 :
1799 : //***** Improvements by A.Dainese
1800 348 : alpha=GetAlpha(); sn=TMath::Sin(alpha); cs=TMath::Cos(alpha);
1801 348 : Double_t s2ylocvtx = cov[0]*sn*sn + cov[2]*cs*cs - 2.*cov[1]*cs*sn;
1802 348 : covar[0] = GetCovariance()[0] + s2ylocvtx; // neglecting correlations
1803 348 : covar[1] = GetCovariance()[1]; // between (x,y) and z
1804 348 : covar[2] = GetCovariance()[2] + cov[5]; // in vertex's covariance matrix
1805 : //*****
1806 :
1807 : return kTRUE;
1808 696 : }
1809 :
1810 : void AliExternalTrackParam::GetDirection(Double_t d[3]) const {
1811 : //----------------------------------------------------------------
1812 : // This function returns a unit vector along the track direction
1813 : // in the global coordinate system.
1814 : //----------------------------------------------------------------
1815 16 : Double_t cs=TMath::Cos(fAlpha), sn=TMath::Sin(fAlpha);
1816 8 : Double_t snp=fP[2];
1817 8 : Double_t csp =TMath::Sqrt((1.-snp)*(1.+snp));
1818 8 : Double_t norm=TMath::Sqrt(1.+ fP[3]*fP[3]);
1819 8 : d[0]=(csp*cs - snp*sn)/norm;
1820 8 : d[1]=(snp*cs + csp*sn)/norm;
1821 8 : d[2]=fP[3]/norm;
1822 8 : }
1823 :
1824 : Bool_t AliExternalTrackParam::GetPxPyPz(Double_t p[3]) const {
1825 : //---------------------------------------------------------------------
1826 : // This function returns the global track momentum components
1827 : // Results for (nearly) straight tracks are meaningless !
1828 : //---------------------------------------------------------------------
1829 614636 : p[0]=fP[4]; p[1]=fP[2]; p[2]=fP[3];
1830 307318 : return Local2GlobalMomentum(p,fAlpha);
1831 : }
1832 :
1833 : Double_t AliExternalTrackParam::Px() const {
1834 : //---------------------------------------------------------------------
1835 : // Returns x-component of momentum
1836 : // Result for (nearly) straight tracks is meaningless !
1837 : //---------------------------------------------------------------------
1838 :
1839 0 : Double_t p[3]={kVeryBig,kVeryBig,kVeryBig};
1840 0 : GetPxPyPz(p);
1841 :
1842 0 : return p[0];
1843 0 : }
1844 :
1845 : Double_t AliExternalTrackParam::Py() const {
1846 : //---------------------------------------------------------------------
1847 : // Returns y-component of momentum
1848 : // Result for (nearly) straight tracks is meaningless !
1849 : //---------------------------------------------------------------------
1850 :
1851 0 : Double_t p[3]={kVeryBig,kVeryBig,kVeryBig};
1852 0 : GetPxPyPz(p);
1853 :
1854 0 : return p[1];
1855 0 : }
1856 :
1857 : Double_t AliExternalTrackParam::Xv() const {
1858 : //---------------------------------------------------------------------
1859 : // Returns x-component of first track point
1860 : //---------------------------------------------------------------------
1861 :
1862 0 : Double_t r[3]={0.,0.,0.};
1863 0 : GetXYZ(r);
1864 :
1865 0 : return r[0];
1866 0 : }
1867 :
1868 : Double_t AliExternalTrackParam::Yv() const {
1869 : //---------------------------------------------------------------------
1870 : // Returns y-component of first track point
1871 : //---------------------------------------------------------------------
1872 :
1873 0 : Double_t r[3]={0.,0.,0.};
1874 0 : GetXYZ(r);
1875 :
1876 0 : return r[1];
1877 0 : }
1878 :
1879 : Double_t AliExternalTrackParam::Theta() const {
1880 : // return theta angle of momentum
1881 :
1882 564 : return 0.5*TMath::Pi() - TMath::ATan(fP[3]);
1883 : }
1884 :
1885 : Double_t AliExternalTrackParam::Phi() const {
1886 : //---------------------------------------------------------------------
1887 : // Returns the azimuthal angle of momentum
1888 : // 0 <= phi < 2*pi
1889 : //---------------------------------------------------------------------
1890 :
1891 2630 : Double_t phi=TMath::ASin(fP[2]) + fAlpha;
1892 1801 : if (phi<0.) phi+=2.*TMath::Pi();
1893 829 : else if (phi>=2.*TMath::Pi()) phi-=2.*TMath::Pi();
1894 :
1895 1315 : return phi;
1896 : }
1897 :
1898 : Double_t AliExternalTrackParam::PhiPos() const {
1899 : //---------------------------------------------------------------------
1900 : // Returns the azimuthal angle of position
1901 : // 0 <= phi < 2*pi
1902 : //---------------------------------------------------------------------
1903 272 : Double_t r[3]={0.,0.,0.};
1904 136 : GetXYZ(r);
1905 136 : Double_t phi=TMath::ATan2(r[1],r[0]);
1906 184 : if (phi<0.) phi+=2.*TMath::Pi();
1907 :
1908 136 : return phi;
1909 136 : }
1910 :
1911 : Double_t AliExternalTrackParam::M() const {
1912 : // return particle mass
1913 :
1914 : // No mass information available so far.
1915 : // Redifine in derived class!
1916 :
1917 0 : return -999.;
1918 : }
1919 :
1920 : Double_t AliExternalTrackParam::E() const {
1921 : // return particle energy
1922 :
1923 : // No PID information available so far.
1924 : // Redifine in derived class!
1925 :
1926 0 : return -999.;
1927 : }
1928 :
1929 : Double_t AliExternalTrackParam::Eta() const {
1930 : // return pseudorapidity
1931 :
1932 504 : return -TMath::Log(TMath::Tan(0.5 * Theta()));
1933 : }
1934 :
1935 : Double_t AliExternalTrackParam::Y() const {
1936 : // return rapidity
1937 :
1938 : // No PID information available so far.
1939 : // Redifine in derived class!
1940 :
1941 0 : return -999.;
1942 : }
1943 :
1944 : Bool_t AliExternalTrackParam::GetXYZ(Double_t *r) const {
1945 : //---------------------------------------------------------------------
1946 : // This function returns the global track position
1947 : //---------------------------------------------------------------------
1948 2185030 : r[0]=fX; r[1]=fP[0]; r[2]=fP[1];
1949 1092515 : return Local2GlobalPosition(r,fAlpha);
1950 : }
1951 :
1952 : Bool_t AliExternalTrackParam::GetCovarianceXYZPxPyPz(Double_t cv[21]) const {
1953 : //---------------------------------------------------------------------
1954 : // This function returns the global covariance matrix of the track params
1955 : //
1956 : // Cov(x,x) ... : cv[0]
1957 : // Cov(y,x) ... : cv[1] cv[2]
1958 : // Cov(z,x) ... : cv[3] cv[4] cv[5]
1959 : // Cov(px,x)... : cv[6] cv[7] cv[8] cv[9]
1960 : // Cov(py,x)... : cv[10] cv[11] cv[12] cv[13] cv[14]
1961 : // Cov(pz,x)... : cv[15] cv[16] cv[17] cv[18] cv[19] cv[20]
1962 : //
1963 : // Results for (nearly) straight tracks are meaningless !
1964 : //---------------------------------------------------------------------
1965 494 : if (TMath::Abs(fP[4])<=kAlmost0) {
1966 0 : for (Int_t i=0; i<21; i++) cv[i]=0.;
1967 0 : return kFALSE;
1968 : }
1969 247 : if (TMath::Abs(fP[2]) > kAlmost1) {
1970 0 : for (Int_t i=0; i<21; i++) cv[i]=0.;
1971 0 : return kFALSE;
1972 : }
1973 247 : Double_t pt=1./TMath::Abs(fP[4]);
1974 247 : Double_t cs=TMath::Cos(fAlpha), sn=TMath::Sin(fAlpha);
1975 247 : Double_t r=TMath::Sqrt((1.-fP[2])*(1.+fP[2]));
1976 :
1977 247 : Double_t m00=-sn, m10=cs;
1978 247 : Double_t m23=-pt*(sn + fP[2]*cs/r), m43=-pt*pt*(r*cs - fP[2]*sn);
1979 247 : Double_t m24= pt*(cs - fP[2]*sn/r), m44=-pt*pt*(r*sn + fP[2]*cs);
1980 247 : Double_t m35=pt, m45=-pt*pt*fP[3];
1981 :
1982 247 : m43*=GetSign();
1983 247 : m44*=GetSign();
1984 247 : m45*=GetSign();
1985 :
1986 247 : cv[0 ] = fC[0]*m00*m00;
1987 247 : cv[1 ] = fC[0]*m00*m10;
1988 247 : cv[2 ] = fC[0]*m10*m10;
1989 247 : cv[3 ] = fC[1]*m00;
1990 247 : cv[4 ] = fC[1]*m10;
1991 247 : cv[5 ] = fC[2];
1992 247 : cv[6 ] = m00*(fC[3]*m23 + fC[10]*m43);
1993 247 : cv[7 ] = m10*(fC[3]*m23 + fC[10]*m43);
1994 247 : cv[8 ] = fC[4]*m23 + fC[11]*m43;
1995 247 : cv[9 ] = m23*(fC[5]*m23 + fC[12]*m43) + m43*(fC[12]*m23 + fC[14]*m43);
1996 247 : cv[10] = m00*(fC[3]*m24 + fC[10]*m44);
1997 247 : cv[11] = m10*(fC[3]*m24 + fC[10]*m44);
1998 247 : cv[12] = fC[4]*m24 + fC[11]*m44;
1999 247 : cv[13] = m23*(fC[5]*m24 + fC[12]*m44) + m43*(fC[12]*m24 + fC[14]*m44);
2000 247 : cv[14] = m24*(fC[5]*m24 + fC[12]*m44) + m44*(fC[12]*m24 + fC[14]*m44);
2001 247 : cv[15] = m00*(fC[6]*m35 + fC[10]*m45);
2002 247 : cv[16] = m10*(fC[6]*m35 + fC[10]*m45);
2003 247 : cv[17] = fC[7]*m35 + fC[11]*m45;
2004 247 : cv[18] = m23*(fC[8]*m35 + fC[12]*m45) + m43*(fC[13]*m35 + fC[14]*m45);
2005 247 : cv[19] = m24*(fC[8]*m35 + fC[12]*m45) + m44*(fC[13]*m35 + fC[14]*m45);
2006 247 : cv[20] = m35*(fC[9]*m35 + fC[13]*m45) + m45*(fC[13]*m35 + fC[14]*m45);
2007 :
2008 : return kTRUE;
2009 247 : }
2010 :
2011 :
2012 : Bool_t
2013 : AliExternalTrackParam::GetPxPyPzAt(Double_t x, Double_t b, Double_t *p) const {
2014 : //---------------------------------------------------------------------
2015 : // This function returns the global track momentum extrapolated to
2016 : // the radial position "x" (cm) in the magnetic field "b" (kG)
2017 : //---------------------------------------------------------------------
2018 360 : p[0]=fP[4];
2019 180 : p[1]=fP[2]+(x-fX)*GetC(b);
2020 180 : p[2]=fP[3];
2021 180 : return Local2GlobalMomentum(p,fAlpha);
2022 : }
2023 :
2024 : Bool_t AliExternalTrackParam::GetYZAt(Double_t x, Double_t b, Double_t *yz) const
2025 : {
2026 : //---------------------------------------------------------------------
2027 : // This function returns the local Y,Z-coordinates of the intersection
2028 : // point between this track and the reference plane "x" (cm).
2029 : // Magnetic field "b" (kG)
2030 : //---------------------------------------------------------------------
2031 0 : double dx=x-fX;
2032 0 : if(TMath::Abs(dx)<=kAlmost0) {
2033 0 : yz[0] = fP[0];
2034 0 : yz[1] = fP[1];
2035 0 : return kTRUE;
2036 : }
2037 0 : double crv=GetC(b);
2038 0 : double f1=fP[2], x2r = crv*dx, f2=f1 + x2r;
2039 0 : if (TMath::Abs(f1) >= kAlmost1 ||
2040 0 : TMath::Abs(f2) >= kAlmost1 ||
2041 0 : TMath::Abs(fP[4])< kAlmost0) return kFALSE;
2042 0 : double r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2));
2043 0 : double dy2dx = (f1+f2)/(r1+r2);
2044 0 : yz[0] = fP[0] + dx*dy2dx;
2045 0 : if (TMath::Abs(x2r)<0.05) yz[1] = fP[1] + dx*(r2 + f2*dy2dx)*fP[3];
2046 : else {
2047 0 : double chord = dx*TMath::Sqrt(1+dy2dx*dy2dx); // distance from old position to new one
2048 0 : double rot = 2*TMath::ASin(0.5*chord*crv); // angular difference seen from the circle center
2049 0 : yz[1] = fP[1] + rot/crv*fP[3];
2050 : }
2051 : return kTRUE;
2052 0 : }
2053 :
2054 :
2055 : Bool_t
2056 : AliExternalTrackParam::GetYAt(Double_t x, Double_t b, Double_t &y) const {
2057 : //---------------------------------------------------------------------
2058 : // This function returns the local Y-coordinate of the intersection
2059 : // point between this track and the reference plane "x" (cm).
2060 : // Magnetic field "b" (kG)
2061 : //---------------------------------------------------------------------
2062 23154 : Double_t dx=x-fX;
2063 11577 : if(TMath::Abs(dx)<=kAlmost0) {y=fP[0]; return kTRUE;}
2064 :
2065 11577 : Double_t f1=fP[2], f2=f1 + dx*GetC(b);
2066 :
2067 11577 : if (TMath::Abs(f1) >= kAlmost1) return kFALSE;
2068 11637 : if (TMath::Abs(f2) >= kAlmost1) return kFALSE;
2069 :
2070 11517 : Double_t r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2));
2071 11517 : y = fP[0] + dx*(f1+f2)/(r1+r2);
2072 : return kTRUE;
2073 11577 : }
2074 :
2075 : Bool_t
2076 : AliExternalTrackParam::GetZAt(Double_t x, Double_t b, Double_t &z) const {
2077 : //---------------------------------------------------------------------
2078 : // This function returns the local Z-coordinate of the intersection
2079 : // point between this track and the reference plane "x" (cm).
2080 : // Magnetic field "b" (kG)
2081 : //---------------------------------------------------------------------
2082 259370 : Double_t dx=x-fX;
2083 129685 : if(TMath::Abs(dx)<=kAlmost0) {z=fP[1]; return kTRUE;}
2084 :
2085 129685 : Double_t crv=GetC(b);
2086 129685 : Double_t x2r = crv*dx;
2087 129685 : Double_t f1=fP[2], f2=f1 + x2r;
2088 :
2089 129685 : if (TMath::Abs(f1) >= kAlmost1) return kFALSE;
2090 129821 : if (TMath::Abs(f2) >= kAlmost1) return kFALSE;
2091 :
2092 129549 : Double_t r1=sqrt((1.-f1)*(1.+f1)), r2=sqrt((1.-f2)*(1.+f2));
2093 129549 : double dy2dx = (f1+f2)/(r1+r2);
2094 129549 : if (TMath::Abs(x2r)<0.05) {
2095 128703 : z = fP[1] + dx*(r2 + f2*dy2dx)*fP[3]; // Many thanks to P.Hristov !
2096 128703 : }
2097 : else {
2098 : // for small dx/R the linear apporximation of the arc by the segment is OK,
2099 : // but at large dx/R the error is very large and leads to incorrect Z propagation
2100 : // angle traversed delta = 2*asin(dist_start_end / R / 2), hence the arc is: R*deltaPhi
2101 : // The dist_start_end is obtained from sqrt(dx^2+dy^2) = x/(r1+r2)*sqrt(2+f1*f2+r1*r2)
2102 : // Similarly, the rotation angle in linear in dx only for dx<<R
2103 846 : double chord = dx*TMath::Sqrt(1+dy2dx*dy2dx); // distance from old position to new one
2104 846 : double rot = 2*TMath::ASin(0.5*chord*crv); // angular difference seen from the circle center
2105 846 : z = fP[1] + rot/crv*fP[3];
2106 : }
2107 : return kTRUE;
2108 129685 : }
2109 :
2110 : Bool_t
2111 : AliExternalTrackParam::GetXYZAt(Double_t x, Double_t b, Double_t *r) const {
2112 : //---------------------------------------------------------------------
2113 : // This function returns the global track position extrapolated to
2114 : // the radial position "x" (cm) in the magnetic field "b" (kG)
2115 : //---------------------------------------------------------------------
2116 86740 : Double_t dx=x-fX;
2117 43468 : if(TMath::Abs(dx)<=kAlmost0) return GetXYZ(r);
2118 :
2119 43272 : Double_t crv=GetC(b);
2120 43272 : Double_t x2r = crv*dx;
2121 43272 : Double_t f1=fP[2], f2=f1 + dx*crv;
2122 :
2123 43272 : if (TMath::Abs(f1) >= kAlmost1) return kFALSE;
2124 45236 : if (TMath::Abs(f2) >= kAlmost1) return kFALSE;
2125 :
2126 41308 : Double_t r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2));
2127 41308 : double dy2dx = (f1+f2)/(r1+r2);
2128 41308 : r[0] = x;
2129 41308 : r[1] = fP[0] + dx*dy2dx;
2130 41308 : if (TMath::Abs(x2r)<0.05) {
2131 38379 : r[2] = fP[1] + dx*(r2 + f2*dy2dx)*fP[3];//Thanks to Andrea & Peter
2132 38379 : }
2133 : else {
2134 : // for small dx/R the linear apporximation of the arc by the segment is OK,
2135 : // but at large dx/R the error is very large and leads to incorrect Z propagation
2136 : // angle traversed delta = 2*asin(dist_start_end / R / 2), hence the arc is: R*deltaPhi
2137 : // The dist_start_end is obtained from sqrt(dx^2+dy^2) = x/(r1+r2)*sqrt(2+f1*f2+r1*r2)
2138 : // Similarly, the rotation angle in linear in dx only for dx<<R
2139 2929 : double chord = dx*TMath::Sqrt(1+dy2dx*dy2dx); // distance from old position to new one
2140 2929 : double rot = 2*TMath::ASin(0.5*chord*crv); // angular difference seen from the circle center
2141 2929 : r[2] = fP[1] + rot/crv*fP[3];
2142 : }
2143 :
2144 41308 : return Local2GlobalPosition(r,fAlpha);
2145 43370 : }
2146 :
2147 : //_____________________________________________________________________________
2148 : void AliExternalTrackParam::Print(Option_t* /*option*/) const
2149 : {
2150 : // print the parameters and the covariance matrix
2151 :
2152 0 : printf("AliExternalTrackParam: x = %-12g alpha = %-12g\n", fX, fAlpha);
2153 0 : printf(" parameters: %12g %12g %12g %12g %12g\n",
2154 0 : fP[0], fP[1], fP[2], fP[3], fP[4]);
2155 0 : printf(" covariance: %12g\n", fC[0]);
2156 0 : printf(" %12g %12g\n", fC[1], fC[2]);
2157 0 : printf(" %12g %12g %12g\n", fC[3], fC[4], fC[5]);
2158 0 : printf(" %12g %12g %12g %12g\n",
2159 0 : fC[6], fC[7], fC[8], fC[9]);
2160 0 : printf(" %12g %12g %12g %12g %12g\n",
2161 0 : fC[10], fC[11], fC[12], fC[13], fC[14]);
2162 0 : }
2163 :
2164 : Double_t AliExternalTrackParam::GetSnpAt(Double_t x,Double_t b) const {
2165 : //
2166 : // Get sinus at given x
2167 : //
2168 0 : Double_t crv=GetC(b);
2169 0 : if (TMath::Abs(b) < kAlmost0Field) crv=0.;
2170 0 : Double_t dx = x-fX;
2171 0 : Double_t res = fP[2]+dx*crv;
2172 0 : return res;
2173 : }
2174 :
2175 : Bool_t AliExternalTrackParam::GetDistance(AliExternalTrackParam *param2, Double_t x, Double_t dist[3], Double_t bz){
2176 : //------------------------------------------------------------------------
2177 : // Get the distance between two tracks at the local position x
2178 : // working in the local frame of this track.
2179 : // Origin : Marian.Ivanov@cern.ch
2180 : //-----------------------------------------------------------------------
2181 0 : Double_t xyz[3];
2182 0 : Double_t xyz2[3];
2183 0 : xyz[0]=x;
2184 0 : if (!GetYAt(x,bz,xyz[1])) return kFALSE;
2185 0 : if (!GetZAt(x,bz,xyz[2])) return kFALSE;
2186 : //
2187 : //
2188 0 : if (TMath::Abs(GetAlpha()-param2->GetAlpha())<kAlmost0){
2189 0 : xyz2[0]=x;
2190 0 : if (!param2->GetYAt(x,bz,xyz2[1])) return kFALSE;
2191 0 : if (!param2->GetZAt(x,bz,xyz2[2])) return kFALSE;
2192 : }else{
2193 : //
2194 0 : Double_t xyz1[3];
2195 0 : Double_t dfi = param2->GetAlpha()-GetAlpha();
2196 0 : Double_t ca = TMath::Cos(dfi), sa = TMath::Sin(dfi);
2197 0 : xyz2[0] = xyz[0]*ca+xyz[1]*sa;
2198 0 : xyz2[1] = -xyz[0]*sa+xyz[1]*ca;
2199 : //
2200 0 : xyz1[0]=xyz2[0];
2201 0 : if (!param2->GetYAt(xyz2[0],bz,xyz1[1])) return kFALSE;
2202 0 : if (!param2->GetZAt(xyz2[0],bz,xyz1[2])) return kFALSE;
2203 : //
2204 0 : xyz2[0] = xyz1[0]*ca-xyz1[1]*sa;
2205 0 : xyz2[1] = +xyz1[0]*sa+xyz1[1]*ca;
2206 0 : xyz2[2] = xyz1[2];
2207 0 : }
2208 0 : dist[0] = xyz[0]-xyz2[0];
2209 0 : dist[1] = xyz[1]-xyz2[1];
2210 0 : dist[2] = xyz[2]-xyz2[2];
2211 :
2212 0 : return kTRUE;
2213 0 : }
2214 :
2215 :
2216 : //
2217 : // Draw functionality.
2218 : // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
2219 : //
2220 :
2221 : void AliExternalTrackParam::DrawTrack(Float_t magf, Float_t minR, Float_t maxR, Float_t stepR){
2222 : //
2223 : // Draw track line
2224 : //
2225 0 : if (minR>maxR) return ;
2226 0 : if (stepR<=0) return ;
2227 0 : Int_t npoints = TMath::Nint((maxR-minR)/stepR)+1;
2228 0 : if (npoints<1) return;
2229 0 : TPolyMarker3D *polymarker = new TPolyMarker3D(npoints);
2230 0 : FillPolymarker(polymarker, magf,minR,maxR,stepR);
2231 0 : polymarker->Draw();
2232 0 : }
2233 :
2234 : //
2235 : void AliExternalTrackParam::FillPolymarker(TPolyMarker3D *pol, Float_t magF, Float_t minR, Float_t maxR, Float_t stepR){
2236 : //
2237 : // Fill points in the polymarker
2238 : //
2239 : Int_t counter=0;
2240 0 : for (Double_t r=minR; r<maxR; r+=stepR){
2241 0 : Double_t point[3];
2242 0 : GetXYZAt(r,magF,point);
2243 0 : pol->SetPoint(counter,point[0],point[1], point[2]);
2244 : // printf("xyz\t%f\t%f\t%f\n",point[0], point[1],point[2]);
2245 0 : counter++;
2246 0 : }
2247 0 : }
2248 :
2249 : Int_t AliExternalTrackParam::GetIndex(Int_t i, Int_t j){
2250 : //
2251 20340 : Int_t min = TMath::Min(i,j);
2252 10170 : Int_t max = TMath::Max(i,j);
2253 :
2254 10170 : return min+(max+1)*max/2;
2255 : }
2256 :
2257 :
2258 : void AliExternalTrackParam::g3helx3(Double_t qfield,
2259 : Double_t step,
2260 : Double_t vect[7]) {
2261 : /******************************************************************
2262 : * *
2263 : * GEANT3 tracking routine in a constant field oriented *
2264 : * along axis 3 *
2265 : * Tracking is performed with a conventional *
2266 : * helix step method *
2267 : * *
2268 : * Authors R.Brun, M.Hansroul ********* *
2269 : * Rewritten V.Perevoztchikov *
2270 : * *
2271 : * Rewritten in C++ by I.Belikov *
2272 : * *
2273 : * qfield (kG) - particle charge times magnetic field *
2274 : * step (cm) - step length along the helix *
2275 : * vect[7](cm,GeV/c) - input/output x, y, z, px/p, py/p ,pz/p, p *
2276 : * *
2277 : ******************************************************************/
2278 : const Int_t ix=0, iy=1, iz=2, ipx=3, ipy=4, ipz=5, ipp=6;
2279 570186 : const Double_t kOvSqSix=TMath::Sqrt(1./6.);
2280 :
2281 285093 : Double_t cosx=vect[ipx], cosy=vect[ipy], cosz=vect[ipz];
2282 :
2283 285093 : Double_t rho = qfield*kB2C/vect[ipp];
2284 285093 : Double_t tet = rho*step;
2285 :
2286 : Double_t tsint, sintt, sint, cos1t;
2287 285093 : if (TMath::Abs(tet) > 0.03) {
2288 3827 : sint = TMath::Sin(tet);
2289 3827 : sintt = sint/tet;
2290 3827 : tsint = (tet - sint)/tet;
2291 3827 : Double_t t=TMath::Sin(0.5*tet);
2292 3827 : cos1t = 2*t*t/tet;
2293 3827 : } else {
2294 281266 : tsint = tet*tet/6.;
2295 281266 : sintt = (1.-tet*kOvSqSix)*(1.+tet*kOvSqSix); // 1.- tsint;
2296 281266 : sint = tet*sintt;
2297 281266 : cos1t = 0.5*tet;
2298 : }
2299 :
2300 285093 : Double_t f1 = step*sintt;
2301 285093 : Double_t f2 = step*cos1t;
2302 285093 : Double_t f3 = step*tsint*cosz;
2303 285093 : Double_t f4 = -tet*cos1t;
2304 : Double_t f5 = sint;
2305 :
2306 285093 : vect[ix] += f1*cosx - f2*cosy;
2307 285093 : vect[iy] += f1*cosy + f2*cosx;
2308 285093 : vect[iz] += f1*cosz + f3;
2309 :
2310 285093 : vect[ipx] += f4*cosx - f5*cosy;
2311 285093 : vect[ipy] += f4*cosy + f5*cosx;
2312 :
2313 285093 : }
2314 :
2315 : Bool_t AliExternalTrackParam::PropagateToBxByBz(Double_t xk, const Double_t b[3]) {
2316 : //----------------------------------------------------------------
2317 : // Extrapolate this track to the plane X=xk in the field b[].
2318 : //
2319 : // X [cm] is in the "tracking coordinate system" of this track.
2320 : // b[]={Bx,By,Bz} [kG] is in the Global coordidate system.
2321 : //----------------------------------------------------------------
2322 :
2323 560416 : Double_t dx=xk-fX;
2324 280492 : if (TMath::Abs(dx)<=kAlmost0) return kTRUE;
2325 279924 : if (TMath::Abs(fP[4])<=kAlmost0) return kFALSE;
2326 : // Do not propagate tracks outside the ALICE detector
2327 559848 : if (TMath::Abs(dx)>1e5 ||
2328 279924 : TMath::Abs(GetY())>1e5 ||
2329 279924 : TMath::Abs(GetZ())>1e5) {
2330 0 : AliWarning(Form("Anomalous track, target X:%f",xk));
2331 0 : Print();
2332 0 : return kFALSE;
2333 : }
2334 :
2335 279924 : Double_t crv=GetC(b[2]);
2336 279926 : if (TMath::Abs(b[2]) < kAlmost0Field) crv=0.;
2337 :
2338 279924 : Double_t x2r = crv*dx;
2339 279924 : Double_t f1=fP[2], f2=f1 + x2r;
2340 279924 : if (TMath::Abs(f1) >= kAlmost1) return kFALSE;
2341 279959 : if (TMath::Abs(f2) >= kAlmost1) return kFALSE;
2342 :
2343 :
2344 : // Estimate the covariance matrix
2345 279889 : Double_t &fP3=fP[3], &fP4=fP[4];
2346 : Double_t
2347 279889 : &fC00=fC[0],
2348 279889 : &fC10=fC[1], &fC11=fC[2],
2349 279889 : &fC20=fC[3], &fC21=fC[4], &fC22=fC[5],
2350 279889 : &fC30=fC[6], &fC31=fC[7], &fC32=fC[8], &fC33=fC[9],
2351 279889 : &fC40=fC[10], &fC41=fC[11], &fC42=fC[12], &fC43=fC[13], &fC44=fC[14];
2352 :
2353 279889 : Double_t r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2));
2354 :
2355 : //f = F - 1
2356 : /*
2357 : Double_t f02= dx/(r1*r1*r1); Double_t cc=crv/fP4;
2358 : Double_t f04=0.5*dx*dx/(r1*r1*r1); f04*=cc;
2359 : Double_t f12= dx*fP3*f1/(r1*r1*r1);
2360 : Double_t f14=0.5*dx*dx*fP3*f1/(r1*r1*r1); f14*=cc;
2361 : Double_t f13= dx/r1;
2362 : Double_t f24= dx; f24*=cc;
2363 : */
2364 279889 : Double_t rinv = 1./r1;
2365 279889 : Double_t r3inv = rinv*rinv*rinv;
2366 279889 : Double_t f24= x2r/fP4;
2367 279889 : Double_t f02= dx*r3inv;
2368 279889 : Double_t f04=0.5*f24*f02;
2369 279889 : Double_t f12= f02*fP3*f1;
2370 279889 : Double_t f14=0.5*f24*f02*fP3*f1;
2371 279889 : Double_t f13= dx*rinv;
2372 :
2373 : //b = C*ft
2374 279889 : Double_t b00=f02*fC20 + f04*fC40, b01=f12*fC20 + f14*fC40 + f13*fC30;
2375 279889 : Double_t b02=f24*fC40;
2376 279889 : Double_t b10=f02*fC21 + f04*fC41, b11=f12*fC21 + f14*fC41 + f13*fC31;
2377 279889 : Double_t b12=f24*fC41;
2378 279889 : Double_t b20=f02*fC22 + f04*fC42, b21=f12*fC22 + f14*fC42 + f13*fC32;
2379 279889 : Double_t b22=f24*fC42;
2380 279889 : Double_t b40=f02*fC42 + f04*fC44, b41=f12*fC42 + f14*fC44 + f13*fC43;
2381 279889 : Double_t b42=f24*fC44;
2382 279889 : Double_t b30=f02*fC32 + f04*fC43, b31=f12*fC32 + f14*fC43 + f13*fC33;
2383 279889 : Double_t b32=f24*fC43;
2384 :
2385 : //a = f*b = f*C*ft
2386 279889 : Double_t a00=f02*b20+f04*b40,a01=f02*b21+f04*b41,a02=f02*b22+f04*b42;
2387 279889 : Double_t a11=f12*b21+f14*b41+f13*b31,a12=f12*b22+f14*b42+f13*b32;
2388 279889 : Double_t a22=f24*b42;
2389 :
2390 : //F*C*Ft = C + (b + bt + a)
2391 279889 : fC00 += b00 + b00 + a00;
2392 279889 : fC10 += b10 + b01 + a01;
2393 279889 : fC20 += b20 + b02 + a02;
2394 279889 : fC30 += b30;
2395 279889 : fC40 += b40;
2396 279889 : fC11 += b11 + b11 + a11;
2397 279889 : fC21 += b21 + b12 + a12;
2398 279889 : fC31 += b31;
2399 279889 : fC41 += b41;
2400 279889 : fC22 += b22 + b22 + a22;
2401 279889 : fC32 += b32;
2402 279889 : fC42 += b42;
2403 :
2404 279889 : CheckCovariance();
2405 :
2406 : // Appoximate step length
2407 279889 : double dy2dx = (f1+f2)/(r1+r2);
2408 839667 : Double_t step = (TMath::Abs(x2r)<0.05) ? dx*TMath::Abs(r2 + f2*dy2dx) // chord
2409 1024 : : 2.*TMath::ASin(0.5*dx*TMath::Sqrt(1.+dy2dx*dy2dx)*crv)/crv; // arc
2410 279889 : step *= TMath::Sqrt(1.+ GetTgl()*GetTgl());
2411 :
2412 : // Get the track's (x,y,z) and (px,py,pz) in the Global System
2413 279889 : Double_t r[3]; GetXYZ(r);
2414 279889 : Double_t p[3]; GetPxPyPz(p);
2415 279889 : Double_t pp=GetP();
2416 279889 : p[0] /= pp;
2417 279889 : p[1] /= pp;
2418 279889 : p[2] /= pp;
2419 :
2420 :
2421 : // Rotate to the system where Bx=By=0.
2422 279889 : Double_t bt=TMath::Sqrt(b[0]*b[0] + b[1]*b[1]);
2423 : Double_t cosphi=1., sinphi=0.;
2424 559776 : if (bt > kAlmost0) {cosphi=b[0]/bt; sinphi=b[1]/bt;}
2425 279889 : Double_t bb=TMath::Sqrt(b[0]*b[0] + b[1]*b[1] + b[2]*b[2]);
2426 : Double_t costet=1., sintet=0.;
2427 559778 : if (bb > kAlmost0) {costet=b[2]/bb; sintet=bt/bb;}
2428 279889 : Double_t vect[7];
2429 :
2430 279889 : vect[0] = costet*cosphi*r[0] + costet*sinphi*r[1] - sintet*r[2];
2431 279889 : vect[1] = -sinphi*r[0] + cosphi*r[1];
2432 279889 : vect[2] = sintet*cosphi*r[0] + sintet*sinphi*r[1] + costet*r[2];
2433 :
2434 279889 : vect[3] = costet*cosphi*p[0] + costet*sinphi*p[1] - sintet*p[2];
2435 279889 : vect[4] = -sinphi*p[0] + cosphi*p[1];
2436 279889 : vect[5] = sintet*cosphi*p[0] + sintet*sinphi*p[1] + costet*p[2];
2437 :
2438 279889 : vect[6] = pp;
2439 :
2440 :
2441 : // Do the helix step
2442 279889 : g3helx3(GetSign()*bb,step,vect);
2443 :
2444 :
2445 : // Rotate back to the Global System
2446 279889 : r[0] = cosphi*costet*vect[0] - sinphi*vect[1] + cosphi*sintet*vect[2];
2447 279889 : r[1] = sinphi*costet*vect[0] + cosphi*vect[1] + sinphi*sintet*vect[2];
2448 279889 : r[2] = -sintet*vect[0] + costet*vect[2];
2449 :
2450 279889 : p[0] = cosphi*costet*vect[3] - sinphi*vect[4] + cosphi*sintet*vect[5];
2451 279889 : p[1] = sinphi*costet*vect[3] + cosphi*vect[4] + sinphi*sintet*vect[5];
2452 279889 : p[2] = -sintet*vect[3] + costet*vect[5];
2453 :
2454 :
2455 : // Rotate back to the Tracking System
2456 279889 : Double_t cosalp = TMath::Cos(fAlpha);
2457 279889 : Double_t sinalp =-TMath::Sin(fAlpha);
2458 :
2459 : Double_t
2460 279889 : t = cosalp*r[0] - sinalp*r[1];
2461 279889 : r[1] = sinalp*r[0] + cosalp*r[1];
2462 279889 : r[0] = t;
2463 :
2464 279889 : t = cosalp*p[0] - sinalp*p[1];
2465 279889 : p[1] = sinalp*p[0] + cosalp*p[1];
2466 279889 : p[0] = t;
2467 :
2468 :
2469 : // Do the final correcting step to the target plane (linear approximation)
2470 279889 : Double_t x=r[0], y=r[1], z=r[2];
2471 279889 : if (TMath::Abs(dx) > kAlmost0) {
2472 279889 : if (TMath::Abs(p[0]) < kAlmost0) return kFALSE;
2473 279889 : dx = xk - r[0];
2474 279889 : x += dx;
2475 279889 : y += p[1]/p[0]*dx;
2476 279889 : z += p[2]/p[0]*dx;
2477 279889 : }
2478 :
2479 :
2480 : // Calculate the track parameters
2481 279889 : t=TMath::Sqrt(p[0]*p[0] + p[1]*p[1]);
2482 279889 : fX = x;
2483 279889 : fP[0] = y;
2484 279889 : fP[1] = z;
2485 279889 : fP[2] = p[1]/t;
2486 279889 : fP[3] = p[2]/t;
2487 279889 : fP[4] = GetSign()/(t*pp);
2488 :
2489 279889 : return kTRUE;
2490 560097 : }
2491 :
2492 : Bool_t AliExternalTrackParam::PropagateParamOnlyBxByBzTo(Double_t xk, const Double_t b[3]) {
2493 : //----------------------------------------------------------------
2494 : // Extrapolate this track params (w/o cov matrix) to the plane X=xk in the field b[].
2495 : //
2496 : // X [cm] is in the "tracking coordinate system" of this track.
2497 : // b[]={Bx,By,Bz} [kG] is in the Global coordidate system.
2498 : //----------------------------------------------------------------
2499 :
2500 10408 : Double_t dx=xk-fX;
2501 5204 : if (TMath::Abs(dx)<=kAlmost0) return kTRUE;
2502 5204 : if (TMath::Abs(fP[4])<=kAlmost0) return kFALSE;
2503 : // Do not propagate tracks outside the ALICE detector
2504 10408 : if (TMath::Abs(dx)>1e5 ||
2505 5204 : TMath::Abs(GetY())>1e5 ||
2506 5204 : TMath::Abs(GetZ())>1e5) {
2507 0 : AliWarning(Form("Anomalous track, target X:%f",xk));
2508 0 : Print();
2509 0 : return kFALSE;
2510 : }
2511 :
2512 5204 : Double_t crv=GetC(b[2]);
2513 5204 : if (TMath::Abs(b[2]) < kAlmost0Field) crv=0.;
2514 :
2515 5204 : Double_t x2r = crv*dx;
2516 5204 : Double_t f1=fP[2], f2=f1 + x2r;
2517 5204 : if (TMath::Abs(f1) >= kAlmost1) return kFALSE;
2518 5204 : if (TMath::Abs(f2) >= kAlmost1) return kFALSE;
2519 : //
2520 5204 : Double_t r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2));
2521 : //
2522 : // Appoximate step length
2523 5204 : double dy2dx = (f1+f2)/(r1+r2);
2524 15612 : Double_t step = (TMath::Abs(x2r)<0.05) ? dx*TMath::Abs(r2 + f2*dy2dx) // chord
2525 280 : : 2.*TMath::ASin(0.5*dx*TMath::Sqrt(1.+dy2dx*dy2dx)*crv)/crv; // arc
2526 5204 : step *= TMath::Sqrt(1.+ GetTgl()*GetTgl());
2527 :
2528 : // Get the track's (x,y,z) and (px,py,pz) in the Global System
2529 5204 : Double_t r[3]; GetXYZ(r);
2530 5204 : Double_t p[3]; GetPxPyPz(p);
2531 5204 : Double_t pp=GetP();
2532 5204 : p[0] /= pp;
2533 5204 : p[1] /= pp;
2534 5204 : p[2] /= pp;
2535 :
2536 : // Rotate to the system where Bx=By=0.
2537 5204 : Double_t bt=TMath::Sqrt(b[0]*b[0] + b[1]*b[1]);
2538 : Double_t cosphi=1., sinphi=0.;
2539 10408 : if (bt > kAlmost0) {cosphi=b[0]/bt; sinphi=b[1]/bt;}
2540 5204 : Double_t bb=TMath::Sqrt(b[0]*b[0] + b[1]*b[1] + b[2]*b[2]);
2541 : Double_t costet=1., sintet=0.;
2542 10408 : if (bb > kAlmost0) {costet=b[2]/bb; sintet=bt/bb;}
2543 5204 : Double_t vect[7];
2544 :
2545 5204 : vect[0] = costet*cosphi*r[0] + costet*sinphi*r[1] - sintet*r[2];
2546 5204 : vect[1] = -sinphi*r[0] + cosphi*r[1];
2547 5204 : vect[2] = sintet*cosphi*r[0] + sintet*sinphi*r[1] + costet*r[2];
2548 :
2549 5204 : vect[3] = costet*cosphi*p[0] + costet*sinphi*p[1] - sintet*p[2];
2550 5204 : vect[4] = -sinphi*p[0] + cosphi*p[1];
2551 5204 : vect[5] = sintet*cosphi*p[0] + sintet*sinphi*p[1] + costet*p[2];
2552 :
2553 5204 : vect[6] = pp;
2554 :
2555 : // Do the helix step
2556 5204 : g3helx3(GetSign()*bb,step,vect);
2557 :
2558 : // Rotate back to the Global System
2559 5204 : r[0] = cosphi*costet*vect[0] - sinphi*vect[1] + cosphi*sintet*vect[2];
2560 5204 : r[1] = sinphi*costet*vect[0] + cosphi*vect[1] + sinphi*sintet*vect[2];
2561 5204 : r[2] = -sintet*vect[0] + costet*vect[2];
2562 :
2563 5204 : p[0] = cosphi*costet*vect[3] - sinphi*vect[4] + cosphi*sintet*vect[5];
2564 5204 : p[1] = sinphi*costet*vect[3] + cosphi*vect[4] + sinphi*sintet*vect[5];
2565 5204 : p[2] = -sintet*vect[3] + costet*vect[5];
2566 :
2567 : // Rotate back to the Tracking System
2568 5204 : Double_t cosalp = TMath::Cos(fAlpha);
2569 5204 : Double_t sinalp =-TMath::Sin(fAlpha);
2570 :
2571 : Double_t
2572 5204 : t = cosalp*r[0] - sinalp*r[1];
2573 5204 : r[1] = sinalp*r[0] + cosalp*r[1];
2574 5204 : r[0] = t;
2575 :
2576 5204 : t = cosalp*p[0] - sinalp*p[1];
2577 5204 : p[1] = sinalp*p[0] + cosalp*p[1];
2578 5204 : p[0] = t;
2579 :
2580 : // Do the final correcting step to the target plane (linear approximation)
2581 5204 : Double_t x=r[0], y=r[1], z=r[2];
2582 5204 : if (TMath::Abs(dx) > kAlmost0) {
2583 5204 : if (TMath::Abs(p[0]) < kAlmost0) return kFALSE;
2584 5204 : dx = xk - r[0];
2585 5204 : x += dx;
2586 5204 : y += p[1]/p[0]*dx;
2587 5204 : z += p[2]/p[0]*dx;
2588 5204 : }
2589 :
2590 :
2591 : // Calculate the track parameters
2592 5204 : t=TMath::Sqrt(p[0]*p[0] + p[1]*p[1]);
2593 5204 : fX = x;
2594 5204 : fP[0] = y;
2595 5204 : fP[1] = z;
2596 5204 : fP[2] = p[1]/t;
2597 5204 : fP[3] = p[2]/t;
2598 5204 : fP[4] = GetSign()/(t*pp);
2599 :
2600 5204 : return kTRUE;
2601 10408 : }
2602 :
2603 :
2604 : Bool_t AliExternalTrackParam::Translate(Double_t *vTrasl,Double_t *covV){
2605 : //
2606 : //Translation: in the event mixing, the tracks can be shifted
2607 : //of the difference among primary vertices (vTrasl) and
2608 : //the covariance matrix is changed accordingly
2609 : //(covV = covariance of the primary vertex).
2610 : //Origin: "Romita, Rossella" <R.Romita@gsi.de>
2611 : //
2612 0 : TVector3 translation;
2613 : // vTrasl coordinates in the local system
2614 0 : translation.SetXYZ(vTrasl[0],vTrasl[1],vTrasl[2]);
2615 0 : translation.RotateZ(-fAlpha);
2616 0 : translation.GetXYZ(vTrasl);
2617 :
2618 : //compute the new x,y,z of the track
2619 0 : Double_t newX=fX-vTrasl[0];
2620 0 : Double_t newY=fP[0]-vTrasl[1];
2621 0 : Double_t newZ=fP[1]-vTrasl[2];
2622 :
2623 : //define the new parameters
2624 0 : Double_t newParam[5];
2625 0 : newParam[0]=newY;
2626 0 : newParam[1]=newZ;
2627 0 : newParam[2]=fP[2];
2628 0 : newParam[3]=fP[3];
2629 0 : newParam[4]=fP[4];
2630 :
2631 : // recompute the covariance matrix:
2632 : // 1. covV in the local system
2633 0 : Double_t cosRot=TMath::Cos(fAlpha), sinRot=TMath::Sin(fAlpha);
2634 0 : TMatrixD qQi(3,3);
2635 0 : qQi(0,0) = cosRot;
2636 0 : qQi(0,1) = sinRot;
2637 0 : qQi(0,2) = 0.;
2638 0 : qQi(1,0) = -sinRot;
2639 0 : qQi(1,1) = cosRot;
2640 0 : qQi(1,2) = 0.;
2641 0 : qQi(2,0) = 0.;
2642 0 : qQi(2,1) = 0.;
2643 0 : qQi(2,2) = 1.;
2644 0 : TMatrixD uUi(3,3);
2645 0 : uUi(0,0) = covV[0];
2646 0 : uUi(0,0) = covV[0];
2647 0 : uUi(1,0) = covV[1];
2648 0 : uUi(0,1) = covV[1];
2649 0 : uUi(2,0) = covV[3];
2650 0 : uUi(0,2) = covV[3];
2651 0 : uUi(1,1) = covV[2];
2652 0 : uUi(2,2) = covV[5];
2653 0 : uUi(1,2) = covV[4];
2654 0 : if(uUi.Determinant() <= 0.) {return kFALSE;}
2655 0 : TMatrixD uUiQi(uUi,TMatrixD::kMult,qQi);
2656 0 : TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiQi);
2657 :
2658 : //2. compute the new covariance matrix of the track
2659 0 : Double_t sigmaXX=m(0,0);
2660 0 : Double_t sigmaXZ=m(2,0);
2661 0 : Double_t sigmaXY=m(1,0);
2662 0 : Double_t sigmaYY=GetSigmaY2()+m(1,1);
2663 0 : Double_t sigmaYZ=fC[1]+m(1,2);
2664 0 : Double_t sigmaZZ=fC[2]+m(2,2);
2665 0 : Double_t covarianceYY=sigmaYY + (-1.)*((sigmaXY*sigmaXY)/sigmaXX);
2666 0 : Double_t covarianceYZ=sigmaYZ-(sigmaXZ*sigmaXY/sigmaXX);
2667 0 : Double_t covarianceZZ=sigmaZZ-((sigmaXZ*sigmaXZ)/sigmaXX);
2668 :
2669 0 : Double_t newCov[15];
2670 0 : newCov[0]=covarianceYY;
2671 0 : newCov[1]=covarianceYZ;
2672 0 : newCov[2]=covarianceZZ;
2673 0 : for(Int_t i=3;i<15;i++){
2674 0 : newCov[i]=fC[i];
2675 : }
2676 :
2677 : // set the new parameters
2678 :
2679 0 : Set(newX,fAlpha,newParam,newCov);
2680 :
2681 : return kTRUE;
2682 0 : }
2683 :
2684 : void AliExternalTrackParam::CheckCovariance() {
2685 :
2686 : // This function forces the diagonal elements of the covariance matrix to be positive.
2687 : // In case the diagonal element is bigger than the maximal allowed value, it is set to
2688 : // the limit and the off-diagonal elements that correspond to it are set to zero.
2689 :
2690 1539102 : fC[0] = TMath::Abs(fC[0]);
2691 769551 : if (fC[0]>kC0max) {
2692 18 : double scl = TMath::Sqrt(kC0max/fC[0]);
2693 18 : fC[0] = kC0max;
2694 18 : fC[1] *= scl;
2695 18 : fC[3] *= scl;
2696 18 : fC[6] *= scl;
2697 18 : fC[10] *= scl;
2698 18 : }
2699 769551 : fC[2] = TMath::Abs(fC[2]);
2700 769551 : if (fC[2]>kC2max) {
2701 4 : double scl = TMath::Sqrt(kC2max/fC[2]);
2702 4 : fC[2] = kC2max;
2703 4 : fC[1] *= scl;
2704 4 : fC[4] *= scl;
2705 4 : fC[7] *= scl;
2706 4 : fC[11] *= scl;
2707 4 : }
2708 769551 : fC[5] = TMath::Abs(fC[5]);
2709 769551 : if (fC[5]>kC5max) {
2710 6 : double scl = TMath::Sqrt(kC5max/fC[5]);
2711 6 : fC[5] = kC5max;
2712 6 : fC[3] *= scl;
2713 6 : fC[4] *= scl;
2714 6 : fC[8] *= scl;
2715 6 : fC[12] *= scl;
2716 6 : }
2717 769551 : fC[9] = TMath::Abs(fC[9]);
2718 769551 : if (fC[9]>kC9max) {
2719 0 : double scl = TMath::Sqrt(kC9max/fC[9]);
2720 0 : fC[9] = kC9max;
2721 0 : fC[6] *= scl;
2722 0 : fC[7] *= scl;
2723 0 : fC[8] *= scl;
2724 0 : fC[13] *= scl;
2725 0 : }
2726 769551 : fC[14] = TMath::Abs(fC[14]);
2727 769551 : if (fC[14]>kC14max) {
2728 0 : double scl = TMath::Sqrt(kC14max/fC[14]);
2729 0 : fC[14] = kC14max;
2730 0 : fC[10] *= scl;
2731 0 : fC[11] *= scl;
2732 0 : fC[12] *= scl;
2733 0 : fC[13] *= scl;
2734 0 : }
2735 :
2736 : // The part below is used for tests and normally is commented out
2737 : // TMatrixDSym m(5);
2738 : // TVectorD eig(5);
2739 :
2740 : // m(0,0)=fC[0];
2741 : // m(1,0)=fC[1]; m(1,1)=fC[2];
2742 : // m(2,0)=fC[3]; m(2,1)=fC[4]; m(2,2)=fC[5];
2743 : // m(3,0)=fC[6]; m(3,1)=fC[7]; m(3,2)=fC[8]; m(3,3)=fC[9];
2744 : // m(4,0)=fC[10]; m(4,1)=fC[11]; m(4,2)=fC[12]; m(4,3)=fC[13]; m(4,4)=fC[14];
2745 :
2746 : // m(0,1)=m(1,0);
2747 : // m(0,2)=m(2,0); m(1,2)=m(2,1);
2748 : // m(0,3)=m(3,0); m(1,3)=m(3,1); m(2,3)=m(3,2);
2749 : // m(0,4)=m(4,0); m(1,4)=m(4,1); m(2,4)=m(4,2); m(3,4)=m(4,3);
2750 : // m.EigenVectors(eig);
2751 :
2752 : // // assert(eig(0)>=0 && eig(1)>=0 && eig(2)>=0 && eig(3)>=0 && eig(4)>=0);
2753 : // if (!(eig(0)>=0 && eig(1)>=0 && eig(2)>=0 && eig(3)>=0 && eig(4)>=0)) {
2754 : // AliWarning("Negative eigenvalues of the covariance matrix!");
2755 : // this->Print();
2756 : // eig.Print();
2757 : // }
2758 769551 : }
2759 :
2760 : Bool_t AliExternalTrackParam::ConstrainToVertex(const AliVVertex* vtx, Double_t b[3])
2761 : {
2762 : // Constrain TPC inner params constrained
2763 : //
2764 0 : if (!vtx)
2765 0 : return kFALSE;
2766 :
2767 0 : Double_t dz[2], cov[3];
2768 0 : if (!PropagateToDCABxByBz(vtx, b, 3, dz, cov))
2769 0 : return kFALSE;
2770 :
2771 0 : Double_t covar[6];
2772 0 : vtx->GetCovarianceMatrix(covar);
2773 :
2774 0 : Double_t p[2]= { fP[0] - dz[0], fP[1] - dz[1] };
2775 0 : Double_t c[3]= { covar[2], 0., covar[5] };
2776 :
2777 0 : Double_t chi2C = GetPredictedChi2(p,c);
2778 0 : if (chi2C>kVeryBig)
2779 0 : return kFALSE;
2780 :
2781 0 : if (!Update(p,c))
2782 0 : return kFALSE;
2783 :
2784 0 : return kTRUE;
2785 0 : }
2786 :
2787 : //___________________________________________________________________________________________
2788 : Bool_t AliExternalTrackParam::GetXatLabR(Double_t r,Double_t &x, Double_t bz, Int_t dir) const
2789 : {
2790 : // Get local X of the track position estimated at the radius lab radius r.
2791 : // The track curvature is accounted exactly
2792 : //
2793 : // The flag "dir" can be used to remove the ambiguity of which intersection to take (out of 2 possible)
2794 : // 0 - take the intersection closest to the current track position
2795 : // >0 - go along the track (increasing fX)
2796 : // <0 - go backward (decreasing fX)
2797 : //
2798 0 : const Double_t &fy=fP[0], &sn = fP[2];
2799 : const double kEps = 1.e-6;
2800 : //
2801 0 : double crv = GetC(bz);
2802 0 : if (TMath::Abs(crv)>kAlmost0) { // helix
2803 : // get center of the track circle
2804 0 : double tR = 1./crv; // track radius (for the moment signed)
2805 0 : double cs = TMath::Sqrt((1-sn)*(1+sn));
2806 0 : double x0 = fX - sn*tR;
2807 0 : double y0 = fy + cs*tR;
2808 0 : double r0 = TMath::Sqrt(x0*x0+y0*y0);
2809 : // printf("Xc:%+e Yc:%+e tR:%e r0:%e\n",x0,y0,tR,r0);
2810 : //
2811 0 : if (r0<=kAlmost0) return kFALSE; // the track is concentric to circle
2812 0 : tR = TMath::Abs(tR);
2813 : double tR2r0=1.,g=0,tmp=0;
2814 0 : if (TMath::Abs(tR-r0)>kEps) {
2815 0 : tR2r0 = tR/r0;
2816 0 : g = 0.5*(r*r/(r0*tR) - tR2r0 - 1./tR2r0);
2817 0 : tmp = 1.+g*tR2r0;
2818 0 : }
2819 : else {
2820 : tR2r0 = 1.0;
2821 0 : g = 0.5*r*r/(r0*tR) - 1;
2822 0 : tmp = 0.5*r*r/(r0*r0);
2823 : }
2824 0 : double det = (1.-g)*(1.+g);
2825 0 : if (det<0) return kFALSE; // does not reach raduis r
2826 0 : det = TMath::Sqrt(det);
2827 : //
2828 : // the intersection happens in 2 points: {x0+tR*C,y0+tR*S}
2829 : // with C=f*c0+-|s0|*det and S=f*s0-+c0 sign(s0)*det
2830 : // where s0 and c0 make direction for the circle center (=x0/r0 and y0/r0)
2831 : //
2832 0 : x = x0*tmp;
2833 0 : double y = y0*tmp;
2834 0 : if (TMath::Abs(y0)>kAlmost0) { // when y0==0 the x,y is unique
2835 0 : double dfx = tR2r0*TMath::Abs(y0)*det;
2836 0 : double dfy = tR2r0*x0*TMath::Sign(det,y0);
2837 0 : if (dir==0) { // chose the one which corresponds to smallest step
2838 0 : double delta = (x-fX)*dfx-(y-fy)*dfy; // the choice of + in C will lead to smaller step if delta<0
2839 0 : if (delta<0) x += dfx;
2840 0 : else x -= dfx;
2841 0 : }
2842 0 : else if (dir>0) { // along track direction: x must be > fX
2843 0 : x -= dfx; // try the smallest step (dfx is positive)
2844 0 : double dfeps = fX-x; // handle special case of very small step
2845 0 : if (dfeps<-kEps) return kTRUE;
2846 0 : if (TMath::Abs(dfeps)<kEps && // are we already in right r?
2847 0 : TMath::Abs(fX*fX+fy*fy - r*r)<kEps) return fX;
2848 0 : x += dfx+dfx;
2849 0 : if (x-fX>0) return kTRUE;
2850 0 : if (x-fX<-kEps) return kFALSE;
2851 0 : x = fX; // don't move
2852 0 : }
2853 : else { // backward: x must be < fX
2854 0 : x += dfx; // try the smallest step (dfx is positive)
2855 0 : double dfeps = x-fX; // handle special case of very small step
2856 0 : if (dfeps<-kEps) return kTRUE;
2857 0 : if (TMath::Abs(dfeps)<kEps && // are we already in right r?
2858 0 : TMath::Abs(fX*fX+fy*fy - r*r)<kEps) return fX;
2859 0 : x-=dfx+dfx;
2860 0 : if (x-fX<0) return kTRUE;
2861 0 : if (x-fX>kEps) return kFALSE;
2862 0 : x = fX; // don't move
2863 0 : }
2864 0 : }
2865 : else { // special case: track touching the circle just in 1 point
2866 0 : if ( (dir>0&&x<fX) || (dir<0&&x>fX) ) return kFALSE;
2867 : }
2868 0 : }
2869 : else { // this is a straight track
2870 0 : if (TMath::Abs(sn)>=kAlmost1) { // || to Y axis
2871 0 : double det = (r-fX)*(r+fX);
2872 0 : if (det<0) return kFALSE; // does not reach raduis r
2873 0 : x = fX;
2874 0 : if (dir==0) return kTRUE;
2875 0 : det = TMath::Sqrt(det);
2876 0 : if (dir>0) { // along the track direction
2877 0 : if (sn>0) {if (fy>det) return kFALSE;} // track is along Y axis and above the circle
2878 0 : else {if (fy<-det) return kFALSE;} // track is against Y axis amd belo the circle
2879 : }
2880 0 : else if(dir>0) { // agains track direction
2881 0 : if (sn>0) {if (fy<-det) return kFALSE;} // track is along Y axis
2882 0 : else if (fy>det) return kFALSE; // track is against Y axis
2883 : }
2884 0 : }
2885 0 : else if (TMath::Abs(sn)<=kAlmost0) { // || to X axis
2886 0 : double det = (r-fy)*(r+fy);
2887 0 : if (det<0) return kFALSE; // does not reach raduis r
2888 0 : det = TMath::Sqrt(det);
2889 0 : if (!dir) {
2890 0 : x = fX>0 ? det : -det; // choose the solution requiring the smalest step
2891 0 : return kTRUE;
2892 : }
2893 0 : else if (dir>0) { // along the track direction
2894 0 : if (fX > det) return kFALSE; // current point is in on the right from the circle
2895 0 : else if (fX <-det) x = -det; // on the left
2896 0 : else x = det; // within the circle
2897 : }
2898 : else { // against the track direction
2899 0 : if (fX <-det) return kFALSE;
2900 0 : else if (fX > det) x = det;
2901 0 : else x = -det;
2902 : }
2903 0 : }
2904 : else { // general case of straight line
2905 0 : double cs = TMath::Sqrt((1-sn)*(1+sn));
2906 0 : double xsyc = fX*sn-fy*cs;
2907 0 : double det = (r-xsyc)*(r+xsyc);
2908 0 : if (det<0) return kFALSE; // does not reach raduis r
2909 0 : det = TMath::Sqrt(det);
2910 0 : double xcys = fX*cs+fy*sn;
2911 0 : double t = -xcys;
2912 0 : if (dir==0) t += t>0 ? -det:det; // chose the solution requiring the smalest step
2913 0 : else if (dir>0) { // go in increasing fX direction. ( t+-det > 0)
2914 0 : if (t>=-det) t += -det; // take minimal step giving t>0
2915 0 : else return kFALSE; // both solutions have negative t
2916 0 : }
2917 : else { // go in increasing fX direction. (t+-det < 0)
2918 0 : if (t<det) t -= det; // take minimal step giving t<0
2919 0 : else return kFALSE; // both solutions have positive t
2920 : }
2921 0 : x = fX + cs*t;
2922 0 : }
2923 : }
2924 : //
2925 0 : return kTRUE;
2926 0 : }
2927 : //_________________________________________________________
2928 : Bool_t AliExternalTrackParam::GetXYZatR(Double_t xr,Double_t bz, Double_t *xyz, Double_t* alpSect) const
2929 : {
2930 : // This method has 3 modes of behaviour
2931 : // 1) xyz[3] array is provided but alpSect pointer is 0: calculate the position of track intersection
2932 : // with circle of radius xr and fill it in xyz array
2933 : // 2) alpSect pointer is provided: find alpha of the sector where the track reaches local coordinate xr
2934 : // Note that in this case xr is NOT the radius but the local coordinate.
2935 : // If the xyz array is provided, it will be filled by track lab coordinates at local X in this sector
2936 : // 3) Neither alpSect nor xyz pointers are provided: just check if the track reaches radius xr
2937 : //
2938 : //
2939 0 : double crv = GetC(bz);
2940 0 : if ( (TMath::Abs(bz))<kAlmost0Field ) crv=0.;
2941 0 : const double &fy = fP[0];
2942 0 : const double &fz = fP[1];
2943 0 : const double &sn = fP[2];
2944 0 : const double &tgl = fP[3];
2945 : //
2946 : // general circle parameterization:
2947 : // x = (r0+tR)cos(phi0) - tR cos(t+phi0)
2948 : // y = (r0+tR)sin(phi0) - tR sin(t+phi0)
2949 : // where qb is the sign of the curvature, tR is the track's signed radius and r0
2950 : // is the DCA of helix to origin
2951 : //
2952 0 : double tR = 1./crv; // track radius signed
2953 0 : double cs = TMath::Sqrt((1-sn)*(1+sn));
2954 0 : double x0 = fX - sn*tR; // helix center coordinates
2955 0 : double y0 = fy + cs*tR;
2956 0 : double phi0 = TMath::ATan2(y0,x0); // angle of PCA wrt to the origin
2957 0 : if (tR<0) phi0 += TMath::Pi();
2958 0 : if (phi0 > TMath::Pi()) phi0 -= 2.*TMath::Pi();
2959 0 : else if (phi0 <-TMath::Pi()) phi0 += 2.*TMath::Pi();
2960 0 : double cs0 = TMath::Cos(phi0);
2961 0 : double sn0 = TMath::Sin(phi0);
2962 0 : double r0 = x0*cs0 + y0*sn0 - tR; // DCA to origin
2963 0 : double r2R = 1.+r0/tR;
2964 : //
2965 : //
2966 0 : if (r2R<kAlmost0) return kFALSE; // helix is centered at the origin, no specific intersection with other concetric circle
2967 0 : if (!xyz && !alpSect) return kTRUE;
2968 0 : double xr2R = xr/tR;
2969 0 : double r2Ri = 1./r2R;
2970 : // the intersection cos(t) = [1 + (r0/tR+1)^2 - (r0/tR)^2]/[2(1+r0/tR)]
2971 0 : double cosT = 0.5*(r2R + (1-xr2R*xr2R)*r2Ri);
2972 0 : if ( TMath::Abs(cosT)>kAlmost1 ) {
2973 : // printf("Does not reach : %f %f\n",r0,tR);
2974 0 : return kFALSE; // track does not reach the radius xr
2975 : }
2976 : //
2977 0 : double t = TMath::ACos(cosT);
2978 0 : if (tR<0) t = -t;
2979 : // intersection point
2980 0 : double xyzi[3];
2981 0 : xyzi[0] = x0 - tR*TMath::Cos(t+phi0);
2982 0 : xyzi[1] = y0 - tR*TMath::Sin(t+phi0);
2983 0 : if (xyz) { // if postition is requested, then z is needed:
2984 0 : double t0 = TMath::ATan2(cs,-sn) - phi0;
2985 0 : double z0 = fz - t0*tR*tgl;
2986 0 : xyzi[2] = z0 + tR*t*tgl;
2987 0 : }
2988 0 : else xyzi[2] = 0;
2989 : //
2990 0 : Local2GlobalPosition(xyzi,fAlpha);
2991 : //
2992 0 : if (xyz) {
2993 0 : xyz[0] = xyzi[0];
2994 0 : xyz[1] = xyzi[1];
2995 0 : xyz[2] = xyzi[2];
2996 0 : }
2997 : //
2998 0 : if (alpSect) {
2999 : double &alp = *alpSect;
3000 : // determine the sector of crossing
3001 0 : double phiPos = TMath::Pi()+TMath::ATan2(-xyzi[1],-xyzi[0]);
3002 0 : int sect = ((Int_t)(phiPos*TMath::RadToDeg()))/20;
3003 0 : alp = TMath::DegToRad()*(20*sect+10);
3004 0 : double x2r,f1,f2,r1,r2,dx,dy2dx,yloc=0, ylocMax = xr*TMath::Tan(TMath::Pi()/18); // min max Y within sector at given X
3005 : //
3006 0 : while(1) {
3007 0 : Double_t ca=TMath::Cos(alp-fAlpha), sa=TMath::Sin(alp-fAlpha);
3008 0 : if ((cs*ca+sn*sa)<0) {
3009 0 : AliDebug(1,Form("Rotation to target sector impossible: local cos(phi) would become %.2f",cs*ca+sn*sa));
3010 0 : return kFALSE;
3011 : }
3012 : //
3013 0 : f1 = sn*ca - cs*sa;
3014 0 : if (TMath::Abs(f1) >= kAlmost1) {
3015 0 : AliDebug(1,Form("Rotation to target sector impossible: local sin(phi) would become %.2f",f1));
3016 0 : return kFALSE;
3017 : }
3018 : //
3019 0 : double tmpX = fX*ca + fy*sa;
3020 0 : double tmpY = -fX*sa + fy*ca;
3021 : //
3022 : // estimate Y at X=xr
3023 0 : dx=xr-tmpX;
3024 0 : x2r = crv*dx;
3025 0 : f2=f1 + x2r;
3026 0 : if (TMath::Abs(f2) >= kAlmost1) {
3027 0 : AliDebug(1,Form("Propagation in target sector failed ! %.10e",f2));
3028 0 : return kFALSE;
3029 : }
3030 0 : r1 = TMath::Sqrt((1.-f1)*(1.+f1));
3031 0 : r2 = TMath::Sqrt((1.-f2)*(1.+f2));
3032 0 : dy2dx = (f1+f2)/(r1+r2);
3033 0 : yloc = tmpY + dx*dy2dx;
3034 0 : if (yloc>ylocMax) {alp += 2*TMath::Pi()/18; sect++;}
3035 0 : else if (yloc<-ylocMax) {alp -= 2*TMath::Pi()/18; sect--;}
3036 0 : else break;
3037 0 : if (alp >= TMath::Pi()) alp -= 2*TMath::Pi();
3038 0 : else if (alp < -TMath::Pi()) alp += 2*TMath::Pi();
3039 : // if (sect>=18) sect = 0;
3040 : // if (sect<=0) sect = 17;
3041 0 : }
3042 : //
3043 : // if alpha was requested, then recalculate the position at intersection in sector
3044 0 : if (xyz) {
3045 0 : xyz[0] = xr;
3046 0 : xyz[1] = yloc;
3047 0 : if (TMath::Abs(x2r)<0.05) xyz[2] = fz + dx*(r2 + f2*dy2dx)*tgl;
3048 : else {
3049 : // for small dx/R the linear apporximation of the arc by the segment is OK,
3050 : // but at large dx/R the error is very large and leads to incorrect Z propagation
3051 : // angle traversed delta = 2*asin(dist_start_end / R / 2), hence the arc is: R*deltaPhi
3052 : // The dist_start_end is obtained from sqrt(dx^2+dy^2) = x/(r1+r2)*sqrt(2+f1*f2+r1*r2)
3053 : // Similarly, the rotation angle in linear in dx only for dx<<R
3054 0 : double chord = dx*TMath::Sqrt(1+dy2dx*dy2dx); // distance from old position to new one
3055 0 : double rot = 2*TMath::ASin(0.5*chord*crv); // angular difference seen from the circle center
3056 0 : xyz[2] = fz + rot/crv*tgl;
3057 : }
3058 0 : Local2GlobalPosition(xyz,alp);
3059 0 : }
3060 0 : }
3061 0 : return kTRUE;
3062 : //
3063 0 : }
3064 :
3065 :
3066 : Double_t AliExternalTrackParam::GetParameterAtRadius(Double_t r, Double_t bz, Int_t parType) const
3067 : {
3068 : //
3069 : // Get track parameters at the radius of interest.
3070 : // Given function is aimed to be used to interactivelly (tree->Draw())
3071 : // access track properties at different radii
3072 : //
3073 : // TO BE USED WITH SPECICAL CARE -
3074 : // it is aimed to be used for rough calculation as constant field and
3075 : // no correction for material is used
3076 : //
3077 : // r - radius of interest
3078 : // bz - magentic field
3079 : // retun values dependens on parType:
3080 : // parType = 0 -gx
3081 : // parType = 1 -gy
3082 : // parType = 2 -gz
3083 : //
3084 : // parType = 3 -pgx
3085 : // parType = 4 -pgy
3086 : // parType = 5 -pgz
3087 : //
3088 : // parType = 6 - r
3089 : // parType = 7 - global position phi
3090 : // parType = 8 - global direction phi
3091 : // parType = 9 - direction phi- positionphi
3092 0 : if (parType<0) {
3093 : parType=-1;
3094 0 : return 0;
3095 : }
3096 0 : Double_t xyz[3];
3097 0 : Double_t pxyz[3];
3098 0 : Double_t localX=0;
3099 0 : Bool_t res = GetXatLabR(r,localX,bz,1);
3100 0 : if (!res) {
3101 : parType=-1;
3102 0 : return 0;
3103 : }
3104 : //
3105 : // position parameters
3106 : //
3107 0 : GetXYZAt(localX,bz,xyz);
3108 0 : if (parType<3) {
3109 0 : return xyz[parType];
3110 : }
3111 :
3112 0 : if (parType==6) return TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
3113 0 : if (parType==7) return TMath::ATan2(xyz[1],xyz[0]);
3114 : //
3115 : // momenta parameters
3116 : //
3117 0 : GetPxPyPzAt(localX,bz,pxyz);
3118 0 : if (parType==8) return TMath::ATan2(pxyz[1],pxyz[0]);
3119 0 : if (parType==9) {
3120 0 : Double_t diff = TMath::ATan2(pxyz[1],pxyz[0])-TMath::ATan2(xyz[1],xyz[0]);
3121 0 : if (diff>TMath::Pi()) diff-=TMath::TwoPi();
3122 0 : if (diff<-TMath::Pi()) diff+=TMath::TwoPi();
3123 : return diff;
3124 : }
3125 0 : if (parType>=3&&parType<6) {
3126 0 : return pxyz[parType%3];
3127 : }
3128 0 : return 0;
3129 0 : }
|
