Fixes for bug #52499: Field polarities inconsistiency
[u/mrichter/AliRoot.git] / STEER / AliKFParticle.cxx
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"
22 #include "AliVTrack.h"
23 #include "AliVVertex.h"
24
25 ClassImp(AliKFParticle)
26
27 Double_t AliKFParticle::fgBz = -5.;  //* Bz compoment of the magnetic field
28
29 void AliKFParticle::Create( const Double_t Param[], const Double_t Cov[], Int_t Charge, Int_t PID )
30 {
31   // Constructor from "cartesian" track, PID hypothesis should be provided
32   //
33   // Param[6] = { X, Y, Z, Px, Py, Pz } - position and momentum
34   // Cov [21] = lower-triangular part of the covariance matrix:
35   //
36   //                (  0  .  .  .  .  . )
37   //                (  1  2  .  .  .  . )
38   //  Cov. matrix = (  3  4  5  .  .  . ) - numbering of covariance elements in Cov[]
39   //                (  6  7  8  9  .  . )
40   //                ( 10 11 12 13 14  . )
41   //                ( 15 16 17 18 19 20 )
42   Double_t C[21];
43   for( int i=0; i<21; i++ ) C[i] = Cov[i];
44   
45   TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(PID);
46   Double_t mass = (particlePDG) ? particlePDG->Mass() :0.13957;
47   
48   AliKFParticleBase::Initialize( Param, C, Charge, mass );
49 }
50
51
52 AliKFParticle::AliKFParticle( const AliVTrack &track, Int_t PID )
53 {
54   // Constructor from ALICE track, PID hypothesis should be provided
55
56   track.XvYvZv(fP);
57   track.PxPyPz(fP+3);
58   fQ = track.Charge();
59   track.GetCovarianceXYZPxPyPz( fC );
60   Create(fP,fC,fQ,PID);
61 }
62
63 AliKFParticle::AliKFParticle( const AliVVertex &vertex )
64 {
65   // Constructor from ALICE vertex
66
67   vertex.GetXYZ( fP );
68   vertex.GetCovarianceMatrix( fC );  
69   fChi2 = vertex.GetChi2();
70   fNDF = 2*vertex.GetNContributors() - 3;
71   fQ = 0;
72   fAtProductionVertex = 0;
73   fIsLinearized = 0;
74   fSFromDecay = 0;
75 }
76
77
78 void AliKFParticle::GetExternalTrackParam( const AliKFParticleBase &p, Double_t &X, Double_t &Alpha, Double_t P[5] ) 
79 {
80   // Conversion to AliExternalTrackParam parameterization
81
82   Double_t cosA = p.GetPx(), sinA = p.GetPy(); 
83   Double_t pt = TMath::Sqrt(cosA*cosA + sinA*sinA);
84   Double_t pti = 0;
85   if( pt<1.e-4 ){
86     cosA = 1;
87     sinA = 0;
88   } else {
89     pti = 1./pt;
90     cosA*=pti;
91     sinA*=pti;
92   }
93   Alpha = TMath::ATan2(sinA,cosA);  
94   X   = p.GetX()*cosA + p.GetY()*sinA;   
95   P[0]= p.GetY()*cosA - p.GetX()*sinA;
96   P[1]= p.GetZ();
97   P[2]= 0;
98   P[3]= p.GetPz()*pti;
99   P[4]= p.GetQ()*pti;
100 }
101
102
103
104 Double_t AliKFParticle::GetDistanceFromVertexXY( const Double_t vtx[] ) const
105 {
106   //* Calculate distance from vertex [cm] in XY-plane
107
108   Double_t mP[8], mC[36];
109   Transport( GetDStoPoint(vtx), mP, mC );
110   Double_t d[2]={ vtx[0]-mP[0], vtx[1]-mP[1] };
111   Double_t dist =  TMath::Sqrt( d[0]*d[0]+d[1]*d[1] );
112   Double_t sign = d[0]*mP[3] - d[1]*mP[4];  
113   return (sign>=0) ?dist :-dist;
114 }
115
116 Double_t AliKFParticle::GetDistanceFromVertexXY( const AliKFParticle &Vtx ) const 
117 {
118   //* Calculate distance from vertex [cm] in XY-plane
119
120   return GetDistanceFromVertexXY( Vtx.fP );
121 }
122
123 Double_t AliKFParticle::GetDistanceFromVertexXY( const AliVVertex &Vtx ) const 
124 {
125   //* Calculate distance from vertex [cm] in XY-plane
126
127   return GetDistanceFromVertexXY( AliKFParticle(Vtx).fP );
128 }
129
130 Double_t AliKFParticle::GetDistanceFromParticleXY( const AliKFParticle &p ) const 
131 {
132   //* Calculate distance to other particle [cm]
133
134   Double_t dS, dS1;
135   GetDStoParticleXY( p, dS, dS1 );   
136   Double_t mP[8], mC[36], mP1[8], mC1[36];
137   Transport( dS, mP, mC ); 
138   p.Transport( dS1, mP1, mC1 ); 
139   Double_t dx = mP[0]-mP1[0]; 
140   Double_t dy = mP[1]-mP1[1]; 
141   return TMath::Sqrt(dx*dx+dy*dy);
142 }
143
144 Double_t AliKFParticle::GetDeviationFromParticleXY( const AliKFParticle &p ) const 
145 {
146   //* Calculate sqrt(Chi2/ndf) deviation from other particle
147
148   Double_t dS, dS1;
149   GetDStoParticleXY( p, dS, dS1 );   
150   Double_t mP1[8], mC1[36];
151   p.Transport( dS1, mP1, mC1 ); 
152
153   Double_t d[2]={ fP[0]-mP1[0], fP[1]-mP1[1] };
154
155   Double_t sigmaS = .1+10.*TMath::Sqrt( (d[0]*d[0]+d[1]*d[1] )/
156                                         (mP1[3]*mP1[3]+mP1[4]*mP1[4] )  );
157
158   Double_t h[2] = { mP1[3]*sigmaS, mP1[4]*sigmaS };       
159   
160   mC1[0] +=h[0]*h[0];
161   mC1[1] +=h[1]*h[0]; 
162   mC1[2] +=h[1]*h[1]; 
163
164   return GetDeviationFromVertexXY( mP1, mC1 )*TMath::Sqrt(2./1.);
165 }
166
167
168 Double_t AliKFParticle::GetDeviationFromVertexXY( const Double_t v[], const Double_t Cv[] ) const 
169 {
170   //* Calculate sqrt(Chi2/ndf) deviation from vertex
171   //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
172
173   Double_t mP[8];
174   Double_t mC[36];
175   
176   Transport( GetDStoPoint(v), mP, mC );  
177
178   Double_t d[2]={ v[0]-mP[0], v[1]-mP[1] };
179
180   Double_t sigmaS = .1+10.*TMath::Sqrt( (d[0]*d[0]+d[1]*d[1] )/
181                                         (mP[3]*mP[3]+mP[4]*mP[4] )  );
182    
183   Double_t h[2] = { mP[3]*sigmaS, mP[4]*sigmaS };       
184   
185   Double_t mSi[3] = { mC[0] +h[0]*h[0], 
186                       mC[1] +h[1]*h[0], mC[2] +h[1]*h[1] };
187
188   if( Cv ){
189     mSi[0]+=Cv[0];
190     mSi[1]+=Cv[1];
191     mSi[2]+=Cv[2];
192   }
193   Double_t s = ( mSi[0]*mSi[2] - mSi[1]*mSi[1] );
194   s = ( s > 1.E-20 )  ?1./s :0;   
195   
196   Double_t mS[3] = { mSi[2], 
197                      -mSi[1], mSi[0] };      
198
199   return TMath::Sqrt( TMath::Abs(s*( ( mS[0]*d[0] + mS[1]*d[1] )*d[0]
200                                      +(mS[1]*d[0] + mS[2]*d[1] )*d[1] ))/1);
201 }
202
203
204 Double_t AliKFParticle::GetDeviationFromVertexXY( const AliKFParticle &Vtx ) const  
205 {
206   //* Calculate sqrt(Chi2/ndf) deviation from vertex
207   //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
208
209   return GetDeviationFromVertexXY( Vtx.fP, Vtx.fC );
210 }
211
212 Double_t AliKFParticle::GetDeviationFromVertexXY( const AliVVertex &Vtx ) const 
213 {
214   //* Calculate sqrt(Chi2/ndf) deviation from vertex
215   //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
216
217   AliKFParticle v(Vtx);
218   return GetDeviationFromVertexXY( v.fP, v.fC );
219 }
220
221
222 Double_t AliKFParticle::GetAngle  ( const AliKFParticle &p ) const 
223 {
224   //* Calculate the opening angle between two particles
225
226   Double_t dS, dS1;
227   GetDStoParticle( p, dS, dS1 );   
228   Double_t mP[8], mC[36], mP1[8], mC1[36];
229   Transport( dS, mP, mC ); 
230   p.Transport( dS1, mP1, mC1 ); 
231   Double_t n = TMath::Sqrt( mP[3]*mP[3] + mP[4]*mP[4] + mP[5]*mP[5] );
232   Double_t n1= TMath::Sqrt( mP1[3]*mP1[3] + mP1[4]*mP1[4] + mP1[5]*mP1[5] );
233   n*=n1;
234   Double_t a = 0;
235   if( n>1.e-8 ) a = ( mP[3]*mP1[3] + mP[4]*mP1[4] + mP[5]*mP1[5] )/n;
236   if (TMath::Abs(a)<1.) a = TMath::ACos(a);
237   else a = (a>=0) ?0 :TMath::Pi();
238   return a;
239 }
240
241 Double_t AliKFParticle::GetAngleXY( const AliKFParticle &p ) const 
242 {
243   //* Calculate the opening angle between two particles in XY plane
244
245   Double_t dS, dS1;
246   GetDStoParticleXY( p, dS, dS1 );   
247   Double_t mP[8], mC[36], mP1[8], mC1[36];
248   Transport( dS, mP, mC ); 
249   p.Transport( dS1, mP1, mC1 ); 
250   Double_t n = TMath::Sqrt( mP[3]*mP[3] + mP[4]*mP[4] );
251   Double_t n1= TMath::Sqrt( mP1[3]*mP1[3] + mP1[4]*mP1[4] );
252   n*=n1;
253   Double_t a = 0;
254   if( n>1.e-8 ) a = ( mP[3]*mP1[3] + mP[4]*mP1[4] )/n;
255   if (TMath::Abs(a)<1.) a = TMath::ACos(a);
256   else a = (a>=0) ?0 :TMath::Pi();
257   return a;
258 }
259
260 Double_t AliKFParticle::GetAngleRZ( const AliKFParticle &p ) const 
261 {
262   //* Calculate the opening angle between two particles in RZ plane
263
264   Double_t dS, dS1;
265   GetDStoParticle( p, dS, dS1 );   
266   Double_t mP[8], mC[36], mP1[8], mC1[36];
267   Transport( dS, mP, mC ); 
268   p.Transport( dS1, mP1, mC1 ); 
269   Double_t nr = TMath::Sqrt( mP[3]*mP[3] + mP[4]*mP[4] );
270   Double_t n1r= TMath::Sqrt( mP1[3]*mP1[3] + mP1[4]*mP1[4]  );
271   Double_t n = TMath::Sqrt( nr*nr + mP[5]*mP[5] );
272   Double_t n1= TMath::Sqrt( n1r*n1r + mP1[5]*mP1[5] );
273   n*=n1;
274   Double_t a = 0;
275   if( n>1.e-8 ) a = ( nr*n1r +mP[5]*mP1[5])/n; 
276   if (TMath::Abs(a)<1.) a = TMath::ACos(a);
277   else a = (a>=0) ?0 :TMath::Pi();
278   return a;
279 }
280
281
282 /*
283
284 #include "AliExternalTrackParam.h"
285
286 void AliKFParticle::GetDStoParticleALICE( const AliKFParticleBase &p, 
287                                            Double_t &DS, Double_t &DS1 ) 
288   const
289
290   DS = DS1 = 0;   
291   Double_t x1, a1, x2, a2;
292   Double_t par1[5], par2[5], cov[15];
293   for(int i=0; i<15; i++) cov[i] = 0;
294   cov[0] = cov[2] = cov[5] = cov[9] = cov[14] = .001;
295
296   GetExternalTrackParam( *this, x1, a1, par1 );
297   GetExternalTrackParam( p, x2, a2, par2 );
298
299   AliExternalTrackParam t1(x1,a1, par1, cov);
300   AliExternalTrackParam t2(x2,a2, par2, cov);
301
302   Double_t xe1=0, xe2=0;
303   t1.GetDCA( &t2, -GetFieldAlice(), xe1, xe2 );
304   t1.PropagateTo( xe1, -GetFieldAlice() );
305   t2.PropagateTo( xe2, -GetFieldAlice() );
306
307   Double_t xyz1[3], xyz2[3];
308   t1.GetXYZ( xyz1 );
309   t2.GetXYZ( xyz2 );
310   
311   DS = GetDStoPoint( xyz1 );
312   DS1 = p.GetDStoPoint( xyz2 );
313
314   return;
315 }
316 */