Corrected protection.
[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 AliKFParticle::AliKFParticle( const AliVTrack &track, Int_t PID )
52 {
53   // Constructor from ALICE track, PID hypothesis should be provided
54
55   track.GetXYZ(fP);
56   track.PxPyPz(fP+3);
57   fQ = track.Charge();
58   track.GetCovarianceXYZPxPyPz( fC );
59   Create(fP,fC,fQ,PID);
60 }
61
62 AliKFParticle::AliKFParticle( const AliVVertex &vertex )
63 {
64   // Constructor from ALICE vertex
65
66   vertex.GetXYZ( fP );
67   vertex.GetCovarianceMatrix( fC );  
68   fChi2 = vertex.GetChi2();
69   fNDF = 2*vertex.GetNContributors() - 3;
70   fQ = 0;
71   fAtProductionVertex = 0;
72   fIsLinearized = 0;
73   fSFromDecay = 0;
74 }
75
76 void AliKFParticle::GetExternalTrackParam( const AliKFParticleBase &p, Double_t &X, Double_t &Alpha, Double_t P[5] ) 
77 {
78   // Conversion to AliExternalTrackParam parameterization
79
80   Double_t cosA = p.GetPx(), sinA = p.GetPy(); 
81   Double_t pt = TMath::Sqrt(cosA*cosA + sinA*sinA);
82   Double_t pti = 0;
83   if( pt<1.e-4 ){
84     cosA = 1;
85     sinA = 0;
86   } else {
87     pti = 1./pt;
88     cosA*=pti;
89     sinA*=pti;
90   }
91   Alpha = TMath::ATan2(sinA,cosA);  
92   X   = p.GetX()*cosA + p.GetY()*sinA;   
93   P[0]= p.GetY()*cosA - p.GetX()*sinA;
94   P[1]= p.GetZ();
95   P[2]= 0;
96   P[3]= p.GetPz()*pti;
97   P[4]= p.GetQ()*pti;
98 }
99
100 Bool_t AliKFParticle::GetDistanceFromVertexXY( const Double_t vtx[], const Double_t Cv[], Double_t &val, Double_t &err ) const
101 {
102   //* Calculate DCA distance from vertex (transverse impact parameter) in XY
103   //* v = [xy], Cv=[Cxx,Cxy,Cyy ]-covariance matrix
104
105   Bool_t ret = 0;
106   
107   Double_t mP[8];
108   Double_t mC[36];
109   
110   Transport( GetDStoPoint(vtx), mP, mC );  
111
112   Double_t dx = mP[0] - vtx[0];
113   Double_t dy = mP[1] - vtx[1];
114   Double_t px = mP[3];
115   Double_t py = mP[4];
116   Double_t pt = TMath::Sqrt(px*px + py*py);
117   Double_t ex=0, ey=0;
118   if( pt<1.e-4 ){
119     ret = 1;
120     pt = 1.;
121     val = 1.e4;
122   } else{
123     ex = px/pt;
124     ey = py/pt;
125     val = dy*ex - dx*ey;
126   }
127
128   Double_t h0 = -ey;
129   Double_t h1 = ex;
130   Double_t h3 = (dy*ey + dx*ex)*ey/pt;
131   Double_t h4 = -(dy*ey + dx*ex)*ex/pt;
132   
133   err = 
134     h0*(h0*GetCovariance(0,0) + h1*GetCovariance(0,1) + h3*GetCovariance(0,3) + h4*GetCovariance(0,4) ) +
135     h1*(h0*GetCovariance(1,0) + h1*GetCovariance(1,1) + h3*GetCovariance(1,3) + h4*GetCovariance(1,4) ) +
136     h3*(h0*GetCovariance(3,0) + h1*GetCovariance(3,1) + h3*GetCovariance(3,3) + h4*GetCovariance(3,4) ) +
137     h4*(h0*GetCovariance(4,0) + h1*GetCovariance(4,1) + h3*GetCovariance(4,3) + h4*GetCovariance(4,4) );
138
139   if( Cv ){
140     err+= h0*(h0*Cv[0] + h1*Cv[1] ) + h1*(h0*Cv[1] + h1*Cv[2] ); 
141   }
142
143   err = TMath::Sqrt(TMath::Abs(err));
144
145   return ret;
146 }
147
148 Bool_t AliKFParticle::GetDistanceFromVertexXY( const Double_t vtx[], Double_t &val, Double_t &err ) const
149 {
150   return GetDistanceFromVertexXY( vtx, 0, val, err );
151 }
152
153
154 Bool_t AliKFParticle::GetDistanceFromVertexXY( const AliKFParticle &Vtx, Double_t &val, Double_t &err ) const 
155 {
156   //* Calculate distance from vertex [cm] in XY-plane
157
158   return GetDistanceFromVertexXY( Vtx.fP, Vtx.fC, val, err );
159 }
160
161 Bool_t AliKFParticle::GetDistanceFromVertexXY( const AliVVertex &Vtx, Double_t &val, Double_t &err ) const 
162 {
163   //* Calculate distance from vertex [cm] in XY-plane
164
165   return GetDistanceFromVertexXY( AliKFParticle(Vtx), val, err );
166 }
167
168 Double_t AliKFParticle::GetDistanceFromVertexXY( const Double_t vtx[] ) const
169 {
170   //* Calculate distance from vertex [cm] in XY-plane
171   Double_t val, err;
172   GetDistanceFromVertexXY( vtx, 0, val, err );
173   return val;
174 }
175
176 Double_t AliKFParticle::GetDistanceFromVertexXY( const AliKFParticle &Vtx ) const 
177 {
178   //* Calculate distance from vertex [cm] in XY-plane
179
180   return GetDistanceFromVertexXY( Vtx.fP );
181 }
182
183 Double_t AliKFParticle::GetDistanceFromVertexXY( const AliVVertex &Vtx ) const 
184 {
185   //* Calculate distance from vertex [cm] in XY-plane
186
187   return GetDistanceFromVertexXY( AliKFParticle(Vtx).fP );
188 }
189
190 Double_t AliKFParticle::GetDistanceFromParticleXY( const AliKFParticle &p ) const 
191 {
192   //* Calculate distance to other particle [cm]
193
194   Double_t dS, dS1;
195   GetDStoParticleXY( p, dS, dS1 );   
196   Double_t mP[8], mC[36], mP1[8], mC1[36];
197   Transport( dS, mP, mC ); 
198   p.Transport( dS1, mP1, mC1 ); 
199   Double_t dx = mP[0]-mP1[0]; 
200   Double_t dy = mP[1]-mP1[1]; 
201   return TMath::Sqrt(dx*dx+dy*dy);
202 }
203
204 Double_t AliKFParticle::GetDeviationFromParticleXY( const AliKFParticle &p ) const 
205 {
206   //* Calculate sqrt(Chi2/ndf) deviation from other particle
207
208   Double_t dS, dS1;
209   GetDStoParticleXY( p, dS, dS1 );   
210   Double_t mP1[8], mC1[36];
211   p.Transport( dS1, mP1, mC1 ); 
212
213   Double_t d[2]={ fP[0]-mP1[0], fP[1]-mP1[1] };
214
215   Double_t sigmaS = .1+10.*TMath::Sqrt( (d[0]*d[0]+d[1]*d[1] )/
216                                         (mP1[3]*mP1[3]+mP1[4]*mP1[4] )  );
217
218   Double_t h[2] = { mP1[3]*sigmaS, mP1[4]*sigmaS };       
219   
220   mC1[0] +=h[0]*h[0];
221   mC1[1] +=h[1]*h[0]; 
222   mC1[2] +=h[1]*h[1]; 
223
224   return GetDeviationFromVertexXY( mP1, mC1 )*TMath::Sqrt(2./1.);
225 }
226
227
228 Double_t AliKFParticle::GetDeviationFromVertexXY( const Double_t vtx[], const Double_t Cv[] ) const 
229 {
230   //* Calculate sqrt(Chi2/ndf) deviation from vertex
231   //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
232
233   Double_t val, err;
234   Bool_t problem = GetDistanceFromVertexXY( vtx, Cv, val, err );
235   if( problem || err<1.e-20 ) return 1.e4;
236   else return val/err;
237 }
238
239
240 Double_t AliKFParticle::GetDeviationFromVertexXY( const AliKFParticle &Vtx ) const  
241 {
242   //* Calculate sqrt(Chi2/ndf) deviation from vertex
243   //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
244
245   return GetDeviationFromVertexXY( Vtx.fP, Vtx.fC );
246 }
247
248 Double_t AliKFParticle::GetDeviationFromVertexXY( const AliVVertex &Vtx ) const 
249 {
250   //* Calculate sqrt(Chi2/ndf) deviation from vertex
251   //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
252
253   AliKFParticle v(Vtx);
254   return GetDeviationFromVertexXY( v.fP, v.fC );
255 }
256
257 Double_t AliKFParticle::GetAngle  ( const AliKFParticle &p ) const 
258 {
259   //* Calculate the opening angle between two particles
260
261   Double_t dS, dS1;
262   GetDStoParticle( p, dS, dS1 );   
263   Double_t mP[8], mC[36], mP1[8], mC1[36];
264   Transport( dS, mP, mC ); 
265   p.Transport( dS1, mP1, mC1 ); 
266   Double_t n = TMath::Sqrt( mP[3]*mP[3] + mP[4]*mP[4] + mP[5]*mP[5] );
267   Double_t n1= TMath::Sqrt( mP1[3]*mP1[3] + mP1[4]*mP1[4] + mP1[5]*mP1[5] );
268   n*=n1;
269   Double_t a = 0;
270   if( n>1.e-8 ) a = ( mP[3]*mP1[3] + mP[4]*mP1[4] + mP[5]*mP1[5] )/n;
271   if (TMath::Abs(a)<1.) a = TMath::ACos(a);
272   else a = (a>=0) ?0 :TMath::Pi();
273   return a;
274 }
275
276 Double_t AliKFParticle::GetAngleXY( const AliKFParticle &p ) const 
277 {
278   //* Calculate the opening angle between two particles in XY plane
279
280   Double_t dS, dS1;
281   GetDStoParticleXY( p, dS, dS1 );   
282   Double_t mP[8], mC[36], mP1[8], mC1[36];
283   Transport( dS, mP, mC ); 
284   p.Transport( dS1, mP1, mC1 ); 
285   Double_t n = TMath::Sqrt( mP[3]*mP[3] + mP[4]*mP[4] );
286   Double_t n1= TMath::Sqrt( mP1[3]*mP1[3] + mP1[4]*mP1[4] );
287   n*=n1;
288   Double_t a = 0;
289   if( n>1.e-8 ) a = ( mP[3]*mP1[3] + mP[4]*mP1[4] )/n;
290   if (TMath::Abs(a)<1.) a = TMath::ACos(a);
291   else a = (a>=0) ?0 :TMath::Pi();
292   return a;
293 }
294
295 Double_t AliKFParticle::GetAngleRZ( const AliKFParticle &p ) const 
296 {
297   //* Calculate the opening angle between two particles in RZ plane
298
299   Double_t dS, dS1;
300   GetDStoParticle( p, dS, dS1 );   
301   Double_t mP[8], mC[36], mP1[8], mC1[36];
302   Transport( dS, mP, mC ); 
303   p.Transport( dS1, mP1, mC1 ); 
304   Double_t nr = TMath::Sqrt( mP[3]*mP[3] + mP[4]*mP[4] );
305   Double_t n1r= TMath::Sqrt( mP1[3]*mP1[3] + mP1[4]*mP1[4]  );
306   Double_t n = TMath::Sqrt( nr*nr + mP[5]*mP[5] );
307   Double_t n1= TMath::Sqrt( n1r*n1r + mP1[5]*mP1[5] );
308   n*=n1;
309   Double_t a = 0;
310   if( n>1.e-8 ) a = ( nr*n1r +mP[5]*mP1[5])/n; 
311   if (TMath::Abs(a)<1.) a = TMath::ACos(a);
312   else a = (a>=0) ?0 :TMath::Pi();
313   return a;
314 }
315
316
317 /*
318
319 #include "AliExternalTrackParam.h"
320
321 void AliKFParticle::GetDStoParticleALICE( const AliKFParticleBase &p, 
322                                            Double_t &DS, Double_t &DS1 ) 
323   const
324
325   DS = DS1 = 0;   
326   Double_t x1, a1, x2, a2;
327   Double_t par1[5], par2[5], cov[15];
328   for(int i=0; i<15; i++) cov[i] = 0;
329   cov[0] = cov[2] = cov[5] = cov[9] = cov[14] = .001;
330
331   GetExternalTrackParam( *this, x1, a1, par1 );
332   GetExternalTrackParam( p, x2, a2, par2 );
333
334   AliExternalTrackParam t1(x1,a1, par1, cov);
335   AliExternalTrackParam t2(x2,a2, par2, cov);
336
337   Double_t xe1=0, xe2=0;
338   t1.GetDCA( &t2, -GetFieldAlice(), xe1, xe2 );
339   t1.PropagateTo( xe1, -GetFieldAlice() );
340   t2.PropagateTo( xe2, -GetFieldAlice() );
341
342   Double_t xyz1[3], xyz2[3];
343   t1.GetXYZ( xyz1 );
344   t2.GetXYZ( xyz2 );
345   
346   DS = GetDStoPoint( xyz1 );
347   DS1 = p.GetDStoPoint( xyz2 );
348
349   return;
350 }
351 */