In AddTaskPHOSPi0Flow.C set Cent. Bin past event buffers/lists,
[u/mrichter/AliRoot.git] / STEER / ESD / 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 AliKFParticle::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
40 void AliKFParticle::Create( const Double_t Param[], const Double_t Cov[], Int_t Charge, Int_t PID )
41 {
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 )
53   Double_t C[21];
54   for( int i=0; i<21; i++ ) C[i] = Cov[i];
55   
56   TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(PID);
57   Double_t mass = (particlePDG) ? particlePDG->Mass() :0.13957;
58   
59   AliKFParticleBase::Initialize( Param, C, Charge, mass );
60 }
61
62 AliKFParticle::AliKFParticle( const AliVTrack &track, Int_t PID )
63 {
64   // Constructor from ALICE track, PID hypothesis should be provided
65
66   track.GetXYZ(fP);
67   track.PxPyPz(fP+3);
68   fQ = track.Charge();
69   track.GetCovarianceXYZPxPyPz( fC );
70   Create(fP,fC,fQ,PID);
71 }
72
73 AliKFParticle::AliKFParticle( const AliVVertex &vertex )
74 {
75   // Constructor from ALICE vertex
76
77   vertex.GetXYZ( fP );
78   vertex.GetCovarianceMatrix( fC );  
79   fChi2 = vertex.GetChi2();
80   fNDF = 2*vertex.GetNContributors() - 3;
81   fQ = 0;
82   fAtProductionVertex = 0;
83   fIsLinearized = 0;
84   fSFromDecay = 0;
85 }
86
87 void AliKFParticle::GetExternalTrackParam( const AliKFParticleBase &p, Double_t &X, Double_t &Alpha, Double_t P[5] ) 
88 {
89   // Conversion to AliExternalTrackParam parameterization
90
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
111 Bool_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;
135     ey = py/pt;
136     val = dy*ex - dx*ey;
137   }
138
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) );
149
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
159 Bool_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
165 Bool_t AliKFParticle::GetDistanceFromVertexXY( const AliKFParticle &Vtx, Double_t &val, Double_t &err ) const 
166 {
167   //* Calculate distance from vertex [cm] in XY-plane
168
169   return GetDistanceFromVertexXY( Vtx.fP, Vtx.fC, val, err );
170 }
171
172 Bool_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
179 Double_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;
185 }
186
187 Double_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
194 Double_t AliKFParticle::GetDistanceFromVertexXY( const AliVVertex &Vtx ) const 
195 {
196   //* Calculate distance from vertex [cm] in XY-plane
197
198   return GetDistanceFromVertexXY( AliKFParticle(Vtx).fP );
199 }
200
201 Double_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
215 Double_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 ); 
223
224   Double_t d[2]={ fP[0]-mP1[0], fP[1]-mP1[1] };
225
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] )  );
228
229   Double_t h[2] = { mP1[3]*sigmaS, mP1[4]*sigmaS };       
230   
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
239 Double_t AliKFParticle::GetDeviationFromVertexXY( const Double_t vtx[], const Double_t Cv[] ) const 
240 {
241   //* Calculate sqrt(Chi2/ndf) deviation from vertex
242   //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
243
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;
248 }
249
250
251 Double_t AliKFParticle::GetDeviationFromVertexXY( const AliKFParticle &Vtx ) const  
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
259 Double_t AliKFParticle::GetDeviationFromVertexXY( const AliVVertex &Vtx ) const 
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
268 Double_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
287 Double_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
306 Double_t AliKFParticle::GetAngleRZ( const AliKFParticle &p ) const 
307 {
308   //* Calculate the opening angle between two particles in RZ plane
309
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;
325 }
326
327
328 /*
329
330 #include "AliExternalTrackParam.h"
331
332 void 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 */
363
364   // * Pseudo Proper Time of decay = (r*pt) / |pt| * M/|pt|
365 Double_t AliKFParticle::GetPseudoProperDecayTime( const AliKFParticle &pV, const Double_t& mass, Double_t* timeErr2 ) const
366 { // TODO optimize with respect to time and stability
367   const Double_t ipt2 = 1/( Px()*Px() + Py()*Py() );
368   const Double_t mipt2 = mass*ipt2;
369   const Double_t dx = X() - pV.X();
370   const Double_t dy = Y() - pV.Y();
371
372   if ( timeErr2 ) {
373       // -- calculate error = sigma(f(r)) = f'Cf'
374       // r = {x,y,px,py,x_pV,y_pV}
375       // df/dr = { px*m/pt^2,
376       //     py*m/pt^2,
377       //    ( x - x_pV )*m*(1/pt^2 - 2(px/pt^2)^2),
378       //    ( y - y_pV )*m*(1/pt^2 - 2(py/pt^2)^2),
379       //     -px*m/pt^2,
380       //     -py*m/pt^2 }
381     const Double_t f0 = Px()*mipt2;
382     const Double_t f1 = Py()*mipt2;
383     const Double_t mipt2derivative = mipt2*(1-2*Px()*Px()*ipt2);
384     const Double_t f2 = dx*mipt2derivative;
385     const Double_t f3 = -dy*mipt2derivative;
386     const Double_t f4 = -f0;
387     const Double_t f5 = -f1;
388
389     const Double_t& mC00 =    GetCovariance(0,0);
390     const Double_t& mC10 =    GetCovariance(0,1);
391     const Double_t& mC11 =    GetCovariance(1,1);
392     const Double_t& mC20 =    GetCovariance(3,0);
393     const Double_t& mC21 =    GetCovariance(3,1);
394     const Double_t& mC22 =    GetCovariance(3,3);
395     const Double_t& mC30 =    GetCovariance(4,0);
396     const Double_t& mC31 =    GetCovariance(4,1);
397     const Double_t& mC32 =    GetCovariance(4,3);
398     const Double_t& mC33 =    GetCovariance(4,4);
399     const Double_t& mC44 = pV.GetCovariance(0,0);
400     const Double_t& mC54 = pV.GetCovariance(1,0);
401     const Double_t& mC55 = pV.GetCovariance(1,1);
402
403     *timeErr2 =
404       f5*mC55*f5 +
405       f5*mC54*f4 +
406       f4*mC44*f4 +
407       f3*mC33*f3 +
408       f3*mC32*f2 +
409       f3*mC31*f1 +
410       f3*mC30*f0 +
411       f2*mC22*f2 +
412       f2*mC21*f1 +
413       f2*mC20*f0 +
414       f1*mC11*f1 +
415       f1*mC10*f0 +
416       f0*mC00*f0;
417   }
418   return ( dx*Px() + dy*Py() )*mipt2;
419 }