Moving the classes that belong to the following libraries: STEERBase, ESD, CDB, AOD...
[u/mrichter/AliRoot.git] / STEER / ESD / AliKFParticle.cxx
CommitLineData
f826d409 1//----------------------------------------------------------------------------
2// Implementation of the AliKFParticle class
3// .
4// @author S.Gorbunov, I.Kisel
5// @version 1.0
6// @since 13.05.07
7//
8// Class to reconstruct and store the decayed particle parameters.
9// The method is described in CBM-SOFT note 2007-003,
10// ``Reconstruction of decayed particles based on the Kalman filter'',
11// http://www.gsi.de/documents/DOC-2007-May-14-1.pdf
12//
13// This class is ALICE interface to general mathematics in AliKFParticleCore
14//
15// -= Copyright &copy ALICE HLT Group =-
16//____________________________________________________________________________
17
18
19#include "AliKFParticle.h"
20#include "TDatabasePDG.h"
21#include "TParticlePDG.h"
706952f5 22#include "AliVTrack.h"
cfb5a2e1 23#include "AliVVertex.h"
f826d409 24
effa6338 25ClassImp(AliKFParticle)
f826d409 26
616ffc76 27Double_t AliKFParticle::fgBz = -5.; //* Bz compoment of the magnetic field
f826d409 28
537b06a4 29AliKFParticle::AliKFParticle( const AliKFParticle &d1, const AliKFParticle &d2, Bool_t gamma )
30{
31 if (!gamma) {
32 AliKFParticle mother;
33 mother+= d1;
34 mother+= d2;
35 *this = mother;
36 } else
37 ConstructGamma(d1, d2);
38}
39
7972a8db 40void AliKFParticle::Create( const Double_t Param[], const Double_t Cov[], Int_t Charge, Int_t PID )
f826d409 41{
e7b09c95 42 // Constructor from "cartesian" track, PID hypothesis should be provided
43 //
44 // Param[6] = { X, Y, Z, Px, Py, Pz } - position and momentum
45 // Cov [21] = lower-triangular part of the covariance matrix:
46 //
47 // ( 0 . . . . . )
48 // ( 1 2 . . . . )
49 // Cov. matrix = ( 3 4 5 . . . ) - numbering of covariance elements in Cov[]
50 // ( 6 7 8 9 . . )
51 // ( 10 11 12 13 14 . )
52 // ( 15 16 17 18 19 20 )
706952f5 53 Double_t C[21];
54 for( int i=0; i<21; i++ ) C[i] = Cov[i];
55
f826d409 56 TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(PID);
57 Double_t mass = (particlePDG) ? particlePDG->Mass() :0.13957;
e7b09c95 58
706952f5 59 AliKFParticleBase::Initialize( Param, C, Charge, mass );
e7b09c95 60}
616ffc76 61
706952f5 62AliKFParticle::AliKFParticle( const AliVTrack &track, Int_t PID )
e7b09c95 63{
64 // Constructor from ALICE track, PID hypothesis should be provided
616ffc76 65
b63c3605 66 track.GetXYZ(fP);
706952f5 67 track.PxPyPz(fP+3);
68 fQ = track.Charge();
f826d409 69 track.GetCovarianceXYZPxPyPz( fC );
7972a8db 70 Create(fP,fC,fQ,PID);
f826d409 71}
72
706952f5 73AliKFParticle::AliKFParticle( const AliVVertex &vertex )
f826d409 74{
75 // Constructor from ALICE vertex
76
77 vertex.GetXYZ( fP );
706952f5 78 vertex.GetCovarianceMatrix( fC );
f826d409 79 fChi2 = vertex.GetChi2();
80 fNDF = 2*vertex.GetNContributors() - 3;
81 fQ = 0;
82 fAtProductionVertex = 0;
83 fIsLinearized = 0;
84 fSFromDecay = 0;
85}
616ffc76 86
616ffc76 87void AliKFParticle::GetExternalTrackParam( const AliKFParticleBase &p, Double_t &X, Double_t &Alpha, Double_t P[5] )
88{
e7b09c95 89 // Conversion to AliExternalTrackParam parameterization
90
616ffc76 91 Double_t cosA = p.GetPx(), sinA = p.GetPy();
92 Double_t pt = TMath::Sqrt(cosA*cosA + sinA*sinA);
93 Double_t pti = 0;
94 if( pt<1.e-4 ){
95 cosA = 1;
96 sinA = 0;
97 } else {
98 pti = 1./pt;
99 cosA*=pti;
100 sinA*=pti;
101 }
102 Alpha = TMath::ATan2(sinA,cosA);
103 X = p.GetX()*cosA + p.GetY()*sinA;
104 P[0]= p.GetY()*cosA - p.GetX()*sinA;
105 P[1]= p.GetZ();
106 P[2]= 0;
107 P[3]= p.GetPz()*pti;
108 P[4]= p.GetQ()*pti;
109}
110
446ce366 111Bool_t AliKFParticle::GetDistanceFromVertexXY( const Double_t vtx[], const Double_t Cv[], Double_t &val, Double_t &err ) const
112{
113 //* Calculate DCA distance from vertex (transverse impact parameter) in XY
114 //* v = [xy], Cv=[Cxx,Cxy,Cyy ]-covariance matrix
115
116 Bool_t ret = 0;
117
118 Double_t mP[8];
119 Double_t mC[36];
120
121 Transport( GetDStoPoint(vtx), mP, mC );
122
123 Double_t dx = mP[0] - vtx[0];
124 Double_t dy = mP[1] - vtx[1];
125 Double_t px = mP[3];
126 Double_t py = mP[4];
127 Double_t pt = TMath::Sqrt(px*px + py*py);
128 Double_t ex=0, ey=0;
129 if( pt<1.e-4 ){
130 ret = 1;
131 pt = 1.;
132 val = 1.e4;
133 } else{
134 ex = px/pt;
924a4a7b 135 ey = py/pt;
446ce366 136 val = dy*ex - dx*ey;
137 }
616ffc76 138
446ce366 139 Double_t h0 = -ey;
140 Double_t h1 = ex;
141 Double_t h3 = (dy*ey + dx*ex)*ey/pt;
142 Double_t h4 = -(dy*ey + dx*ex)*ex/pt;
143
144 err =
145 h0*(h0*GetCovariance(0,0) + h1*GetCovariance(0,1) + h3*GetCovariance(0,3) + h4*GetCovariance(0,4) ) +
146 h1*(h0*GetCovariance(1,0) + h1*GetCovariance(1,1) + h3*GetCovariance(1,3) + h4*GetCovariance(1,4) ) +
147 h3*(h0*GetCovariance(3,0) + h1*GetCovariance(3,1) + h3*GetCovariance(3,3) + h4*GetCovariance(3,4) ) +
148 h4*(h0*GetCovariance(4,0) + h1*GetCovariance(4,1) + h3*GetCovariance(4,3) + h4*GetCovariance(4,4) );
616ffc76 149
446ce366 150 if( Cv ){
151 err+= h0*(h0*Cv[0] + h1*Cv[1] ) + h1*(h0*Cv[1] + h1*Cv[2] );
152 }
153
154 err = TMath::Sqrt(TMath::Abs(err));
155
156 return ret;
157}
158
159Bool_t AliKFParticle::GetDistanceFromVertexXY( const Double_t vtx[], Double_t &val, Double_t &err ) const
160{
161 return GetDistanceFromVertexXY( vtx, 0, val, err );
162}
163
164
165Bool_t AliKFParticle::GetDistanceFromVertexXY( const AliKFParticle &Vtx, Double_t &val, Double_t &err ) const
e7b09c95 166{
167 //* Calculate distance from vertex [cm] in XY-plane
168
446ce366 169 return GetDistanceFromVertexXY( Vtx.fP, Vtx.fC, val, err );
170}
171
172Bool_t AliKFParticle::GetDistanceFromVertexXY( const AliVVertex &Vtx, Double_t &val, Double_t &err ) const
173{
174 //* Calculate distance from vertex [cm] in XY-plane
175
176 return GetDistanceFromVertexXY( AliKFParticle(Vtx), val, err );
177}
178
179Double_t AliKFParticle::GetDistanceFromVertexXY( const Double_t vtx[] ) const
180{
181 //* Calculate distance from vertex [cm] in XY-plane
182 Double_t val, err;
183 GetDistanceFromVertexXY( vtx, 0, val, err );
184 return val;
e7b09c95 185}
186
187Double_t AliKFParticle::GetDistanceFromVertexXY( const AliKFParticle &Vtx ) const
188{
189 //* Calculate distance from vertex [cm] in XY-plane
190
191 return GetDistanceFromVertexXY( Vtx.fP );
192}
193
706952f5 194Double_t AliKFParticle::GetDistanceFromVertexXY( const AliVVertex &Vtx ) const
e7b09c95 195{
196 //* Calculate distance from vertex [cm] in XY-plane
197
198 return GetDistanceFromVertexXY( AliKFParticle(Vtx).fP );
199}
200
201Double_t AliKFParticle::GetDistanceFromParticleXY( const AliKFParticle &p ) const
202{
203 //* Calculate distance to other particle [cm]
204
205 Double_t dS, dS1;
206 GetDStoParticleXY( p, dS, dS1 );
207 Double_t mP[8], mC[36], mP1[8], mC1[36];
208 Transport( dS, mP, mC );
209 p.Transport( dS1, mP1, mC1 );
210 Double_t dx = mP[0]-mP1[0];
211 Double_t dy = mP[1]-mP1[1];
212 return TMath::Sqrt(dx*dx+dy*dy);
213}
214
215Double_t AliKFParticle::GetDeviationFromParticleXY( const AliKFParticle &p ) const
216{
217 //* Calculate sqrt(Chi2/ndf) deviation from other particle
218
219 Double_t dS, dS1;
220 GetDStoParticleXY( p, dS, dS1 );
221 Double_t mP1[8], mC1[36];
222 p.Transport( dS1, mP1, mC1 );
616ffc76 223
e7b09c95 224 Double_t d[2]={ fP[0]-mP1[0], fP[1]-mP1[1] };
616ffc76 225
e7b09c95 226 Double_t sigmaS = .1+10.*TMath::Sqrt( (d[0]*d[0]+d[1]*d[1] )/
227 (mP1[3]*mP1[3]+mP1[4]*mP1[4] ) );
616ffc76 228
e7b09c95 229 Double_t h[2] = { mP1[3]*sigmaS, mP1[4]*sigmaS };
616ffc76 230
e7b09c95 231 mC1[0] +=h[0]*h[0];
232 mC1[1] +=h[1]*h[0];
233 mC1[2] +=h[1]*h[1];
234
235 return GetDeviationFromVertexXY( mP1, mC1 )*TMath::Sqrt(2./1.);
236}
237
238
446ce366 239Double_t AliKFParticle::GetDeviationFromVertexXY( const Double_t vtx[], const Double_t Cv[] ) const
e7b09c95 240{
241 //* Calculate sqrt(Chi2/ndf) deviation from vertex
242 //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
243
446ce366 244 Double_t val, err;
245 Bool_t problem = GetDistanceFromVertexXY( vtx, Cv, val, err );
246 if( problem || err<1.e-20 ) return 1.e4;
247 else return val/err;
e7b09c95 248}
249
250
706952f5 251Double_t AliKFParticle::GetDeviationFromVertexXY( const AliKFParticle &Vtx ) const
e7b09c95 252{
253 //* Calculate sqrt(Chi2/ndf) deviation from vertex
254 //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
255
256 return GetDeviationFromVertexXY( Vtx.fP, Vtx.fC );
257}
258
706952f5 259Double_t AliKFParticle::GetDeviationFromVertexXY( const AliVVertex &Vtx ) const
e7b09c95 260{
261 //* Calculate sqrt(Chi2/ndf) deviation from vertex
262 //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
263
264 AliKFParticle v(Vtx);
265 return GetDeviationFromVertexXY( v.fP, v.fC );
266}
267
e7b09c95 268Double_t AliKFParticle::GetAngle ( const AliKFParticle &p ) const
269{
270 //* Calculate the opening angle between two particles
271
272 Double_t dS, dS1;
273 GetDStoParticle( p, dS, dS1 );
274 Double_t mP[8], mC[36], mP1[8], mC1[36];
275 Transport( dS, mP, mC );
276 p.Transport( dS1, mP1, mC1 );
277 Double_t n = TMath::Sqrt( mP[3]*mP[3] + mP[4]*mP[4] + mP[5]*mP[5] );
278 Double_t n1= TMath::Sqrt( mP1[3]*mP1[3] + mP1[4]*mP1[4] + mP1[5]*mP1[5] );
279 n*=n1;
280 Double_t a = 0;
281 if( n>1.e-8 ) a = ( mP[3]*mP1[3] + mP[4]*mP1[4] + mP[5]*mP1[5] )/n;
282 if (TMath::Abs(a)<1.) a = TMath::ACos(a);
283 else a = (a>=0) ?0 :TMath::Pi();
284 return a;
285}
286
287Double_t AliKFParticle::GetAngleXY( const AliKFParticle &p ) const
288{
289 //* Calculate the opening angle between two particles in XY plane
290
291 Double_t dS, dS1;
292 GetDStoParticleXY( p, dS, dS1 );
293 Double_t mP[8], mC[36], mP1[8], mC1[36];
294 Transport( dS, mP, mC );
295 p.Transport( dS1, mP1, mC1 );
296 Double_t n = TMath::Sqrt( mP[3]*mP[3] + mP[4]*mP[4] );
297 Double_t n1= TMath::Sqrt( mP1[3]*mP1[3] + mP1[4]*mP1[4] );
298 n*=n1;
299 Double_t a = 0;
300 if( n>1.e-8 ) a = ( mP[3]*mP1[3] + mP[4]*mP1[4] )/n;
301 if (TMath::Abs(a)<1.) a = TMath::ACos(a);
302 else a = (a>=0) ?0 :TMath::Pi();
303 return a;
304}
305
306Double_t AliKFParticle::GetAngleRZ( const AliKFParticle &p ) const
307{
308 //* Calculate the opening angle between two particles in RZ plane
616ffc76 309
e7b09c95 310 Double_t dS, dS1;
311 GetDStoParticle( p, dS, dS1 );
312 Double_t mP[8], mC[36], mP1[8], mC1[36];
313 Transport( dS, mP, mC );
314 p.Transport( dS1, mP1, mC1 );
315 Double_t nr = TMath::Sqrt( mP[3]*mP[3] + mP[4]*mP[4] );
316 Double_t n1r= TMath::Sqrt( mP1[3]*mP1[3] + mP1[4]*mP1[4] );
317 Double_t n = TMath::Sqrt( nr*nr + mP[5]*mP[5] );
318 Double_t n1= TMath::Sqrt( n1r*n1r + mP1[5]*mP1[5] );
319 n*=n1;
320 Double_t a = 0;
321 if( n>1.e-8 ) a = ( nr*n1r +mP[5]*mP1[5])/n;
322 if (TMath::Abs(a)<1.) a = TMath::ACos(a);
323 else a = (a>=0) ?0 :TMath::Pi();
324 return a;
616ffc76 325}
706952f5 326
327
328/*
329
330#include "AliExternalTrackParam.h"
331
332void AliKFParticle::GetDStoParticleALICE( const AliKFParticleBase &p,
333 Double_t &DS, Double_t &DS1 )
334 const
335{
336 DS = DS1 = 0;
337 Double_t x1, a1, x2, a2;
338 Double_t par1[5], par2[5], cov[15];
339 for(int i=0; i<15; i++) cov[i] = 0;
340 cov[0] = cov[2] = cov[5] = cov[9] = cov[14] = .001;
341
342 GetExternalTrackParam( *this, x1, a1, par1 );
343 GetExternalTrackParam( p, x2, a2, par2 );
344
345 AliExternalTrackParam t1(x1,a1, par1, cov);
346 AliExternalTrackParam t2(x2,a2, par2, cov);
347
348 Double_t xe1=0, xe2=0;
349 t1.GetDCA( &t2, -GetFieldAlice(), xe1, xe2 );
350 t1.PropagateTo( xe1, -GetFieldAlice() );
351 t2.PropagateTo( xe2, -GetFieldAlice() );
352
353 Double_t xyz1[3], xyz2[3];
354 t1.GetXYZ( xyz1 );
355 t2.GetXYZ( xyz2 );
356
357 DS = GetDStoPoint( xyz1 );
358 DS1 = p.GetDStoPoint( xyz2 );
359
360 return;
361}
362*/