Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
KFParticle.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file KFParticle.cxx
1 /*
2  * This file is part of KFParticle package
3  * Copyright (C) 2007-2019 FIAS Frankfurt Institute for Advanced Studies
4  * 2007-2019 Goethe University of Frankfurt
5  * 2007-2019 Ivan Kisel <I.Kisel@compeng.uni-frankfurt.de>
6  * 2007-2019 Maksym Zyzak
7  * 2007-2019 Sergey Gorbunov
8  *
9  * KFParticle is free software: you can redistribute it and/or modify
10  * it under the terms of the GNU General Public License as published by
11  * the Free Software Foundation, either version 3 of the License, or
12  * (at your option) any later version.
13  *
14  * KFParticle is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  * GNU General Public License for more details.
18  *
19  * You should have received a copy of the GNU General Public License
20  * along with this program. If not, see <https://www.gnu.org/licenses/>.
21  */
22 
23 
24 #include "KFParticle.h"
25 #include "KFParticleDatabase.h"
26 
27 #include "KFPTrack.h"
28 #include "KFPVertex.h"
29 
30 #ifndef KFParticleStandalone
31 ClassImp(KFParticle);
32 #endif
33 
34 
35 #ifdef HomogeneousField
36 float KFParticle::fgBz = -5.; //* Bz compoment of the magnetic field
37 #endif
38 
39 KFParticle::KFParticle( const KFParticle &d1, const KFParticle &d2 ): KFParticleBase()
40 {
45  KFParticle mother;
46  mother+= d1;
47  mother+= d2;
48  *this = mother;
49 }
50 
51 void KFParticle::Create( const float Param[], const float Cov[], Int_t Charge, float mass )
52 {
67  float C[21];
68  for( int i=0; i<21; i++ ) C[i] = Cov[i];
69 
70  KFParticleBase::Initialize( Param, C, Charge, mass );
71 }
72 
73 void KFParticle::Create( const Double_t Param[], const Double_t Cov[], Int_t Charge, float mass )
74 {
89  float P[6];
90  for(int i=0; i<6; i++ ) P[i] = Param[i];
91  float C[21];
92  for( int i=0; i<21; i++ ) C[i] = Cov[i];
93 
94  KFParticleBase::Initialize( P, C, Charge, mass );
95 }
96 
97 KFParticle::KFParticle( const KFPTrack &track, const int PID ): KFParticleBase()
98 {
104  track.XvYvZv(fP);
105  track.PxPyPz(fP+3);
106  fQ = track.Charge();
107  track.GetCovarianceXYZPxPyPz( fC );
108 
109  float mass = KFParticleDatabase::Instance()->GetMass(PID);
110 
111  Create(fP,fC,fQ,mass);
112  fChi2 = track.GetChi2();
113  fNDF = track.GetNDF();
114  SetPDG(PID);
115 #ifdef NonhomogeneousField
116  for(int iF=0; iF<10; iF++)
117  SetFieldCoeff( track.GetFieldCoeff()[iF], iF);
118 #endif
119 }
120 
121 KFParticle::KFParticle( const KFPVertex &vertex ): KFParticleBase()
122 {
127  vertex.GetXYZ( fP );
128  vertex.GetCovarianceMatrix( fC );
129  fChi2 = vertex.GetChi2();
130  fNDF = 2*vertex.GetNContributors() - 3;
131  fQ = 0;
132  fAtProductionVertex = 0;
133  fSFromDecay = 0;
134 }
135 
136 Bool_t KFParticle::GetDistanceFromVertexXY( const float vtx[], const float Cv[], float &val, float &err ) const
137 {
146  Bool_t ret = 0;
147 
148  float mP[8];
149  float mC[36];
150  float dsdr[6] = {0.f};
151  const float dS = GetDStoPoint(vtx, dsdr);
152  Transport( dS, dsdr, mP, mC );
153 
154  float dx = mP[0] - vtx[0];
155  float dy = mP[1] - vtx[1];
156  float px = mP[3];
157  float py = mP[4];
158  float pt = sqrt(px*px + py*py);
159  float ex=0, ey=0;
160  if( pt<1.e-4 ){
161  ret = 1;
162  pt = 1.;
163  val = 1.e4;
164  } else{
165  ex = px/pt;
166  ey = py/pt;
167  val = dy*ex - dx*ey;
168  }
169 
170  float h0 = -ey;
171  float h1 = ex;
172  float h3 = (dy*ey + dx*ex)*ey/pt;
173  float h4 = -(dy*ey + dx*ex)*ex/pt;
174 
175  err =
176  h0*(h0*GetCovariance(0,0) + h1*GetCovariance(0,1) + h3*GetCovariance(0,3) + h4*GetCovariance(0,4) ) +
177  h1*(h0*GetCovariance(1,0) + h1*GetCovariance(1,1) + h3*GetCovariance(1,3) + h4*GetCovariance(1,4) ) +
178  h3*(h0*GetCovariance(3,0) + h1*GetCovariance(3,1) + h3*GetCovariance(3,3) + h4*GetCovariance(3,4) ) +
179  h4*(h0*GetCovariance(4,0) + h1*GetCovariance(4,1) + h3*GetCovariance(4,3) + h4*GetCovariance(4,4) );
180 
181  if( Cv ){
182  err+= h0*(h0*Cv[0] + h1*Cv[1] ) + h1*(h0*Cv[1] + h1*Cv[2] );
183  }
184 
185  err = sqrt(fabs(err));
186 
187  return ret;
188 }
189 
190 Bool_t KFParticle::GetDistanceFromVertexXY( const float vtx[], float &val, float &err ) const
191 {
198  return GetDistanceFromVertexXY( vtx, 0, val, err );
199 }
200 
201 
202 Bool_t KFParticle::GetDistanceFromVertexXY( const KFParticle &Vtx, float &val, float &err ) const
203 {
211  return GetDistanceFromVertexXY( Vtx.fP, Vtx.fC, val, err );
212 }
213 
214 #ifdef HomogeneousField
215 Bool_t KFParticle::GetDistanceFromVertexXY( const KFPVertex &Vtx, float &val, float &err ) const
216 {
224  return GetDistanceFromVertexXY( KFParticle(Vtx), val, err );
225 }
226 #endif
227 
228 float KFParticle::GetDistanceFromVertexXY( const float vtx[] ) const
229 {
234  float val, err;
235  GetDistanceFromVertexXY( vtx, 0, val, err );
236  return val;
237 }
238 
239 float KFParticle::GetDistanceFromVertexXY( const KFParticle &Vtx ) const
240 {
245  return GetDistanceFromVertexXY( Vtx.fP );
246 }
247 
248 #ifdef HomogeneousField
249 float KFParticle::GetDistanceFromVertexXY( const KFPVertex &Vtx ) const
250 {
255  return GetDistanceFromVertexXY( KFParticle(Vtx).fP );
256 }
257 #endif
258 
259 float KFParticle::GetDistanceFromParticleXY( const KFParticle &p ) const
260 {
265  float dsdr[4][6];
266  float dS[2];
267  GetDStoParticle( p, dS, dsdr );
268  float mP[8], mC[36], mP1[8], mC1[36];
269  Transport( dS[0], dsdr[0], mP, mC );
270  p.Transport( dS[1], dsdr[3], mP1, mC1 );
271  float dx = mP[0]-mP1[0];
272  float dy = mP[1]-mP1[1];
273  return sqrt(dx*dx+dy*dy);
274 }
275 
276 float KFParticle::GetDeviationFromParticleXY( const KFParticle &p ) const
277 {
282  float dsdr[4][6];
283  float dS[2];
284  GetDStoParticle( p, dS, dsdr );
285  float mP[8], mC[36], mP1[8], mC1[36];
286  Transport( dS[0], dsdr[0], mP, mC );
287  p.Transport( dS[1], dsdr[3], mP1, mC1 );
288 
289  float d[2]={ mP[0]-mP1[0], mP[1]-mP1[1] };
290 
291  float sigmaS = .1+10.*sqrt( (d[0]*d[0]+d[1]*d[1] )/
292  (mP1[3]*mP1[3]+mP1[4]*mP1[4] ) );
293 
294  float h[2] = { mP1[3]*sigmaS, mP1[4]*sigmaS };
295 
296  mC1[0] +=h[0]*h[0];
297  mC1[1] +=h[1]*h[0];
298  mC1[2] +=h[1]*h[1];
299 
300  return GetDeviationFromVertexXY( mP1, mC1 )*sqrt(2./1.);
301 }
302 
303 
304 float KFParticle::GetDeviationFromVertexXY( const float vtx[], const float Cv[] ) const
305 {
311  float val, err;
312  Bool_t problem = GetDistanceFromVertexXY( vtx, Cv, val, err );
313  if( problem || err<1.e-20 ) return 1.e4;
314  else return val/err;
315 }
316 
317 
318 float KFParticle::GetDeviationFromVertexXY( const KFParticle &Vtx ) const
319 {
324  return GetDeviationFromVertexXY( Vtx.fP, Vtx.fC );
325 }
326 
327 #ifdef HomogeneousField
328 float KFParticle::GetDeviationFromVertexXY( const KFPVertex &Vtx ) const
329 {
334  KFParticle v(Vtx);
335  return GetDeviationFromVertexXY( v.fP, v.fC );
336 }
337 #endif
338 
339 void KFParticle::GetParametersAtPoint(const float* point, const float* pointCov, float* m, float* mV)
340 {
348  float dsdr[6] = {0.f, 0.f, 0.f, 0.f, 0.f, 0.f};
349  float dS = GetDStoPoint(point, dsdr);
350  float dsdp[6] = {-dsdr[0], -dsdr[1], -dsdr[2], 0, 0, 0};
351 
352  float F[36], F1[36];
353  for(int i2=0; i2<36; i2++){
354  mV[i2] = 0.f;
355  F[i2] = 0.f;
356  F1[i2] = 0.f;
357  }
358  Transport(dS, dsdr, m, mV, dsdp, F, F1);
359 
360  float V1Tmp[36];
361  for(int i=0; i<36; i++)
362  V1Tmp[i] = 0.f;
363  KFParticle::MultQSQt(F1, pointCov, V1Tmp, 6);
364 
365  for(int iC=0; iC<21; iC++)
366  mV[iC] += V1Tmp[iC];
367 }
368 
369 float KFParticle::GetAngle ( const KFParticle &p ) const
370 {
375  float dsdr[4][6];
376  float dS[2];
377  GetDStoParticle( p, dS, dsdr );
378  float mP[8], mC[36], mP1[8], mC1[36];
379  Transport( dS[0], dsdr[0], mP, mC );
380  p.Transport( dS[1], dsdr[3], mP1, mC1 );
381  float n = sqrt( mP[3]*mP[3] + mP[4]*mP[4] + mP[5]*mP[5] );
382  float n1= sqrt( mP1[3]*mP1[3] + mP1[4]*mP1[4] + mP1[5]*mP1[5] );
383  n*=n1;
384  float a = 0;
385  if( n>1.e-8 ) a = ( mP[3]*mP1[3] + mP[4]*mP1[4] + mP[5]*mP1[5] )/n;
386  if (fabs(a)<1.) a = acos(a);
387  else a = (a>=0) ?0 :3.14;
388  return a;
389 }
390 
391 float KFParticle::GetAngleXY( const KFParticle &p ) const
392 {
397  float dsdr[4][6];
398  float dS[2];
399  GetDStoParticle( p, dS, dsdr );
400  float mP[8], mC[36], mP1[8], mC1[36];
401  Transport( dS[0], dsdr[0], mP, mC );
402  p.Transport( dS[1], dsdr[3], mP1, mC1 );
403  float n = sqrt( mP[3]*mP[3] + mP[4]*mP[4] );
404  float n1= sqrt( mP1[3]*mP1[3] + mP1[4]*mP1[4] );
405  n*=n1;
406  float a = 0;
407  if( n>1.e-8 ) a = ( mP[3]*mP1[3] + mP[4]*mP1[4] )/n;
408  if (fabs(a)<1.) a = acos(a);
409  else a = (a>=0) ?0 :3.14;
410  return a;
411 }
412 
413 float KFParticle::GetAngleRZ( const KFParticle &p ) const
414 {
419  float dsdr[4][6];
420  float dS[2];
421  GetDStoParticle( p, dS, dsdr );
422  float mP[8], mC[36], mP1[8], mC1[36];
423  Transport( dS[0], dsdr[0], mP, mC );
424  p.Transport( dS[1], dsdr[3], mP1, mC1 );
425  float nr = sqrt( mP[3]*mP[3] + mP[4]*mP[4] );
426  float n1r= sqrt( mP1[3]*mP1[3] + mP1[4]*mP1[4] );
427  float n = sqrt( nr*nr + mP[5]*mP[5] );
428  float n1= sqrt( n1r*n1r + mP1[5]*mP1[5] );
429  n*=n1;
430  float a = 0;
431  if( n>1.e-8 ) a = ( nr*n1r +mP[5]*mP1[5])/n;
432  if (fabs(a)<1.) a = acos(a);
433  else a = (a>=0) ?0 :3.14;
434  return a;
435 }
436 
437 float KFParticle::GetPseudoProperDecayTime( const KFParticle &pV, const float& mass, float* timeErr2 ) const
438 {
445  const float ipt2 = 1/( Px()*Px() + Py()*Py() );
446  const float mipt2 = mass*ipt2;
447  const float dx = X() - pV.X();
448  const float dy = Y() - pV.Y();
449 
450  if ( timeErr2 ) {
451  // -- calculate error = sigma(f(r)) = f'Cf'
452  // r = {x,y,px,py,x_pV,y_pV}
453  // df/dr = { px*m/pt^2,
454  // py*m/pt^2,
455  // ( x - x_pV )*m*(1/pt^2 - 2(px/pt^2)^2),
456  // ( y - y_pV )*m*(1/pt^2 - 2(py/pt^2)^2),
457  // -px*m/pt^2,
458  // -py*m/pt^2 }
459  const float f0 = Px()*mipt2;
460  const float f1 = Py()*mipt2;
461  const float mipt2derivative = mipt2*(1-2*Px()*Px()*ipt2);
462  const float f2 = dx*mipt2derivative;
463  const float f3 = -dy*mipt2derivative;
464  const float f4 = -f0;
465  const float f5 = -f1;
466 
467  const float& mC00 = GetCovariance(0,0);
468  const float& mC10 = GetCovariance(0,1);
469  const float& mC11 = GetCovariance(1,1);
470  const float& mC20 = GetCovariance(3,0);
471  const float& mC21 = GetCovariance(3,1);
472  const float& mC22 = GetCovariance(3,3);
473  const float& mC30 = GetCovariance(4,0);
474  const float& mC31 = GetCovariance(4,1);
475  const float& mC32 = GetCovariance(4,3);
476  const float& mC33 = GetCovariance(4,4);
477  const float& mC44 = pV.GetCovariance(0,0);
478  const float& mC54 = pV.GetCovariance(1,0);
479  const float& mC55 = pV.GetCovariance(1,1);
480 
481  *timeErr2 =
482  f5*mC55*f5 +
483  f5*mC54*f4 +
484  f4*mC44*f4 +
485  f3*mC33*f3 +
486  f3*mC32*f2 +
487  f3*mC31*f1 +
488  f3*mC30*f0 +
489  f2*mC22*f2 +
490  f2*mC21*f1 +
491  f2*mC20*f0 +
492  f1*mC11*f1 +
493  f1*mC10*f0 +
494  f0*mC00*f0;
495  }
496  return ( dx*Px() + dy*Py() )*mipt2;
497 }