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