Line data Source code
1 : //---------------------------------------------------------------------------------
2 : // Implementation of the AliKFParticleBase class
3 : // .
4 : // @author S.Gorbunov, I.Kisel, I.Kulakov, M.Zyzak
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 describes general mathematics which is used by AliKFParticle class
14 : //
15 : // -= Copyright © ALICE HLT Group =-
16 : //_________________________________________________________________________________
17 :
18 :
19 : #include "AliKFParticleBase.h"
20 : #include "TMath.h"
21 :
22 : #include <iostream>
23 172 : ClassImp(AliKFParticleBase)
24 :
25 :
26 604 : AliKFParticleBase::AliKFParticleBase() :fQ(0), fNDF(-3), fChi2(0), fSFromDecay(0), fAtProductionVertex(0), fIsLinearized(0),
27 302 : fConstructMethod(2), SumDaughterMass(0), fMassHypo(-1)
28 906 : {
29 : //* Constructor
30 :
31 302 : Initialize();
32 302 : }
33 :
34 : void AliKFParticleBase::Initialize( const Double_t Param[], const Double_t Cov[], Int_t Charge, Double_t Mass )
35 : {
36 : // Constructor from "cartesian" track, particle mass hypothesis should be provided
37 : //
38 : // Param[6] = { X, Y, Z, Px, Py, Pz } - position and momentum
39 : // Cov [21] = lower-triangular part of the covariance matrix:
40 : //
41 : // ( 0 . . . . . )
42 : // ( 1 2 . . . . )
43 : // Cov. matrix = ( 3 4 5 . . . ) - numbering of covariance elements in Cov[]
44 : // ( 6 7 8 9 . . )
45 : // ( 10 11 12 13 14 . )
46 : // ( 15 16 17 18 19 20 )
47 :
48 :
49 1680 : for( Int_t i=0; i<6 ; i++ ) fP[i] = Param[i];
50 4928 : for( Int_t i=0; i<21; i++ ) fC[i] = Cov[i];
51 :
52 112 : Double_t energy = TMath::Sqrt( Mass*Mass + fP[3]*fP[3] + fP[4]*fP[4] + fP[5]*fP[5]);
53 112 : fP[6] = energy;
54 112 : fP[7] = 0;
55 112 : fQ = Charge;
56 112 : fNDF = 0;
57 112 : fChi2 = 0;
58 112 : fAtProductionVertex = 0;
59 112 : fIsLinearized = 0;
60 112 : fSFromDecay = 0;
61 :
62 112 : Double_t energyInv = 1./energy;
63 : Double_t
64 112 : h0 = fP[3]*energyInv,
65 112 : h1 = fP[4]*energyInv,
66 112 : h2 = fP[5]*energyInv;
67 :
68 112 : fC[21] = h0*fC[ 6] + h1*fC[10] + h2*fC[15];
69 112 : fC[22] = h0*fC[ 7] + h1*fC[11] + h2*fC[16];
70 112 : fC[23] = h0*fC[ 8] + h1*fC[12] + h2*fC[17];
71 112 : fC[24] = h0*fC[ 9] + h1*fC[13] + h2*fC[18];
72 112 : fC[25] = h0*fC[13] + h1*fC[14] + h2*fC[19];
73 112 : fC[26] = h0*fC[18] + h1*fC[19] + h2*fC[20];
74 224 : fC[27] = ( h0*h0*fC[ 9] + h1*h1*fC[14] + h2*h2*fC[20]
75 112 : + 2*(h0*h1*fC[13] + h0*h2*fC[18] + h1*h2*fC[19] ) );
76 2016 : for( Int_t i=28; i<36; i++ ) fC[i] = 0;
77 112 : fC[35] = 1.;
78 :
79 112 : SumDaughterMass = Mass;
80 112 : fMassHypo = Mass;
81 112 : }
82 :
83 : void AliKFParticleBase::Initialize()
84 : {
85 : //* Initialise covariance matrix and set current parameters to 0.0
86 :
87 6080 : for( Int_t i=0; i<8; i++) fP[i] = 0;
88 23680 : for(Int_t i=0;i<36;++i) fC[i]=0.;
89 320 : fC[0] = fC[2] = fC[5] = 100.;
90 320 : fC[35] = 1.;
91 320 : fNDF = -3;
92 320 : fChi2 = 0.;
93 320 : fQ = 0;
94 320 : fSFromDecay = 0;
95 320 : fAtProductionVertex = 0;
96 320 : fVtxGuess[0]=fVtxGuess[1]=fVtxGuess[2]=0.;
97 320 : fIsLinearized = 0;
98 320 : SumDaughterMass = 0;
99 320 : fMassHypo = -1;
100 320 : }
101 :
102 : void AliKFParticleBase::SetVtxGuess( Double_t x, Double_t y, Double_t z )
103 : {
104 : //* Set decay vertex parameters for linearisation
105 :
106 0 : fVtxGuess[0] = x;
107 0 : fVtxGuess[1] = y;
108 0 : fVtxGuess[2] = z;
109 0 : fIsLinearized = 1;
110 0 : }
111 :
112 : Int_t AliKFParticleBase::GetMomentum( Double_t &p, Double_t &error ) const
113 : {
114 : //* Calculate particle momentum
115 :
116 0 : Double_t x = fP[3];
117 0 : Double_t y = fP[4];
118 0 : Double_t z = fP[5];
119 :
120 0 : Double_t x2 = x*x;
121 0 : Double_t y2 = y*y;
122 0 : Double_t z2 = z*z;
123 0 : Double_t p2 = x2+y2+z2;
124 0 : p = TMath::Sqrt(p2);
125 :
126 0 : error = (x2*fC[9]+y2*fC[14]+z2*fC[20] + 2*(x*y*fC[13]+x*z*fC[18]+y*z*fC[19]) );
127 0 : if( error>1.e-16 && p>1.e-4 ){
128 0 : error = TMath::Sqrt(error)/p;
129 0 : return 0;
130 : }
131 0 : error = 1.e8;
132 0 : return 1;
133 0 : }
134 :
135 : Int_t AliKFParticleBase::GetPt( Double_t &pt, Double_t &error ) const
136 : {
137 : //* Calculate particle transverse momentum
138 :
139 0 : Double_t px = fP[3];
140 0 : Double_t py = fP[4];
141 0 : Double_t px2 = px*px;
142 0 : Double_t py2 = py*py;
143 0 : Double_t pt2 = px2+py2;
144 0 : pt = TMath::Sqrt(pt2);
145 0 : error = (px2*fC[9] + py2*fC[14] + 2*px*py*fC[13] );
146 0 : if( error>0 && pt>1.e-4 ){
147 0 : error = TMath::Sqrt(error)/pt;
148 0 : return 0;
149 : }
150 0 : error = 1.e10;
151 0 : return 1;
152 0 : }
153 :
154 : Int_t AliKFParticleBase::GetEta( Double_t &eta, Double_t &error ) const
155 : {
156 : //* Calculate particle pseudorapidity
157 :
158 0 : Double_t px = fP[3];
159 0 : Double_t py = fP[4];
160 0 : Double_t pz = fP[5];
161 0 : Double_t pt2 = px*px + py*py;
162 0 : Double_t p2 = pt2 + pz*pz;
163 0 : Double_t p = TMath::Sqrt(p2);
164 0 : Double_t a = p + pz;
165 0 : Double_t b = p - pz;
166 0 : eta = 1.e10;
167 0 : if( b > 1.e-8 ){
168 0 : Double_t c = a/b;
169 0 : if( c>1.e-8 ) eta = 0.5*TMath::Log(c);
170 0 : }
171 0 : Double_t h3 = -px*pz;
172 0 : Double_t h4 = -py*pz;
173 0 : Double_t pt4 = pt2*pt2;
174 0 : Double_t p2pt4 = p2*pt4;
175 0 : error = (h3*h3*fC[9] + h4*h4*fC[14] + pt4*fC[20] + 2*( h3*(h4*fC[13] + fC[18]*pt2) + pt2*h4*fC[19] ) );
176 :
177 0 : if( error>0 && p2pt4>1.e-10 ){
178 0 : error = TMath::Sqrt(error/p2pt4);
179 0 : return 0;
180 : }
181 :
182 0 : error = 1.e10;
183 0 : return 1;
184 0 : }
185 :
186 : Int_t AliKFParticleBase::GetPhi( Double_t &phi, Double_t &error ) const
187 : {
188 : //* Calculate particle polar angle
189 :
190 0 : Double_t px = fP[3];
191 0 : Double_t py = fP[4];
192 0 : Double_t px2 = px*px;
193 0 : Double_t py2 = py*py;
194 0 : Double_t pt2 = px2 + py2;
195 0 : phi = TMath::ATan2(py,px);
196 0 : error = (py2*fC[9] + px2*fC[14] - 2*px*py*fC[13] );
197 0 : if( error>0 && pt2>1.e-4 ){
198 0 : error = TMath::Sqrt(error)/pt2;
199 0 : return 0;
200 : }
201 0 : error = 1.e10;
202 0 : return 1;
203 0 : }
204 :
205 : Int_t AliKFParticleBase::GetR( Double_t &r, Double_t &error ) const
206 : {
207 : //* Calculate distance to the origin
208 :
209 0 : Double_t x = fP[0];
210 0 : Double_t y = fP[1];
211 0 : Double_t x2 = x*x;
212 0 : Double_t y2 = y*y;
213 0 : r = TMath::Sqrt(x2 + y2);
214 0 : error = (x2*fC[0] + y2*fC[2] - 2*x*y*fC[1] );
215 0 : if( error>0 && r>1.e-4 ){
216 0 : error = TMath::Sqrt(error)/r;
217 0 : return 0;
218 : }
219 0 : error = 1.e10;
220 0 : return 1;
221 0 : }
222 :
223 : Int_t AliKFParticleBase::GetMass( Double_t &m, Double_t &error ) const
224 : {
225 : //* Calculate particle mass
226 :
227 : // s = sigma^2 of m2/2
228 :
229 36 : Double_t s = ( fP[3]*fP[3]*fC[9] + fP[4]*fP[4]*fC[14] + fP[5]*fP[5]*fC[20]
230 18 : + fP[6]*fP[6]*fC[27]
231 36 : +2*( + fP[3]*fP[4]*fC[13] + fP[5]*(fP[3]*fC[18] + fP[4]*fC[19])
232 18 : - fP[6]*( fP[3]*fC[24] + fP[4]*fC[25] + fP[5]*fC[26] ) )
233 : );
234 : // Double_t m2 = TMath::Abs(fP[6]*fP[6] - fP[3]*fP[3] - fP[4]*fP[4] - fP[5]*fP[5]);
235 : // m = TMath::Sqrt(m2);
236 : // if( m>1.e-10 ){
237 : // if( s>=0 ){
238 : // error = TMath::Sqrt(s)/m;
239 : // return 0;
240 : // }
241 : // }
242 : // error = 1.e20;
243 : // return 1;
244 18 : Double_t m2 = (fP[6]*fP[6] - fP[3]*fP[3] - fP[4]*fP[4] - fP[5]*fP[5]);
245 :
246 18 : if(m2<0.)
247 : {
248 4 : error = 1.e20;
249 4 : m = -TMath::Sqrt(-m2);
250 4 : return 1;
251 : }
252 :
253 14 : m = TMath::Sqrt(m2);
254 14 : if( m>1.e-6 ){
255 14 : if( s >= 0 ) {
256 14 : error = TMath::Sqrt(s)/m;
257 14 : return 0;
258 : }
259 : }
260 : else {
261 0 : error = 1.e20;
262 0 : return 0;
263 : }
264 0 : error = 1.e20;
265 :
266 0 : return 1;
267 18 : }
268 :
269 :
270 : Int_t AliKFParticleBase::GetDecayLength( Double_t &l, Double_t &error ) const
271 : {
272 : //* Calculate particle decay length [cm]
273 :
274 0 : Double_t x = fP[3];
275 0 : Double_t y = fP[4];
276 0 : Double_t z = fP[5];
277 0 : Double_t t = fP[7];
278 0 : Double_t x2 = x*x;
279 0 : Double_t y2 = y*y;
280 0 : Double_t z2 = z*z;
281 0 : Double_t p2 = x2+y2+z2;
282 0 : l = t*TMath::Sqrt(p2);
283 0 : if( p2>1.e-4){
284 0 : error = p2*fC[35] + t*t/p2*(x2*fC[9]+y2*fC[14]+z2*fC[20]
285 0 : + 2*(x*y*fC[13]+x*z*fC[18]+y*z*fC[19]) )
286 0 : + 2*t*(x*fC[31]+y*fC[32]+z*fC[33]);
287 0 : error = TMath::Sqrt(TMath::Abs(error));
288 0 : return 0;
289 : }
290 0 : error = 1.e20;
291 0 : return 1;
292 0 : }
293 :
294 : Int_t AliKFParticleBase::GetDecayLengthXY( Double_t &l, Double_t &error ) const
295 : {
296 : //* Calculate particle decay length in XY projection [cm]
297 :
298 0 : Double_t x = fP[3];
299 0 : Double_t y = fP[4];
300 0 : Double_t t = fP[7];
301 0 : Double_t x2 = x*x;
302 0 : Double_t y2 = y*y;
303 0 : Double_t pt2 = x2+y2;
304 0 : l = t*TMath::Sqrt(pt2);
305 0 : if( pt2>1.e-4){
306 0 : error = pt2*fC[35] + t*t/pt2*(x2*fC[9]+y2*fC[14] + 2*x*y*fC[13] )
307 0 : + 2*t*(x*fC[31]+y*fC[32]);
308 0 : error = TMath::Sqrt(TMath::Abs(error));
309 0 : return 0;
310 : }
311 0 : error = 1.e20;
312 0 : return 1;
313 0 : }
314 :
315 :
316 : Int_t AliKFParticleBase::GetLifeTime( Double_t &tauC, Double_t &error ) const
317 : {
318 : //* Calculate particle decay time [s]
319 :
320 0 : Double_t m, dm;
321 0 : GetMass( m, dm );
322 0 : Double_t cTM = (-fP[3]*fC[31] - fP[4]*fC[32] - fP[5]*fC[33] + fP[6]*fC[34]);
323 0 : tauC = fP[7]*m;
324 0 : error = m*m*fC[35] + 2*fP[7]*cTM + fP[7]*fP[7]*dm*dm;
325 0 : if( error > 0 ){
326 0 : error = TMath::Sqrt( error );
327 0 : return 0;
328 : }
329 0 : error = 1.e20;
330 0 : return 1;
331 0 : }
332 :
333 :
334 : void AliKFParticleBase::operator +=( const AliKFParticleBase &Daughter )
335 : {
336 : //* Add daughter via operator+=
337 :
338 296 : AddDaughter( Daughter );
339 148 : }
340 :
341 : Double_t AliKFParticleBase::GetSCorrection( const Double_t Part[], const Double_t XYZ[] )
342 : {
343 : //* Get big enough correction for S error to let the particle Part be fitted to XYZ point
344 :
345 888 : Double_t d[3] = { XYZ[0]-Part[0], XYZ[1]-Part[1], XYZ[2]-Part[2] };
346 444 : Double_t p2 = Part[3]*Part[3]+Part[4]*Part[4]+Part[5]*Part[5];
347 1332 : Double_t sigmaS = (p2>1.e-4) ? ( 10.1+3.*TMath::Sqrt( d[0]*d[0]+d[1]*d[1]+d[2]*d[2]) )/TMath::Sqrt(p2) : 1.;
348 444 : return sigmaS;
349 : }
350 :
351 : void AliKFParticleBase::GetMeasurement( const Double_t XYZ[], Double_t m[], Double_t V[] ) const
352 : {
353 : //* Get additional covariances V used during measurement
354 :
355 888 : Double_t b[3];
356 444 : GetFieldValue( XYZ, b );
357 : const Double_t kCLight = 0.000299792458;
358 444 : b[0]*=kCLight; b[1]*=kCLight; b[2]*=kCLight;
359 :
360 444 : Transport( GetDStoPoint(XYZ), m, V );
361 :
362 444 : Double_t sigmaS = GetSCorrection( m, XYZ );
363 :
364 : Double_t h[6];
365 :
366 444 : h[0] = m[3]*sigmaS;
367 444 : h[1] = m[4]*sigmaS;
368 444 : h[2] = m[5]*sigmaS;
369 444 : h[3] = ( h[1]*b[2]-h[2]*b[1] )*GetQ();
370 444 : h[4] = ( h[2]*b[0]-h[0]*b[2] )*GetQ();
371 444 : h[5] = ( h[0]*b[1]-h[1]*b[0] )*GetQ();
372 :
373 444 : V[ 0]+= h[0]*h[0];
374 444 : V[ 1]+= h[1]*h[0];
375 444 : V[ 2]+= h[1]*h[1];
376 444 : V[ 3]+= h[2]*h[0];
377 444 : V[ 4]+= h[2]*h[1];
378 444 : V[ 5]+= h[2]*h[2];
379 :
380 444 : V[ 6]+= h[3]*h[0];
381 444 : V[ 7]+= h[3]*h[1];
382 444 : V[ 8]+= h[3]*h[2];
383 444 : V[ 9]+= h[3]*h[3];
384 :
385 444 : V[10]+= h[4]*h[0];
386 444 : V[11]+= h[4]*h[1];
387 444 : V[12]+= h[4]*h[2];
388 444 : V[13]+= h[4]*h[3];
389 444 : V[14]+= h[4]*h[4];
390 :
391 444 : V[15]+= h[5]*h[0];
392 444 : V[16]+= h[5]*h[1];
393 444 : V[17]+= h[5]*h[2];
394 444 : V[18]+= h[5]*h[3];
395 444 : V[19]+= h[5]*h[4];
396 444 : V[20]+= h[5]*h[5];
397 444 : }
398 :
399 : void AliKFParticleBase::AddDaughter( const AliKFParticleBase &Daughter )
400 : {
401 296 : if( fNDF<-1 ){ // first daughter -> just copy
402 74 : fNDF = -1;
403 74 : fQ = Daughter.GetQ();
404 1184 : for( Int_t i=0; i<7; i++) fP[i] = Daughter.fP[i];
405 4292 : for( Int_t i=0; i<28; i++) fC[i] = Daughter.fC[i];
406 74 : fSFromDecay = 0;
407 74 : fMassHypo = Daughter.fMassHypo;
408 74 : SumDaughterMass = Daughter.SumDaughterMass;
409 74 : return;
410 : }
411 :
412 74 : if(fConstructMethod == 0)
413 0 : AddDaughterWithEnergyFit(Daughter);
414 74 : else if(fConstructMethod == 1)
415 0 : AddDaughterWithEnergyCalc(Daughter);
416 74 : else if(fConstructMethod == 2)
417 74 : AddDaughterWithEnergyFitMC(Daughter);
418 :
419 74 : SumDaughterMass += Daughter.SumDaughterMass;
420 74 : fMassHypo = -1;
421 222 : }
422 :
423 : void AliKFParticleBase::AddDaughterWithEnergyFit( const AliKFParticleBase &Daughter )
424 : {
425 : //* Energy considered as an independent veriable, fitted independently from momentum, without any constraints on mass
426 :
427 : //* Add daughter
428 :
429 0 : TransportToDecayVertex();
430 :
431 0 : Double_t b[3];
432 : Int_t maxIter = 1;
433 :
434 0 : if( !fIsLinearized ){
435 0 : if( fNDF==-1 ){
436 0 : Double_t ds, ds1;
437 0 : GetDStoParticle(Daughter, ds, ds1);
438 0 : TransportToDS( ds );
439 0 : Double_t m[8];
440 0 : Double_t mCd[36];
441 0 : Daughter.Transport( ds1, m, mCd );
442 0 : fVtxGuess[0] = .5*( fP[0] + m[0] );
443 0 : fVtxGuess[1] = .5*( fP[1] + m[1] );
444 0 : fVtxGuess[2] = .5*( fP[2] + m[2] );
445 0 : } else {
446 0 : fVtxGuess[0] = fP[0];
447 0 : fVtxGuess[1] = fP[1];
448 0 : fVtxGuess[2] = fP[2];
449 : }
450 : maxIter = 3;
451 0 : }
452 :
453 0 : for( Int_t iter=0; iter<maxIter; iter++ ){
454 :
455 : {
456 0 : GetFieldValue( fVtxGuess, b );
457 : const Double_t kCLight = 0.000299792458;
458 0 : b[0]*=kCLight; b[1]*=kCLight; b[2]*=kCLight;
459 : }
460 :
461 0 : Double_t *ffP = fP, *ffC = fC, tmpP[8], tmpC[36];
462 0 : if( fNDF==-1 ){
463 0 : GetMeasurement( fVtxGuess, tmpP, tmpC );
464 : ffP = tmpP;
465 : ffC = tmpC;
466 0 : }
467 :
468 0 : Double_t m[8], mV[36];
469 :
470 0 : if( Daughter.fC[35]>0 ){
471 0 : Daughter.GetMeasurement( fVtxGuess, m, mV );
472 0 : } else {
473 0 : for( Int_t i=0; i<8; i++ ) m[i] = Daughter.fP[i];
474 0 : for( Int_t i=0; i<36; i++ ) mV[i] = Daughter.fC[i];
475 : }
476 : //*
477 :
478 : Double_t mS[6];
479 : {
480 0 : Double_t mSi[6] = { ffC[0]+mV[0],
481 0 : ffC[1]+mV[1], ffC[2]+mV[2],
482 0 : ffC[3]+mV[3], ffC[4]+mV[4], ffC[5]+mV[5] };
483 :
484 0 : mS[0] = mSi[2]*mSi[5] - mSi[4]*mSi[4];
485 0 : mS[1] = mSi[3]*mSi[4] - mSi[1]*mSi[5];
486 0 : mS[2] = mSi[0]*mSi[5] - mSi[3]*mSi[3];
487 0 : mS[3] = mSi[1]*mSi[4] - mSi[2]*mSi[3];
488 0 : mS[4] = mSi[1]*mSi[3] - mSi[0]*mSi[4];
489 0 : mS[5] = mSi[0]*mSi[2] - mSi[1]*mSi[1];
490 :
491 0 : Double_t s = ( mSi[0]*mS[0] + mSi[1]*mS[1] + mSi[3]*mS[3] );
492 0 : s = ( TMath::Abs(s) > 1.E-20 ) ?1./s :0;
493 0 : mS[0]*=s;
494 0 : mS[1]*=s;
495 0 : mS[2]*=s;
496 0 : mS[3]*=s;
497 0 : mS[4]*=s;
498 0 : mS[5]*=s;
499 : }
500 : //* Residual (measured - estimated)
501 :
502 0 : Double_t zeta[3] = { m[0]-ffP[0], m[1]-ffP[1], m[2]-ffP[2] };
503 :
504 : //* CHt = CH' - D'
505 :
506 0 : Double_t mCHt0[7], mCHt1[7], mCHt2[7];
507 :
508 0 : mCHt0[0]=ffC[ 0] ; mCHt1[0]=ffC[ 1] ; mCHt2[0]=ffC[ 3] ;
509 0 : mCHt0[1]=ffC[ 1] ; mCHt1[1]=ffC[ 2] ; mCHt2[1]=ffC[ 4] ;
510 0 : mCHt0[2]=ffC[ 3] ; mCHt1[2]=ffC[ 4] ; mCHt2[2]=ffC[ 5] ;
511 0 : mCHt0[3]=ffC[ 6]-mV[ 6]; mCHt1[3]=ffC[ 7]-mV[ 7]; mCHt2[3]=ffC[ 8]-mV[ 8];
512 0 : mCHt0[4]=ffC[10]-mV[10]; mCHt1[4]=ffC[11]-mV[11]; mCHt2[4]=ffC[12]-mV[12];
513 0 : mCHt0[5]=ffC[15]-mV[15]; mCHt1[5]=ffC[16]-mV[16]; mCHt2[5]=ffC[17]-mV[17];
514 0 : mCHt0[6]=ffC[21]-mV[21]; mCHt1[6]=ffC[22]-mV[22]; mCHt2[6]=ffC[23]-mV[23];
515 :
516 : //* Kalman gain K = mCH'*S
517 :
518 0 : Double_t k0[7], k1[7], k2[7];
519 :
520 0 : for(Int_t i=0;i<7;++i){
521 0 : k0[i] = mCHt0[i]*mS[0] + mCHt1[i]*mS[1] + mCHt2[i]*mS[3];
522 0 : k1[i] = mCHt0[i]*mS[1] + mCHt1[i]*mS[2] + mCHt2[i]*mS[4];
523 0 : k2[i] = mCHt0[i]*mS[3] + mCHt1[i]*mS[4] + mCHt2[i]*mS[5];
524 : }
525 :
526 : //* New estimation of the vertex position
527 :
528 0 : if( iter<maxIter-1 ){
529 0 : for(Int_t i=0; i<3; ++i)
530 0 : fVtxGuess[i]= ffP[i] + k0[i]*zeta[0]+k1[i]*zeta[1]+k2[i]*zeta[2];
531 0 : continue;
532 : }
533 :
534 : // last itearation -> update the particle
535 :
536 : //* Add the daughter momentum to the particle momentum
537 :
538 0 : ffP[ 3] += m[ 3];
539 0 : ffP[ 4] += m[ 4];
540 0 : ffP[ 5] += m[ 5];
541 0 : ffP[ 6] += m[ 6];
542 :
543 0 : ffC[ 9] += mV[ 9];
544 0 : ffC[13] += mV[13];
545 0 : ffC[14] += mV[14];
546 0 : ffC[18] += mV[18];
547 0 : ffC[19] += mV[19];
548 0 : ffC[20] += mV[20];
549 0 : ffC[24] += mV[24];
550 0 : ffC[25] += mV[25];
551 0 : ffC[26] += mV[26];
552 0 : ffC[27] += mV[27];
553 :
554 :
555 : //* New estimation of the vertex position r += K*zeta
556 :
557 0 : for(Int_t i=0;i<7;++i)
558 0 : fP[i] = ffP[i] + k0[i]*zeta[0] + k1[i]*zeta[1] + k2[i]*zeta[2];
559 :
560 : //* New covariance matrix C -= K*(mCH')'
561 :
562 0 : for(Int_t i=0, k=0;i<7;++i){
563 0 : for(Int_t j=0;j<=i;++j,++k){
564 0 : fC[k] = ffC[k] - (k0[i]*mCHt0[j] + k1[i]*mCHt1[j] + k2[i]*mCHt2[j] );
565 : }
566 : }
567 :
568 : //* Calculate Chi^2
569 :
570 0 : fNDF += 2;
571 0 : fQ += Daughter.GetQ();
572 0 : fSFromDecay = 0;
573 0 : fChi2 += (mS[0]*zeta[0] + mS[1]*zeta[1] + mS[3]*zeta[2])*zeta[0]
574 0 : + (mS[1]*zeta[0] + mS[2]*zeta[1] + mS[4]*zeta[2])*zeta[1]
575 0 : + (mS[3]*zeta[0] + mS[4]*zeta[1] + mS[5]*zeta[2])*zeta[2];
576 :
577 0 : }
578 0 : }
579 :
580 : void AliKFParticleBase::AddDaughterWithEnergyCalc( const AliKFParticleBase &Daughter )
581 : {
582 : //* Energy considered as a dependent variable, calculated from the momentum and mass hypothesis
583 :
584 : //* Add daughter
585 :
586 0 : TransportToDecayVertex();
587 :
588 0 : Double_t b[3];
589 : Int_t maxIter = 1;
590 :
591 0 : if( !fIsLinearized ){
592 0 : if( fNDF==-1 ){
593 0 : Double_t ds, ds1;
594 0 : GetDStoParticle(Daughter, ds, ds1);
595 0 : TransportToDS( ds );
596 0 : Double_t m[8];
597 0 : Double_t mCd[36];
598 0 : Daughter.Transport( ds1, m, mCd );
599 0 : fVtxGuess[0] = .5*( fP[0] + m[0] );
600 0 : fVtxGuess[1] = .5*( fP[1] + m[1] );
601 0 : fVtxGuess[2] = .5*( fP[2] + m[2] );
602 0 : } else {
603 0 : fVtxGuess[0] = fP[0];
604 0 : fVtxGuess[1] = fP[1];
605 0 : fVtxGuess[2] = fP[2];
606 : }
607 : maxIter = 3;
608 0 : }
609 :
610 0 : for( Int_t iter=0; iter<maxIter; iter++ ){
611 :
612 : {
613 0 : GetFieldValue( fVtxGuess, b );
614 : const Double_t kCLight = 0.000299792458;
615 0 : b[0]*=kCLight; b[1]*=kCLight; b[2]*=kCLight;
616 : }
617 :
618 0 : Double_t *ffP = fP, *ffC = fC, tmpP[8], tmpC[36];
619 0 : if( fNDF==-1 ){
620 0 : GetMeasurement( fVtxGuess, tmpP, tmpC );
621 : ffP = tmpP;
622 : ffC = tmpC;
623 0 : }
624 :
625 0 : Double_t m[8], mV[36];
626 :
627 0 : if( Daughter.fC[35]>0 ){
628 0 : Daughter.GetMeasurement( fVtxGuess, m, mV );
629 0 : } else {
630 0 : for( Int_t i=0; i<8; i++ ) m[i] = Daughter.fP[i];
631 0 : for( Int_t i=0; i<36; i++ ) mV[i] = Daughter.fC[i];
632 : }
633 :
634 0 : double massMf2 = m[6]*m[6] - (m[3]*m[3] + m[4]*m[4] + m[5]*m[5]);
635 0 : double massRf2 = fP[6]*fP[6] - (fP[3]*fP[3] + fP[4]*fP[4] + fP[5]*fP[5]);
636 :
637 : //*
638 :
639 : Double_t mS[6];
640 : {
641 0 : Double_t mSi[6] = { ffC[0]+mV[0],
642 0 : ffC[1]+mV[1], ffC[2]+mV[2],
643 0 : ffC[3]+mV[3], ffC[4]+mV[4], ffC[5]+mV[5] };
644 :
645 0 : mS[0] = mSi[2]*mSi[5] - mSi[4]*mSi[4];
646 0 : mS[1] = mSi[3]*mSi[4] - mSi[1]*mSi[5];
647 0 : mS[2] = mSi[0]*mSi[5] - mSi[3]*mSi[3];
648 0 : mS[3] = mSi[1]*mSi[4] - mSi[2]*mSi[3];
649 0 : mS[4] = mSi[1]*mSi[3] - mSi[0]*mSi[4];
650 0 : mS[5] = mSi[0]*mSi[2] - mSi[1]*mSi[1];
651 :
652 0 : Double_t s = ( mSi[0]*mS[0] + mSi[1]*mS[1] + mSi[3]*mS[3] );
653 :
654 0 : s = ( s > 1.E-20 ) ?1./s :0;
655 0 : mS[0]*=s;
656 0 : mS[1]*=s;
657 0 : mS[2]*=s;
658 0 : mS[3]*=s;
659 0 : mS[4]*=s;
660 0 : mS[5]*=s;
661 : }
662 :
663 : //* Residual (measured - estimated)
664 :
665 0 : Double_t zeta[3] = { m[0]-ffP[0], m[1]-ffP[1], m[2]-ffP[2] };
666 :
667 : //* CHt = CH' - D'
668 :
669 0 : Double_t mCHt0[6], mCHt1[6], mCHt2[6];
670 :
671 0 : mCHt0[0]=ffC[ 0] ; mCHt1[0]=ffC[ 1] ; mCHt2[0]=ffC[ 3] ;
672 0 : mCHt0[1]=ffC[ 1] ; mCHt1[1]=ffC[ 2] ; mCHt2[1]=ffC[ 4] ;
673 0 : mCHt0[2]=ffC[ 3] ; mCHt1[2]=ffC[ 4] ; mCHt2[2]=ffC[ 5] ;
674 0 : mCHt0[3]=ffC[ 6]-mV[ 6]; mCHt1[3]=ffC[ 7]-mV[ 7]; mCHt2[3]=ffC[ 8]-mV[ 8];
675 0 : mCHt0[4]=ffC[10]-mV[10]; mCHt1[4]=ffC[11]-mV[11]; mCHt2[4]=ffC[12]-mV[12];
676 0 : mCHt0[5]=ffC[15]-mV[15]; mCHt1[5]=ffC[16]-mV[16]; mCHt2[5]=ffC[17]-mV[17];
677 :
678 : //* Kalman gain K = mCH'*S
679 :
680 0 : Double_t k0[6], k1[6], k2[6];
681 :
682 0 : for(Int_t i=0;i<6;++i){
683 0 : k0[i] = mCHt0[i]*mS[0] + mCHt1[i]*mS[1] + mCHt2[i]*mS[3];
684 0 : k1[i] = mCHt0[i]*mS[1] + mCHt1[i]*mS[2] + mCHt2[i]*mS[4];
685 0 : k2[i] = mCHt0[i]*mS[3] + mCHt1[i]*mS[4] + mCHt2[i]*mS[5];
686 : }
687 :
688 : //* New estimation of the vertex position
689 :
690 0 : if( iter<maxIter-1 ){
691 0 : for(Int_t i=0; i<3; ++i)
692 0 : fVtxGuess[i]= ffP[i] + k0[i]*zeta[0]+k1[i]*zeta[1]+k2[i]*zeta[2];
693 0 : continue;
694 : }
695 :
696 : //* find mf and mVf - optimum value of the measurement and its covariance matrix
697 : //* mVHt = V*H'
698 0 : Double_t mVHt0[6], mVHt1[6], mVHt2[6];
699 :
700 0 : mVHt0[0]= mV[ 0] ; mVHt1[0]= mV[ 1] ; mVHt2[0]= mV[ 3] ;
701 0 : mVHt0[1]= mV[ 1] ; mVHt1[1]= mV[ 2] ; mVHt2[1]= mV[ 4] ;
702 0 : mVHt0[2]= mV[ 3] ; mVHt1[2]= mV[ 4] ; mVHt2[2]= mV[ 5] ;
703 0 : mVHt0[3]= mV[ 6] ; mVHt1[3]= mV[ 7] ; mVHt2[3]= mV[ 8] ;
704 0 : mVHt0[4]= mV[10] ; mVHt1[4]= mV[11] ; mVHt2[4]= mV[12] ;
705 0 : mVHt0[5]= mV[15] ; mVHt1[5]= mV[16] ; mVHt2[5]= mV[17] ;
706 :
707 : //* Kalman gain Km = mCH'*S
708 :
709 0 : Double_t km0[6], km1[6], km2[6];
710 :
711 0 : for(Int_t i=0;i<6;++i){
712 0 : km0[i] = mVHt0[i]*mS[0] + mVHt1[i]*mS[1] + mVHt2[i]*mS[3];
713 0 : km1[i] = mVHt0[i]*mS[1] + mVHt1[i]*mS[2] + mVHt2[i]*mS[4];
714 0 : km2[i] = mVHt0[i]*mS[3] + mVHt1[i]*mS[4] + mVHt2[i]*mS[5];
715 : }
716 :
717 0 : Double_t mf[7] = { m[0], m[1], m[2], m[3], m[4], m[5], m[6] };
718 :
719 0 : for(Int_t i=0;i<6;++i)
720 0 : mf[i] = mf[i] - km0[i]*zeta[0] - km1[i]*zeta[1] - km2[i]*zeta[2];
721 :
722 0 : Double_t energyMf = TMath::Sqrt( massMf2 + (mf[3]*mf[3] + mf[4]*mf[4] + mf[5]*mf[5]) );
723 :
724 0 : Double_t mVf[28];
725 0 : for(Int_t iC=0; iC<28; iC++)
726 0 : mVf[iC] = mV[iC];
727 :
728 : //* hmf = d(energyMf)/d(mf)
729 : Double_t hmf[7];
730 0 : if( TMath::Abs(energyMf) < 1.e-10) hmf[3] = 0; else hmf[3] = mf[3]/energyMf;
731 0 : if( TMath::Abs(energyMf) < 1.e-10) hmf[4] = 0; else hmf[4] = mf[4]/energyMf;
732 0 : if( TMath::Abs(energyMf) < 1.e-10) hmf[5] = 0; else hmf[5] = mf[5]/energyMf;
733 : // if( TMath::Abs(energyMf) < 1.e-10) hmf[6] = 0; else hmf[6] = mf[6]/energyMf;
734 : hmf[6] = 0;
735 :
736 0 : for(Int_t i=0, k=0;i<6;++i){
737 0 : for(Int_t j=0;j<=i;++j,++k){
738 0 : mVf[k] = mVf[k] - (km0[i]*mVHt0[j] + km1[i]*mVHt1[j] + km2[i]*mVHt2[j] );
739 : }
740 : }
741 0 : Double_t mVf24 = mVf[24], mVf25 = mVf[25], mVf26 = mVf[26];
742 0 : mVf[21] = mVf[6 ]*hmf[3] + mVf[10]*hmf[4] + mVf[15]*hmf[5] + mVf[21]*hmf[6];
743 0 : mVf[22] = mVf[7 ]*hmf[3] + mVf[11]*hmf[4] + mVf[16]*hmf[5] + mVf[22]*hmf[6];
744 0 : mVf[23] = mVf[8 ]*hmf[3] + mVf[12]*hmf[4] + mVf[17]*hmf[5] + mVf[23]*hmf[6];
745 0 : mVf[24] = mVf[9 ]*hmf[3] + mVf[13]*hmf[4] + mVf[18]*hmf[5] + mVf[24]*hmf[6];
746 0 : mVf[25] = mVf[13]*hmf[3] + mVf[14]*hmf[4] + mVf[19]*hmf[5] + mVf[25]*hmf[6];
747 0 : mVf[26] = mVf[18]*hmf[3] + mVf[19]*hmf[4] + mVf[20]*hmf[5] + mVf[26]*hmf[6];
748 0 : mVf[27] = mVf[24]*hmf[3] + mVf[25]*hmf[4] + mVf[26]*hmf[5] + (mVf24*hmf[3] + mVf25*hmf[4] + mVf26*hmf[5] + mVf[27]*hmf[6])*hmf[6]; //here mVf[] are already modified
749 :
750 0 : mf[6] = energyMf;
751 :
752 : //* find rf and mCf - optimum value of the measurement and its covariance matrix
753 :
754 : //* mCCHt = C*H'
755 0 : Double_t mCCHt0[6], mCCHt1[6], mCCHt2[6];
756 :
757 0 : mCCHt0[0]=ffC[ 0]; mCCHt1[0]=ffC[ 1]; mCCHt2[0]=ffC[ 3];
758 0 : mCCHt0[1]=ffC[ 1]; mCCHt1[1]=ffC[ 2]; mCCHt2[1]=ffC[ 4];
759 0 : mCCHt0[2]=ffC[ 3]; mCCHt1[2]=ffC[ 4]; mCCHt2[2]=ffC[ 5];
760 0 : mCCHt0[3]=ffC[ 6]; mCCHt1[3]=ffC[ 7]; mCCHt2[3]=ffC[ 8];
761 0 : mCCHt0[4]=ffC[10]; mCCHt1[4]=ffC[11]; mCCHt2[4]=ffC[12];
762 0 : mCCHt0[5]=ffC[15]; mCCHt1[5]=ffC[16]; mCCHt2[5]=ffC[17];
763 :
764 : //* Kalman gain Krf = mCH'*S
765 :
766 0 : Double_t krf0[6], krf1[6], krf2[6];
767 :
768 0 : for(Int_t i=0;i<6;++i){
769 0 : krf0[i] = mCCHt0[i]*mS[0] + mCCHt1[i]*mS[1] + mCCHt2[i]*mS[3];
770 0 : krf1[i] = mCCHt0[i]*mS[1] + mCCHt1[i]*mS[2] + mCCHt2[i]*mS[4];
771 0 : krf2[i] = mCCHt0[i]*mS[3] + mCCHt1[i]*mS[4] + mCCHt2[i]*mS[5];
772 : }
773 0 : Double_t rf[7] = { ffP[0], ffP[1], ffP[2], ffP[3], ffP[4], ffP[5], ffP[6] };
774 :
775 0 : for(Int_t i=0;i<6;++i)
776 0 : rf[i] = rf[i] + krf0[i]*zeta[0] + krf1[i]*zeta[1] + krf2[i]*zeta[2];
777 :
778 0 : Double_t energyRf = TMath::Sqrt( massRf2 + (rf[3]*rf[3] + rf[4]*rf[4] + rf[5]*rf[5]) );
779 :
780 0 : Double_t mCf[28];
781 0 : for(Int_t iC=0; iC<28; iC++)
782 0 : mCf[iC] = ffC[iC];
783 : //* hrf = d(Erf)/d(rf)
784 : Double_t hrf[7];
785 0 : if( TMath::Abs(energyRf) < 1.e-10) hrf[3] = 0; else hrf[3] = rf[3]/energyRf;
786 0 : if( TMath::Abs(energyRf) < 1.e-10) hrf[4] = 0; else hrf[4] = rf[4]/energyRf;
787 0 : if( TMath::Abs(energyRf) < 1.e-10) hrf[5] = 0; else hrf[5] = rf[5]/energyRf;
788 : // if( TMath::Abs(energyRf) < 1.e-10) hrf[6] = 0; else hrf[6] = rf[6]/energyRf;
789 : hrf[6] = 0;
790 :
791 0 : for(Int_t i=0, k=0;i<6;++i){
792 0 : for(Int_t j=0;j<=i;++j,++k){
793 0 : mCf[k] = mCf[k] - (krf0[i]*mCCHt0[j] + krf1[i]*mCCHt1[j] + krf2[i]*mCCHt2[j] );
794 : }
795 : }
796 0 : Double_t mCf24 = mCf[24], mCf25 = mCf[25], mCf26 = mCf[26];
797 0 : mCf[21] = mCf[6 ]*hrf[3] + mCf[10]*hrf[4] + mCf[15]*hrf[5] + mCf[21]*hrf[6];
798 0 : mCf[22] = mCf[7 ]*hrf[3] + mCf[11]*hrf[4] + mCf[16]*hrf[5] + mCf[22]*hrf[6];
799 0 : mCf[23] = mCf[8 ]*hrf[3] + mCf[12]*hrf[4] + mCf[17]*hrf[5] + mCf[23]*hrf[6];
800 0 : mCf[24] = mCf[9 ]*hrf[3] + mCf[13]*hrf[4] + mCf[18]*hrf[5] + mCf[24]*hrf[6];
801 0 : mCf[25] = mCf[13]*hrf[3] + mCf[14]*hrf[4] + mCf[19]*hrf[5] + mCf[25]*hrf[6];
802 0 : mCf[26] = mCf[18]*hrf[3] + mCf[19]*hrf[4] + mCf[20]*hrf[5] + mCf[26]*hrf[6];
803 0 : mCf[27] = mCf[24]*hrf[3] + mCf[25]*hrf[4] + mCf[26]*hrf[5] + (mCf24*hrf[3] + mCf25*hrf[4] + mCf26*hrf[5] + mCf[27]*hrf[6])*hrf[6]; //here mCf[] are already modified
804 :
805 0 : for(Int_t iC=21; iC<28; iC++)
806 : {
807 0 : ffC[iC] = mCf[iC];
808 0 : mV[iC] = mVf[iC];
809 : }
810 :
811 0 : fP[6] = energyRf + energyMf;
812 0 : rf[6] = energyRf;
813 :
814 : //Double_t Dvv[3][3]; do not need this
815 0 : Double_t mDvp[3][3];
816 : // Double_t mDpv[3][3];
817 0 : Double_t mDpp[3][3];
818 : Double_t mDe[7];
819 :
820 0 : for(int i=0; i<3; i++)
821 : {
822 0 : for(int j=0; j<3; j++)
823 : {
824 0 : mDvp[i][j] = km0[i+3]*mCCHt0[j] + km1[i+3]*mCCHt1[j] + km2[i+3]*mCCHt2[j];
825 : // mDpv[i][j] = km0[i]*mCCHt0[j+3] + km1[i]*mCCHt1[j+3] + km2[i]*mCCHt2[j+3];
826 0 : mDpp[i][j] = km0[i+3]*mCCHt0[j+3] + km1[i+3]*mCCHt1[j+3] + km2[i+3]*mCCHt2[j+3];
827 : }
828 : }
829 :
830 0 : mDe[0] = hmf[3]*mDvp[0][0] + hmf[4]*mDvp[1][0] + hmf[5]*mDvp[2][0];
831 0 : mDe[1] = hmf[3]*mDvp[0][1] + hmf[4]*mDvp[1][1] + hmf[5]*mDvp[2][1];
832 0 : mDe[2] = hmf[3]*mDvp[0][2] + hmf[4]*mDvp[1][2] + hmf[5]*mDvp[2][2];
833 0 : mDe[3] = hmf[3]*mDpp[0][0] + hmf[4]*mDpp[1][0] + hmf[5]*mDpp[2][0];
834 0 : mDe[4] = hmf[3]*mDpp[0][1] + hmf[4]*mDpp[1][1] + hmf[5]*mDpp[2][1];
835 0 : mDe[5] = hmf[3]*mDpp[0][2] + hmf[4]*mDpp[1][2] + hmf[5]*mDpp[2][2];
836 0 : mDe[6] = 2*(mDe[3]*hrf[3] + mDe[4]*hrf[4] + mDe[5]*hrf[5]);
837 :
838 : // last itearation -> update the particle
839 :
840 : //* Add the daughter momentum to the particle momentum
841 :
842 0 : ffP[ 3] += m[ 3];
843 0 : ffP[ 4] += m[ 4];
844 0 : ffP[ 5] += m[ 5];
845 :
846 0 : ffC[ 9] += mV[ 9];
847 0 : ffC[13] += mV[13];
848 0 : ffC[14] += mV[14];
849 0 : ffC[18] += mV[18];
850 0 : ffC[19] += mV[19];
851 0 : ffC[20] += mV[20];
852 0 : ffC[24] += mV[24];
853 0 : ffC[25] += mV[25];
854 0 : ffC[26] += mV[26];
855 0 : ffC[27] += mV[27];
856 :
857 0 : ffC[21] += mDe[0];
858 0 : ffC[22] += mDe[1];
859 0 : ffC[23] += mDe[2];
860 0 : ffC[24] += mDe[3];
861 0 : ffC[25] += mDe[4];
862 0 : ffC[26] += mDe[5];
863 0 : ffC[27] += mDe[6];
864 :
865 : //* New estimation of the vertex position r += K*zeta
866 :
867 0 : for(Int_t i=0;i<6;++i)
868 0 : fP[i] = ffP[i] + k0[i]*zeta[0] + k1[i]*zeta[1] + k2[i]*zeta[2];
869 :
870 : //* New covariance matrix C -= K*(mCH')'
871 :
872 0 : for(Int_t i=0, k=0;i<6;++i){
873 0 : for(Int_t j=0;j<=i;++j,++k){
874 0 : fC[k] = ffC[k] - (k0[i]*mCHt0[j] + k1[i]*mCHt1[j] + k2[i]*mCHt2[j] );
875 : }
876 : }
877 :
878 0 : for(int i=21; i<28; i++) fC[i] = ffC[i];
879 :
880 : //* Calculate Chi^2
881 :
882 0 : fNDF += 2;
883 0 : fQ += Daughter.GetQ();
884 0 : fSFromDecay = 0;
885 0 : fChi2 += (mS[0]*zeta[0] + mS[1]*zeta[1] + mS[3]*zeta[2])*zeta[0]
886 0 : + (mS[1]*zeta[0] + mS[2]*zeta[1] + mS[4]*zeta[2])*zeta[1]
887 0 : + (mS[3]*zeta[0] + mS[4]*zeta[1] + mS[5]*zeta[2])*zeta[2];
888 0 : }
889 0 : }
890 :
891 : void AliKFParticleBase::AddDaughterWithEnergyFitMC( const AliKFParticleBase &Daughter )
892 : {
893 : //* Energy considered as an independent variable, fitted independently from momentum, without any constraints on mass
894 :
895 : //* Add daughter
896 :
897 148 : TransportToDecayVertex();
898 :
899 74 : Double_t b[3];
900 : Int_t maxIter = 1;
901 :
902 74 : if( !fIsLinearized ){
903 74 : if( fNDF==-1 ){
904 74 : Double_t ds, ds1;
905 74 : GetDStoParticle(Daughter, ds, ds1);
906 74 : TransportToDS( ds );
907 74 : Double_t m[8];
908 74 : Double_t mCd[36];
909 74 : Daughter.Transport( ds1, m, mCd );
910 74 : fVtxGuess[0] = .5*( fP[0] + m[0] );
911 74 : fVtxGuess[1] = .5*( fP[1] + m[1] );
912 74 : fVtxGuess[2] = .5*( fP[2] + m[2] );
913 74 : } else {
914 0 : fVtxGuess[0] = fP[0];
915 0 : fVtxGuess[1] = fP[1];
916 0 : fVtxGuess[2] = fP[2];
917 : }
918 : maxIter = 3;
919 74 : }
920 :
921 592 : for( Int_t iter=0; iter<maxIter; iter++ ){
922 :
923 : {
924 222 : GetFieldValue( fVtxGuess, b );
925 : const Double_t kCLight = 0.000299792458;
926 222 : b[0]*=kCLight; b[1]*=kCLight; b[2]*=kCLight;
927 : }
928 :
929 222 : Double_t *ffP = fP, *ffC = fC, tmpP[8], tmpC[36];
930 222 : if( fNDF==-1 ){
931 222 : GetMeasurement( fVtxGuess, tmpP, tmpC );
932 : ffP = tmpP;
933 : ffC = tmpC;
934 222 : }
935 222 : Double_t m[8], mV[36];
936 :
937 222 : if( Daughter.fC[35]>0 ){
938 222 : Daughter.GetMeasurement( fVtxGuess, m, mV );
939 222 : } else {
940 0 : for( Int_t i=0; i<8; i++ ) m[i] = Daughter.fP[i];
941 0 : for( Int_t i=0; i<36; i++ ) mV[i] = Daughter.fC[i];
942 : }
943 : //*
944 :
945 : Double_t mS[6];
946 : {
947 222 : Double_t mSi[6] = { ffC[0]+mV[0],
948 222 : ffC[1]+mV[1], ffC[2]+mV[2],
949 222 : ffC[3]+mV[3], ffC[4]+mV[4], ffC[5]+mV[5] };
950 :
951 222 : mS[0] = mSi[2]*mSi[5] - mSi[4]*mSi[4];
952 222 : mS[1] = mSi[3]*mSi[4] - mSi[1]*mSi[5];
953 222 : mS[2] = mSi[0]*mSi[5] - mSi[3]*mSi[3];
954 222 : mS[3] = mSi[1]*mSi[4] - mSi[2]*mSi[3];
955 222 : mS[4] = mSi[1]*mSi[3] - mSi[0]*mSi[4];
956 222 : mS[5] = mSi[0]*mSi[2] - mSi[1]*mSi[1];
957 :
958 222 : Double_t s = ( mSi[0]*mS[0] + mSi[1]*mS[1] + mSi[3]*mS[3] );
959 :
960 666 : s = ( s > 1.E-20 ) ?1./s :0;
961 222 : mS[0]*=s;
962 222 : mS[1]*=s;
963 222 : mS[2]*=s;
964 222 : mS[3]*=s;
965 222 : mS[4]*=s;
966 222 : mS[5]*=s;
967 : }
968 : //* Residual (measured - estimated)
969 :
970 222 : Double_t zeta[3] = { m[0]-ffP[0], m[1]-ffP[1], m[2]-ffP[2] };
971 :
972 :
973 : //* CHt = CH'
974 :
975 222 : Double_t mCHt0[7], mCHt1[7], mCHt2[7];
976 :
977 222 : mCHt0[0]=ffC[ 0] ; mCHt1[0]=ffC[ 1] ; mCHt2[0]=ffC[ 3] ;
978 222 : mCHt0[1]=ffC[ 1] ; mCHt1[1]=ffC[ 2] ; mCHt2[1]=ffC[ 4] ;
979 222 : mCHt0[2]=ffC[ 3] ; mCHt1[2]=ffC[ 4] ; mCHt2[2]=ffC[ 5] ;
980 222 : mCHt0[3]=ffC[ 6] ; mCHt1[3]=ffC[ 7] ; mCHt2[3]=ffC[ 8] ;
981 222 : mCHt0[4]=ffC[10] ; mCHt1[4]=ffC[11] ; mCHt2[4]=ffC[12] ;
982 222 : mCHt0[5]=ffC[15] ; mCHt1[5]=ffC[16] ; mCHt2[5]=ffC[17] ;
983 222 : mCHt0[6]=ffC[21] ; mCHt1[6]=ffC[22] ; mCHt2[6]=ffC[23] ;
984 :
985 : //* Kalman gain K = mCH'*S
986 :
987 222 : Double_t k0[7], k1[7], k2[7];
988 :
989 3552 : for(Int_t i=0;i<7;++i){
990 1554 : k0[i] = mCHt0[i]*mS[0] + mCHt1[i]*mS[1] + mCHt2[i]*mS[3];
991 1554 : k1[i] = mCHt0[i]*mS[1] + mCHt1[i]*mS[2] + mCHt2[i]*mS[4];
992 1554 : k2[i] = mCHt0[i]*mS[3] + mCHt1[i]*mS[4] + mCHt2[i]*mS[5];
993 : }
994 :
995 : //* New estimation of the vertex position
996 :
997 222 : if( iter<maxIter-1 ){
998 1184 : for(Int_t i=0; i<3; ++i)
999 444 : fVtxGuess[i]= ffP[i] + k0[i]*zeta[0]+k1[i]*zeta[1]+k2[i]*zeta[2];
1000 148 : continue;
1001 : }
1002 :
1003 : // last itearation -> update the particle
1004 :
1005 : //* VHt = VH'
1006 :
1007 74 : Double_t mVHt0[7], mVHt1[7], mVHt2[7];
1008 :
1009 74 : mVHt0[0]=mV[ 0] ; mVHt1[0]=mV[ 1] ; mVHt2[0]=mV[ 3] ;
1010 74 : mVHt0[1]=mV[ 1] ; mVHt1[1]=mV[ 2] ; mVHt2[1]=mV[ 4] ;
1011 74 : mVHt0[2]=mV[ 3] ; mVHt1[2]=mV[ 4] ; mVHt2[2]=mV[ 5] ;
1012 74 : mVHt0[3]=mV[ 6] ; mVHt1[3]=mV[ 7] ; mVHt2[3]=mV[ 8] ;
1013 74 : mVHt0[4]=mV[10] ; mVHt1[4]=mV[11] ; mVHt2[4]=mV[12] ;
1014 74 : mVHt0[5]=mV[15] ; mVHt1[5]=mV[16] ; mVHt2[5]=mV[17] ;
1015 74 : mVHt0[6]=mV[21] ; mVHt1[6]=mV[22] ; mVHt2[6]=mV[23] ;
1016 :
1017 : //* Kalman gain Km = mCH'*S
1018 :
1019 74 : Double_t km0[7], km1[7], km2[7];
1020 :
1021 1184 : for(Int_t i=0;i<7;++i){
1022 518 : km0[i] = mVHt0[i]*mS[0] + mVHt1[i]*mS[1] + mVHt2[i]*mS[3];
1023 518 : km1[i] = mVHt0[i]*mS[1] + mVHt1[i]*mS[2] + mVHt2[i]*mS[4];
1024 518 : km2[i] = mVHt0[i]*mS[3] + mVHt1[i]*mS[4] + mVHt2[i]*mS[5];
1025 : }
1026 :
1027 1184 : for(Int_t i=0;i<7;++i)
1028 518 : ffP[i] = ffP[i] + k0[i]*zeta[0] + k1[i]*zeta[1] + k2[i]*zeta[2];
1029 :
1030 1184 : for(Int_t i=0;i<7;++i)
1031 518 : m[i] = m[i] - km0[i]*zeta[0] - km1[i]*zeta[1] - km2[i]*zeta[2];
1032 :
1033 1184 : for(Int_t i=0, k=0;i<7;++i){
1034 5180 : for(Int_t j=0;j<=i;++j,++k){
1035 2072 : ffC[k] = ffC[k] - (k0[i]*mCHt0[j] + k1[i]*mCHt1[j] + k2[i]*mCHt2[j] );
1036 : }
1037 : }
1038 :
1039 1184 : for(Int_t i=0, k=0;i<7;++i){
1040 5180 : for(Int_t j=0;j<=i;++j,++k){
1041 2072 : mV[k] = mV[k] - (km0[i]*mVHt0[j] + km1[i]*mVHt1[j] + km2[i]*mVHt2[j] );
1042 : }
1043 : }
1044 :
1045 74 : Double_t mDf[7][7];
1046 :
1047 1184 : for(Int_t i=0;i<7;++i){
1048 8288 : for(Int_t j=0;j<7;++j){
1049 3626 : mDf[i][j] = (km0[i]*mCHt0[j] + km1[i]*mCHt1[j] + km2[i]*mCHt2[j] );
1050 : }
1051 : }
1052 :
1053 74 : Double_t mJ1[7][7], mJ2[7][7];
1054 1184 : for(Int_t iPar1=0; iPar1<7; iPar1++)
1055 : {
1056 8288 : for(Int_t iPar2=0; iPar2<7; iPar2++)
1057 : {
1058 3626 : mJ1[iPar1][iPar2] = 0;
1059 3626 : mJ2[iPar1][iPar2] = 0;
1060 : }
1061 : }
1062 :
1063 74 : Double_t mMassParticle = ffP[6]*ffP[6] - (ffP[3]*ffP[3] + ffP[4]*ffP[4] + ffP[5]*ffP[5]);
1064 74 : Double_t mMassDaughter = m[6]*m[6] - (m[3]*m[3] + m[4]*m[4] + m[5]*m[5]);
1065 136 : if(mMassParticle > 0) mMassParticle = TMath::Sqrt(mMassParticle);
1066 130 : if(mMassDaughter > 0) mMassDaughter = TMath::Sqrt(mMassDaughter);
1067 :
1068 74 : if( fMassHypo > -0.5)
1069 74 : SetMassConstraint(ffP,ffC,mJ1,fMassHypo);
1070 0 : else if((mMassParticle < SumDaughterMass) || (ffP[6]<0) )
1071 0 : SetMassConstraint(ffP,ffC,mJ1,SumDaughterMass);
1072 :
1073 74 : if(Daughter.fMassHypo > -0.5)
1074 74 : SetMassConstraint(m,mV,mJ2,Daughter.fMassHypo);
1075 0 : else if((mMassDaughter < Daughter.SumDaughterMass) || (m[6] < 0) )
1076 0 : SetMassConstraint(m,mV,mJ2,Daughter.SumDaughterMass);
1077 :
1078 74 : Double_t mDJ[7][7];
1079 :
1080 1184 : for(Int_t i=0; i<7; i++) {
1081 8288 : for(Int_t j=0; j<7; j++) {
1082 3626 : mDJ[i][j] = 0;
1083 58016 : for(Int_t k=0; k<7; k++) {
1084 25382 : mDJ[i][j] += mDf[i][k]*mJ1[j][k];
1085 : }
1086 : }
1087 : }
1088 :
1089 1184 : for(Int_t i=0; i<7; ++i){
1090 8288 : for(Int_t j=0; j<7; ++j){
1091 3626 : mDf[i][j]=0;
1092 58016 : for(Int_t l=0; l<7; l++){
1093 25382 : mDf[i][j] += mJ2[i][l]*mDJ[l][j];
1094 : }
1095 : }
1096 : }
1097 :
1098 : //* Add the daughter momentum to the particle momentum
1099 :
1100 74 : ffP[ 3] += m[ 3];
1101 74 : ffP[ 4] += m[ 4];
1102 74 : ffP[ 5] += m[ 5];
1103 74 : ffP[ 6] += m[ 6];
1104 :
1105 74 : ffC[ 9] += mV[ 9];
1106 74 : ffC[13] += mV[13];
1107 74 : ffC[14] += mV[14];
1108 74 : ffC[18] += mV[18];
1109 74 : ffC[19] += mV[19];
1110 74 : ffC[20] += mV[20];
1111 74 : ffC[24] += mV[24];
1112 74 : ffC[25] += mV[25];
1113 74 : ffC[26] += mV[26];
1114 74 : ffC[27] += mV[27];
1115 :
1116 74 : ffC[6 ] += mDf[3][0]; ffC[7 ] += mDf[3][1]; ffC[8 ] += mDf[3][2];
1117 74 : ffC[10] += mDf[4][0]; ffC[11] += mDf[4][1]; ffC[12] += mDf[4][2];
1118 74 : ffC[15] += mDf[5][0]; ffC[16] += mDf[5][1]; ffC[17] += mDf[5][2];
1119 74 : ffC[21] += mDf[6][0]; ffC[22] += mDf[6][1]; ffC[23] += mDf[6][2];
1120 :
1121 74 : ffC[9 ] += mDf[3][3] + mDf[3][3];
1122 74 : ffC[13] += mDf[4][3] + mDf[3][4]; ffC[14] += mDf[4][4] + mDf[4][4];
1123 74 : ffC[18] += mDf[5][3] + mDf[3][5]; ffC[19] += mDf[5][4] + mDf[4][5]; ffC[20] += mDf[5][5] + mDf[5][5];
1124 74 : ffC[24] += mDf[6][3] + mDf[3][6]; ffC[25] += mDf[6][4] + mDf[4][6]; ffC[26] += mDf[6][5] + mDf[5][6]; ffC[27] += mDf[6][6] + mDf[6][6];
1125 :
1126 : //* New estimation of the vertex position r += K*zeta
1127 :
1128 1184 : for(Int_t i=0;i<7;++i)
1129 518 : fP[i] = ffP[i];
1130 :
1131 : //* New covariance matrix C -= K*(mCH')'
1132 :
1133 1184 : for(Int_t i=0, k=0;i<7;++i){
1134 5180 : for(Int_t j=0;j<=i;++j,++k){
1135 2072 : fC[k] = ffC[k];
1136 : }
1137 : }
1138 : //* Calculate Chi^2
1139 :
1140 74 : fNDF += 2;
1141 74 : fQ += Daughter.GetQ();
1142 74 : fSFromDecay = 0;
1143 148 : fChi2 += (mS[0]*zeta[0] + mS[1]*zeta[1] + mS[3]*zeta[2])*zeta[0]
1144 74 : + (mS[1]*zeta[0] + mS[2]*zeta[1] + mS[4]*zeta[2])*zeta[1]
1145 74 : + (mS[3]*zeta[0] + mS[4]*zeta[1] + mS[5]*zeta[2])*zeta[2];
1146 296 : }
1147 74 : }
1148 :
1149 : void AliKFParticleBase::SetProductionVertex( const AliKFParticleBase &Vtx )
1150 : {
1151 : //* Set production vertex for the particle, when the particle was not used in the vertex fit
1152 :
1153 148 : const Double_t *m = Vtx.fP, *mV = Vtx.fC;
1154 :
1155 74 : Bool_t noS = ( fC[35]<=0 ); // no decay length allowed
1156 :
1157 74 : if( noS ){
1158 0 : TransportToDecayVertex();
1159 0 : fP[7] = 0;
1160 0 : fC[28] = fC[29] = fC[30] = fC[31] = fC[32] = fC[33] = fC[34] = fC[35] = 0;
1161 0 : } else {
1162 74 : TransportToDS( GetDStoPoint( m ) );
1163 74 : fP[7] = -fSFromDecay;
1164 74 : fC[28] = fC[29] = fC[30] = fC[31] = fC[32] = fC[33] = fC[34] = 0;
1165 74 : fC[35] = 0.1;
1166 :
1167 74 : Convert(1);
1168 : }
1169 :
1170 74 : Double_t mAi[6];
1171 :
1172 74 : InvertSym3( fC, mAi );
1173 :
1174 : Double_t mB[5][3];
1175 :
1176 74 : mB[0][0] = fC[ 6]*mAi[0] + fC[ 7]*mAi[1] + fC[ 8]*mAi[3];
1177 74 : mB[0][1] = fC[ 6]*mAi[1] + fC[ 7]*mAi[2] + fC[ 8]*mAi[4];
1178 74 : mB[0][2] = fC[ 6]*mAi[3] + fC[ 7]*mAi[4] + fC[ 8]*mAi[5];
1179 :
1180 74 : mB[1][0] = fC[10]*mAi[0] + fC[11]*mAi[1] + fC[12]*mAi[3];
1181 74 : mB[1][1] = fC[10]*mAi[1] + fC[11]*mAi[2] + fC[12]*mAi[4];
1182 74 : mB[1][2] = fC[10]*mAi[3] + fC[11]*mAi[4] + fC[12]*mAi[5];
1183 :
1184 74 : mB[2][0] = fC[15]*mAi[0] + fC[16]*mAi[1] + fC[17]*mAi[3];
1185 74 : mB[2][1] = fC[15]*mAi[1] + fC[16]*mAi[2] + fC[17]*mAi[4];
1186 74 : mB[2][2] = fC[15]*mAi[3] + fC[16]*mAi[4] + fC[17]*mAi[5];
1187 :
1188 74 : mB[3][0] = fC[21]*mAi[0] + fC[22]*mAi[1] + fC[23]*mAi[3];
1189 74 : mB[3][1] = fC[21]*mAi[1] + fC[22]*mAi[2] + fC[23]*mAi[4];
1190 74 : mB[3][2] = fC[21]*mAi[3] + fC[22]*mAi[4] + fC[23]*mAi[5];
1191 :
1192 74 : mB[4][0] = fC[28]*mAi[0] + fC[29]*mAi[1] + fC[30]*mAi[3];
1193 74 : mB[4][1] = fC[28]*mAi[1] + fC[29]*mAi[2] + fC[30]*mAi[4];
1194 74 : mB[4][2] = fC[28]*mAi[3] + fC[29]*mAi[4] + fC[30]*mAi[5];
1195 :
1196 74 : Double_t z[3] = { m[0]-fP[0], m[1]-fP[1], m[2]-fP[2] };
1197 :
1198 : {
1199 296 : Double_t mAVi[6] = { fC[0]-mV[0], fC[1]-mV[1], fC[2]-mV[2],
1200 222 : fC[3]-mV[3], fC[4]-mV[4], fC[5]-mV[5] };
1201 :
1202 74 : if( !InvertSym3( mAVi, mAVi ) ){
1203 :
1204 74 : Double_t dChi2 = ( +(mAVi[0]*z[0] + mAVi[1]*z[1] + mAVi[3]*z[2])*z[0]
1205 74 : +(mAVi[1]*z[0] + mAVi[2]*z[1] + mAVi[4]*z[2])*z[1]
1206 74 : +(mAVi[3]*z[0] + mAVi[4]*z[1] + mAVi[5]*z[2])*z[2] );
1207 :
1208 : // Take Abs(dChi2) here. Negative value of 'det' or 'dChi2' shows that the particle
1209 : // was not used in the production vertex fit
1210 :
1211 74 : fChi2+= TMath::Abs( dChi2 );
1212 74 : }
1213 74 : fNDF += 2;
1214 74 : }
1215 :
1216 74 : fP[0] = m[0];
1217 74 : fP[1] = m[1];
1218 74 : fP[2] = m[2];
1219 74 : fP[3]+= mB[0][0]*z[0] + mB[0][1]*z[1] + mB[0][2]*z[2];
1220 74 : fP[4]+= mB[1][0]*z[0] + mB[1][1]*z[1] + mB[1][2]*z[2];
1221 74 : fP[5]+= mB[2][0]*z[0] + mB[2][1]*z[1] + mB[2][2]*z[2];
1222 74 : fP[6]+= mB[3][0]*z[0] + mB[3][1]*z[1] + mB[3][2]*z[2];
1223 74 : fP[7]+= mB[4][0]*z[0] + mB[4][1]*z[1] + mB[4][2]*z[2];
1224 :
1225 : Double_t d0, d1, d2;
1226 :
1227 74 : fC[0] = mV[0];
1228 74 : fC[1] = mV[1];
1229 74 : fC[2] = mV[2];
1230 74 : fC[3] = mV[3];
1231 74 : fC[4] = mV[4];
1232 74 : fC[5] = mV[5];
1233 :
1234 74 : d0= mB[0][0]*mV[0] + mB[0][1]*mV[1] + mB[0][2]*mV[3] - fC[ 6];
1235 74 : d1= mB[0][0]*mV[1] + mB[0][1]*mV[2] + mB[0][2]*mV[4] - fC[ 7];
1236 74 : d2= mB[0][0]*mV[3] + mB[0][1]*mV[4] + mB[0][2]*mV[5] - fC[ 8];
1237 :
1238 74 : fC[ 6]+= d0;
1239 74 : fC[ 7]+= d1;
1240 74 : fC[ 8]+= d2;
1241 74 : fC[ 9]+= d0*mB[0][0] + d1*mB[0][1] + d2*mB[0][2];
1242 :
1243 74 : d0= mB[1][0]*mV[0] + mB[1][1]*mV[1] + mB[1][2]*mV[3] - fC[10];
1244 74 : d1= mB[1][0]*mV[1] + mB[1][1]*mV[2] + mB[1][2]*mV[4] - fC[11];
1245 74 : d2= mB[1][0]*mV[3] + mB[1][1]*mV[4] + mB[1][2]*mV[5] - fC[12];
1246 :
1247 74 : fC[10]+= d0;
1248 74 : fC[11]+= d1;
1249 74 : fC[12]+= d2;
1250 74 : fC[13]+= d0*mB[0][0] + d1*mB[0][1] + d2*mB[0][2];
1251 74 : fC[14]+= d0*mB[1][0] + d1*mB[1][1] + d2*mB[1][2];
1252 :
1253 74 : d0= mB[2][0]*mV[0] + mB[2][1]*mV[1] + mB[2][2]*mV[3] - fC[15];
1254 74 : d1= mB[2][0]*mV[1] + mB[2][1]*mV[2] + mB[2][2]*mV[4] - fC[16];
1255 74 : d2= mB[2][0]*mV[3] + mB[2][1]*mV[4] + mB[2][2]*mV[5] - fC[17];
1256 :
1257 74 : fC[15]+= d0;
1258 74 : fC[16]+= d1;
1259 74 : fC[17]+= d2;
1260 74 : fC[18]+= d0*mB[0][0] + d1*mB[0][1] + d2*mB[0][2];
1261 74 : fC[19]+= d0*mB[1][0] + d1*mB[1][1] + d2*mB[1][2];
1262 74 : fC[20]+= d0*mB[2][0] + d1*mB[2][1] + d2*mB[2][2];
1263 :
1264 74 : d0= mB[3][0]*mV[0] + mB[3][1]*mV[1] + mB[3][2]*mV[3] - fC[21];
1265 74 : d1= mB[3][0]*mV[1] + mB[3][1]*mV[2] + mB[3][2]*mV[4] - fC[22];
1266 74 : d2= mB[3][0]*mV[3] + mB[3][1]*mV[4] + mB[3][2]*mV[5] - fC[23];
1267 :
1268 74 : fC[21]+= d0;
1269 74 : fC[22]+= d1;
1270 74 : fC[23]+= d2;
1271 74 : fC[24]+= d0*mB[0][0] + d1*mB[0][1] + d2*mB[0][2];
1272 74 : fC[25]+= d0*mB[1][0] + d1*mB[1][1] + d2*mB[1][2];
1273 74 : fC[26]+= d0*mB[2][0] + d1*mB[2][1] + d2*mB[2][2];
1274 74 : fC[27]+= d0*mB[3][0] + d1*mB[3][1] + d2*mB[3][2];
1275 :
1276 74 : d0= mB[4][0]*mV[0] + mB[4][1]*mV[1] + mB[4][2]*mV[3] - fC[28];
1277 74 : d1= mB[4][0]*mV[1] + mB[4][1]*mV[2] + mB[4][2]*mV[4] - fC[29];
1278 74 : d2= mB[4][0]*mV[3] + mB[4][1]*mV[4] + mB[4][2]*mV[5] - fC[30];
1279 :
1280 74 : fC[28]+= d0;
1281 74 : fC[29]+= d1;
1282 74 : fC[30]+= d2;
1283 74 : fC[31]+= d0*mB[0][0] + d1*mB[0][1] + d2*mB[0][2];
1284 74 : fC[32]+= d0*mB[1][0] + d1*mB[1][1] + d2*mB[1][2];
1285 74 : fC[33]+= d0*mB[2][0] + d1*mB[2][1] + d2*mB[2][2];
1286 74 : fC[34]+= d0*mB[3][0] + d1*mB[3][1] + d2*mB[3][2];
1287 74 : fC[35]+= d0*mB[4][0] + d1*mB[4][1] + d2*mB[4][2];
1288 :
1289 74 : if( noS ){
1290 0 : fP[7] = 0;
1291 0 : fC[28] = fC[29] = fC[30] = fC[31] = fC[32] = fC[33] = fC[34] = fC[35] = 0;
1292 0 : } else {
1293 74 : TransportToDS( fP[7] );
1294 74 : Convert(0);
1295 : }
1296 :
1297 74 : fSFromDecay = 0;
1298 74 : }
1299 :
1300 : void AliKFParticleBase::SetMassConstraint( Double_t *mP, Double_t *mC, Double_t mJ[7][7], Double_t mass )
1301 : {
1302 : //* Set nonlinear mass constraint (Mass) on the state vector mP with a covariance matrix mC.
1303 :
1304 296 : const Double_t energy2 = mP[6]*mP[6], p2 = mP[3]*mP[3]+mP[4]*mP[4]+mP[5]*mP[5], mass2 = mass*mass;
1305 :
1306 148 : const Double_t a = energy2 - p2 + 2.*mass2;
1307 148 : const Double_t b = -2.*(energy2 + p2);
1308 148 : const Double_t c = energy2 - p2 - mass2;
1309 :
1310 : Double_t lambda = 0;
1311 296 : if(TMath::Abs(b) > 1.e-10) lambda = -c / b;
1312 :
1313 148 : Double_t d = 4.*energy2*p2 - mass2*(energy2-p2-2.*mass2);
1314 444 : if(d>=0 && TMath::Abs(a) > 1.e-10) lambda = (energy2 + p2 - sqrt(d))/a;
1315 :
1316 148 : if(mP[6] < 0) //If energy < 0 we need a lambda < 0
1317 0 : lambda = -1000000.; //Empirical, a better solution should be found
1318 :
1319 : Int_t iIter=0;
1320 306 : for(iIter=0; iIter<100; iIter++)
1321 : {
1322 153 : Double_t lambda2 = lambda*lambda;
1323 153 : Double_t lambda4 = lambda2*lambda2;
1324 :
1325 : Double_t lambda0 = lambda;
1326 :
1327 153 : Double_t f = -mass2 * lambda4 + a*lambda2 + b*lambda + c;
1328 153 : Double_t df = -4.*mass2 * lambda2*lambda + 2.*a*lambda + b;
1329 153 : if(TMath::Abs(df) < 1.e-10) break;
1330 153 : lambda -= f/df;
1331 301 : if(TMath::Abs(lambda0 - lambda) < 1.e-8) break;
1332 5 : }
1333 :
1334 148 : const Double_t lpi = 1./(1. + lambda);
1335 148 : const Double_t lmi = 1./(1. - lambda);
1336 148 : const Double_t lp2i = lpi*lpi;
1337 148 : const Double_t lm2i = lmi*lmi;
1338 :
1339 148 : Double_t lambda2 = lambda*lambda;
1340 :
1341 148 : Double_t dfl = -4.*mass2 * lambda2*lambda + 2.*a*lambda + b;
1342 148 : Double_t dfx[7] = {0};//,0,0,0};
1343 148 : dfx[0] = -2.*(1. + lambda)*(1. + lambda)*mP[3];
1344 148 : dfx[1] = -2.*(1. + lambda)*(1. + lambda)*mP[4];
1345 148 : dfx[2] = -2.*(1. + lambda)*(1. + lambda)*mP[5];
1346 148 : dfx[3] = 2.*(1. - lambda)*(1. - lambda)*mP[6];
1347 148 : Double_t dlx[4] = {1,1,1,1};
1348 148 : if(TMath::Abs(dfl) > 1.e-10 )
1349 : {
1350 1480 : for(int i=0; i<4; i++)
1351 592 : dlx[i] = -dfx[i] / dfl;
1352 148 : }
1353 :
1354 148 : Double_t dxx[4] = {mP[3]*lm2i, mP[4]*lm2i, mP[5]*lm2i, -mP[6]*lp2i};
1355 :
1356 2368 : for(Int_t i=0; i<7; i++)
1357 16576 : for(Int_t j=0; j<7; j++)
1358 7252 : mJ[i][j]=0;
1359 148 : mJ[0][0] = 1.;
1360 148 : mJ[1][1] = 1.;
1361 148 : mJ[2][2] = 1.;
1362 :
1363 1480 : for(Int_t i=3; i<7; i++)
1364 5920 : for(Int_t j=3; j<7; j++)
1365 2368 : mJ[i][j] = dlx[j-3]*dxx[i-3];
1366 :
1367 1184 : for(Int_t i=3; i<6; i++)
1368 444 : mJ[i][i] += lmi;
1369 148 : mJ[6][6] += lpi;
1370 :
1371 148 : Double_t mCJ[7][7];
1372 :
1373 2368 : for(Int_t i=0; i<7; i++) {
1374 16576 : for(Int_t j=0; j<7; j++) {
1375 7252 : mCJ[i][j] = 0;
1376 116032 : for(Int_t k=0; k<7; k++) {
1377 50764 : mCJ[i][j] += mC[IJ(i,k)]*mJ[j][k];
1378 : }
1379 : }
1380 : }
1381 :
1382 2368 : for(Int_t i=0; i<7; ++i){
1383 10360 : for(Int_t j=0; j<=i; ++j){
1384 4144 : mC[IJ(i,j)]=0;
1385 66304 : for(Int_t l=0; l<7; l++){
1386 29008 : mC[IJ(i,j)] += mJ[i][l]*mCJ[l][j];
1387 : }
1388 : }
1389 : }
1390 :
1391 148 : mP[3] *= lmi;
1392 148 : mP[4] *= lmi;
1393 148 : mP[5] *= lmi;
1394 148 : mP[6] *= lpi;
1395 148 : }
1396 :
1397 : void AliKFParticleBase::SetNonlinearMassConstraint( Double_t mass )
1398 : {
1399 : //* Set nonlinear mass constraint (mass)
1400 :
1401 0 : Double_t mJ[7][7];
1402 0 : SetMassConstraint( fP, fC, mJ, mass );
1403 0 : fMassHypo = mass;
1404 0 : SumDaughterMass = mass;
1405 0 : }
1406 :
1407 : void AliKFParticleBase::SetMassConstraint( Double_t Mass, Double_t SigmaMass )
1408 : {
1409 : //* Set hard( SigmaMass=0 ) or soft (SigmaMass>0) mass constraint
1410 :
1411 36 : fMassHypo = Mass;
1412 18 : SumDaughterMass = Mass;
1413 :
1414 18 : Double_t m2 = Mass*Mass; // measurement, weighted by Mass
1415 18 : Double_t s2 = m2*SigmaMass*SigmaMass; // sigma^2
1416 :
1417 18 : Double_t p2 = fP[3]*fP[3] + fP[4]*fP[4] + fP[5]*fP[5];
1418 18 : Double_t e0 = TMath::Sqrt(m2+p2);
1419 :
1420 18 : Double_t mH[8];
1421 18 : mH[0] = mH[1] = mH[2] = 0.;
1422 18 : mH[3] = -2*fP[3];
1423 18 : mH[4] = -2*fP[4];
1424 18 : mH[5] = -2*fP[5];
1425 18 : mH[6] = 2*fP[6];//e0;
1426 18 : mH[7] = 0;
1427 :
1428 18 : Double_t zeta = e0*e0 - e0*fP[6];
1429 18 : zeta = m2 - (fP[6]*fP[6]-p2);
1430 :
1431 18 : Double_t mCHt[8], s2_est=0;
1432 324 : for( Int_t i=0; i<8; ++i ){
1433 144 : mCHt[i] = 0.0;
1434 2592 : for (Int_t j=0;j<8;++j) mCHt[i] += Cij(i,j)*mH[j];
1435 144 : s2_est += mH[i]*mCHt[i];
1436 : }
1437 :
1438 18 : if( s2_est<1.e-20 ) return; // calculated mass error is already 0,
1439 : // the particle can not be constrained on mass
1440 :
1441 18 : Double_t w2 = 1./( s2 + s2_est );
1442 18 : fChi2 += zeta*zeta*w2;
1443 18 : fNDF += 1;
1444 324 : for( Int_t i=0, ii=0; i<8; ++i ){
1445 144 : Double_t ki = mCHt[i]*w2;
1446 144 : fP[i]+= ki*zeta;
1447 1584 : for(Int_t j=0;j<=i;++j) fC[ii++] -= ki*mCHt[j];
1448 : }
1449 36 : }
1450 :
1451 :
1452 : void AliKFParticleBase::SetNoDecayLength()
1453 : {
1454 : //* Set no decay length for resonances
1455 :
1456 0 : TransportToDecayVertex();
1457 :
1458 0 : Double_t h[8];
1459 0 : h[0] = h[1] = h[2] = h[3] = h[4] = h[5] = h[6] = 0;
1460 0 : h[7] = 1;
1461 :
1462 0 : Double_t zeta = 0 - fP[7];
1463 0 : for(Int_t i=0;i<8;++i) zeta -= h[i]*(fP[i]-fP[i]);
1464 :
1465 0 : Double_t s = fC[35];
1466 0 : if( s>1.e-20 ){
1467 0 : s = 1./s;
1468 0 : fChi2 += zeta*zeta*s;
1469 0 : fNDF += 1;
1470 0 : for( Int_t i=0, ii=0; i<7; ++i ){
1471 0 : Double_t ki = fC[28+i]*s;
1472 0 : fP[i]+= ki*zeta;
1473 0 : for(Int_t j=0;j<=i;++j) fC[ii++] -= ki*fC[28+j];
1474 : }
1475 0 : }
1476 0 : fP[7] = 0;
1477 0 : fC[28] = fC[29] = fC[30] = fC[31] = fC[32] = fC[33] = fC[34] = fC[35] = 0;
1478 0 : }
1479 :
1480 :
1481 : void AliKFParticleBase::Construct( const AliKFParticleBase* vDaughters[], Int_t NDaughters,
1482 : const AliKFParticleBase *Parent, Double_t Mass, Bool_t IsConstrained )
1483 : {
1484 : //* Full reconstruction in one go
1485 :
1486 : Int_t maxIter = 1;
1487 0 : bool wasLinearized = fIsLinearized;
1488 0 : if( !fIsLinearized || IsConstrained ){
1489 : //fVtxGuess[0] = fVtxGuess[1] = fVtxGuess[2] = 0; //!!!!
1490 0 : fVtxGuess[0] = GetX();
1491 0 : fVtxGuess[1] = GetY();
1492 0 : fVtxGuess[2] = GetZ();
1493 0 : fIsLinearized = 1;
1494 : maxIter = 3;
1495 0 : }
1496 :
1497 0 : Double_t constraintC[6];
1498 :
1499 0 : if( IsConstrained ){
1500 0 : for(Int_t i=0;i<6;++i) constraintC[i]=fC[i];
1501 0 : } else {
1502 0 : for(Int_t i=0;i<6;++i) constraintC[i]=0.;
1503 0 : constraintC[0] = constraintC[2] = constraintC[5] = 100.;
1504 : }
1505 :
1506 :
1507 0 : for( Int_t iter=0; iter<maxIter; iter++ ){
1508 0 : fAtProductionVertex = 0;
1509 0 : fSFromDecay = 0;
1510 0 : fP[0] = fVtxGuess[0];
1511 0 : fP[1] = fVtxGuess[1];
1512 0 : fP[2] = fVtxGuess[2];
1513 0 : fP[3] = 0;
1514 0 : fP[4] = 0;
1515 0 : fP[5] = 0;
1516 0 : fP[6] = 0;
1517 0 : fP[7] = 0;
1518 0 : SumDaughterMass = 0;
1519 :
1520 0 : for(Int_t i=0;i<6; ++i) fC[i]=constraintC[i];
1521 0 : for(Int_t i=6;i<36;++i) fC[i]=0.;
1522 0 : fC[35] = 1.;
1523 :
1524 0 : fNDF = IsConstrained ?0 :-3;
1525 0 : fChi2 = 0.;
1526 0 : fQ = 0;
1527 :
1528 0 : for( Int_t itr =0; itr<NDaughters; itr++ ){
1529 0 : AddDaughter( *vDaughters[itr] );
1530 : }
1531 0 : if( iter<maxIter-1){
1532 0 : for( Int_t i=0; i<3; i++ ) fVtxGuess[i] = fP[i];
1533 0 : }
1534 : }
1535 0 : fIsLinearized = wasLinearized;
1536 :
1537 0 : if( Mass>=0 ) SetMassConstraint( Mass );
1538 0 : if( Parent ) SetProductionVertex( *Parent );
1539 0 : }
1540 :
1541 :
1542 : void AliKFParticleBase::Convert( bool ToProduction )
1543 : {
1544 : //* Tricky function - convert the particle error along its trajectory to
1545 : //* the value which corresponds to its production/decay vertex
1546 : //* It is done by combination of the error of decay length with the position errors
1547 :
1548 148 : Double_t fld[3];
1549 : {
1550 148 : GetFieldValue( fP, fld );
1551 148 : const Double_t kCLight = fQ*0.000299792458;
1552 148 : fld[0]*=kCLight; fld[1]*=kCLight; fld[2]*=kCLight;
1553 : }
1554 :
1555 : Double_t h[6];
1556 :
1557 148 : h[0] = fP[3];
1558 148 : h[1] = fP[4];
1559 148 : h[2] = fP[5];
1560 222 : if( ToProduction ){ h[0]=-h[0]; h[1]=-h[1]; h[2]=-h[2]; }
1561 148 : h[3] = h[1]*fld[2]-h[2]*fld[1];
1562 148 : h[4] = h[2]*fld[0]-h[0]*fld[2];
1563 148 : h[5] = h[0]*fld[1]-h[1]*fld[0];
1564 :
1565 : Double_t c;
1566 :
1567 148 : c = fC[28]+h[0]*fC[35];
1568 148 : fC[ 0]+= h[0]*(c+fC[28]);
1569 148 : fC[28] = c;
1570 :
1571 148 : fC[ 1]+= h[1]*fC[28] + h[0]*fC[29];
1572 148 : c = fC[29]+h[1]*fC[35];
1573 148 : fC[ 2]+= h[1]*(c+fC[29]);
1574 148 : fC[29] = c;
1575 :
1576 148 : fC[ 3]+= h[2]*fC[28] + h[0]*fC[30];
1577 148 : fC[ 4]+= h[2]*fC[29] + h[1]*fC[30];
1578 148 : c = fC[30]+h[2]*fC[35];
1579 148 : fC[ 5]+= h[2]*(c+fC[30]);
1580 148 : fC[30] = c;
1581 :
1582 148 : fC[ 6]+= h[3]*fC[28] + h[0]*fC[31];
1583 148 : fC[ 7]+= h[3]*fC[29] + h[1]*fC[31];
1584 148 : fC[ 8]+= h[3]*fC[30] + h[2]*fC[31];
1585 148 : c = fC[31]+h[3]*fC[35];
1586 148 : fC[ 9]+= h[3]*(c+fC[31]);
1587 148 : fC[31] = c;
1588 :
1589 148 : fC[10]+= h[4]*fC[28] + h[0]*fC[32];
1590 148 : fC[11]+= h[4]*fC[29] + h[1]*fC[32];
1591 148 : fC[12]+= h[4]*fC[30] + h[2]*fC[32];
1592 148 : fC[13]+= h[4]*fC[31] + h[3]*fC[32];
1593 148 : c = fC[32]+h[4]*fC[35];
1594 148 : fC[14]+= h[4]*(c+fC[32]);
1595 148 : fC[32] = c;
1596 :
1597 148 : fC[15]+= h[5]*fC[28] + h[0]*fC[33];
1598 148 : fC[16]+= h[5]*fC[29] + h[1]*fC[33];
1599 148 : fC[17]+= h[5]*fC[30] + h[2]*fC[33];
1600 148 : fC[18]+= h[5]*fC[31] + h[3]*fC[33];
1601 148 : fC[19]+= h[5]*fC[32] + h[4]*fC[33];
1602 148 : c = fC[33]+h[5]*fC[35];
1603 148 : fC[20]+= h[5]*(c+fC[33]);
1604 148 : fC[33] = c;
1605 :
1606 148 : fC[21]+= h[0]*fC[34];
1607 148 : fC[22]+= h[1]*fC[34];
1608 148 : fC[23]+= h[2]*fC[34];
1609 148 : fC[24]+= h[3]*fC[34];
1610 148 : fC[25]+= h[4]*fC[34];
1611 148 : fC[26]+= h[5]*fC[34];
1612 148 : }
1613 :
1614 :
1615 : void AliKFParticleBase::TransportToDecayVertex()
1616 : {
1617 : //* Transport the particle to its decay vertex
1618 :
1619 148 : if( fSFromDecay != 0 ) TransportToDS( -fSFromDecay );
1620 74 : if( fAtProductionVertex ) Convert(0);
1621 74 : fAtProductionVertex = 0;
1622 74 : }
1623 :
1624 : void AliKFParticleBase::TransportToProductionVertex()
1625 : {
1626 : //* Transport the particle to its production vertex
1627 :
1628 0 : if( fSFromDecay != -fP[7] ) TransportToDS( -fSFromDecay-fP[7] );
1629 0 : if( !fAtProductionVertex ) Convert( 1 );
1630 0 : fAtProductionVertex = 1;
1631 0 : }
1632 :
1633 :
1634 : void AliKFParticleBase::TransportToDS( Double_t dS )
1635 : {
1636 : //* Transport the particle on dS parameter (SignedPath/Momentum)
1637 :
1638 444 : Transport( dS, fP, fC );
1639 222 : fSFromDecay+= dS;
1640 222 : }
1641 :
1642 :
1643 : Double_t AliKFParticleBase::GetDStoPointLine( const Double_t xyz[] ) const
1644 : {
1645 : //* Get dS to a certain space point without field
1646 :
1647 0 : Double_t p2 = fP[3]*fP[3] + fP[4]*fP[4] + fP[5]*fP[5];
1648 0 : if( p2<1.e-4 ) p2 = 1;
1649 0 : return ( fP[3]*(xyz[0]-fP[0]) + fP[4]*(xyz[1]-fP[1]) + fP[5]*(xyz[2]-fP[2]) )/p2;
1650 : }
1651 :
1652 :
1653 : Double_t AliKFParticleBase::GetDStoPointBz( Double_t B, const Double_t xyz[] )
1654 : const
1655 : {
1656 :
1657 : //* Get dS to a certain space point for Bz field
1658 : const Double_t kCLight = 0.000299792458;
1659 1036 : Double_t bq = B*fQ*kCLight;
1660 518 : Double_t pt2 = fP[3]*fP[3] + fP[4]*fP[4];
1661 518 : if( pt2<1.e-4 ) return 0;
1662 518 : Double_t dx = xyz[0] - fP[0];
1663 518 : Double_t dy = xyz[1] - fP[1];
1664 518 : Double_t a = dx*fP[3]+dy*fP[4];
1665 : Double_t dS;
1666 :
1667 564 : if( TMath::Abs(bq)<1.e-8 ) dS = a/pt2;
1668 472 : else dS = TMath::ATan2( bq*a, pt2 + bq*(dy*fP[3] -dx*fP[4]) )/bq;
1669 :
1670 : if(0){
1671 :
1672 : Double_t px = fP[3];
1673 : Double_t py = fP[4];
1674 : Double_t pz = fP[5];
1675 : Double_t ss[2], g[2][5];
1676 :
1677 : ss[0] = dS;
1678 : ss[1] = -dS;
1679 : for( Int_t i=0; i<2; i++){
1680 : Double_t bs = bq*ss[i];
1681 : Double_t c = TMath::Cos(bs), s = TMath::Sin(bs);
1682 : Double_t cB,sB;
1683 : if( TMath::Abs(bq)>1.e-8){
1684 : cB= (1-c)/bq;
1685 : sB= s/bq;
1686 : }else{
1687 : const Double_t kOvSqr6 = 1./TMath::Sqrt(6.);
1688 : sB = (1.-bs*kOvSqr6)*(1.+bs*kOvSqr6)*ss[i];
1689 : cB = .5*sB*bs;
1690 : }
1691 : g[i][0] = fP[0] + sB*px + cB*py;
1692 : g[i][1] = fP[1] - cB*px + sB*py;
1693 : g[i][2] = fP[2] + ss[i]*pz;
1694 : g[i][3] = + c*px + s*py;
1695 : g[i][4] = - s*px + c*py;
1696 : }
1697 :
1698 : Int_t i=0;
1699 :
1700 : Double_t dMin = 1.e10;
1701 : for( Int_t j=0; j<2; j++){
1702 : Double_t xx = g[j][0]-xyz[0];
1703 : Double_t yy = g[j][1]-xyz[1];
1704 : Double_t zz = g[j][2]-xyz[2];
1705 : Double_t d = xx*xx + yy*yy + zz*zz;
1706 : if( d<dMin ){
1707 : dMin = d;
1708 : i = j;
1709 : }
1710 : }
1711 :
1712 : dS = ss[i];
1713 :
1714 : Double_t x= g[i][0], y= g[i][1], z= g[i][2], ppx= g[i][3], ppy= g[i][4];
1715 : Double_t ddx = x-xyz[0];
1716 : Double_t ddy = y-xyz[1];
1717 : Double_t ddz = z-xyz[2];
1718 : Double_t c = ddx*ppx + ddy*ppy + ddz*pz ;
1719 : Double_t pp2 = ppx*ppx + ppy*ppy + pz*pz;
1720 : if( TMath::Abs(pp2)>1.e-8 ){
1721 : dS+=c/pp2;
1722 : }
1723 : }
1724 : return dS;
1725 518 : }
1726 :
1727 :
1728 : void AliKFParticleBase::GetDStoParticleBz( Double_t B, const AliKFParticleBase &p,
1729 : Double_t &DS, Double_t &DS1 )
1730 : const
1731 : {
1732 : //* Get dS to another particle for Bz field
1733 148 : Double_t px = fP[3];
1734 74 : Double_t py = fP[4];
1735 74 : Double_t pz = fP[5];
1736 :
1737 74 : Double_t px1 = p.fP[3];
1738 74 : Double_t py1 = p.fP[4];
1739 74 : Double_t pz1 = p.fP[5];
1740 :
1741 : const Double_t kCLight = 0.000299792458;
1742 :
1743 74 : Double_t bq = B*fQ*kCLight;
1744 74 : Double_t bq1 = B*p.fQ*kCLight;
1745 : Double_t s=0, ds=0, s1=0, ds1=0;
1746 :
1747 74 : if( TMath::Abs(bq)>1.e-8 || TMath::Abs(bq1)>1.e-8 ){
1748 :
1749 74 : Double_t dx = (p.fP[0] - fP[0]);
1750 74 : Double_t dy = (p.fP[1] - fP[1]);
1751 74 : Double_t d2 = (dx*dx+dy*dy);
1752 :
1753 74 : Double_t p2 = (px *px + py *py);
1754 74 : Double_t p21 = (px1*px1 + py1*py1);
1755 :
1756 148 : if( TMath::Abs(p2) < 1.e-8 || TMath::Abs(p21) < 1.e-8 )
1757 : {
1758 0 : DS=0.;
1759 0 : DS1=0.;
1760 0 : return;
1761 : }
1762 :
1763 74 : Double_t a = (px*py1 - py*px1);
1764 74 : Double_t b = (px*px1 + py*py1);
1765 :
1766 74 : Double_t ldx = bq*bq1*dx - bq1*py + bq*py1 ;
1767 74 : Double_t ldy = bq*bq1*dy + bq1*px - bq*px1 ;
1768 74 : Double_t l2 = ldx*ldx + ldy*ldy;
1769 :
1770 74 : Double_t cS = bq1*p2 + bq*bq1*(dy* px - dx* py) - bq*b;
1771 74 : Double_t cS1= bq*p21 - bq*bq1*(dy*px1 - dx*py1) - bq1*b;
1772 :
1773 74 : Double_t ca = bq*bq*bq1*d2 +2*( cS + bq*bq*(py1*dx-px1*dy)) ;
1774 74 : Double_t ca1 = bq*bq1*bq1*d2 +2*( cS1 - bq1*bq1*(py*dx-px*dy)) ;
1775 :
1776 74 : Double_t sa = 4*l2*p2 - ca*ca;
1777 74 : Double_t sa1 = 4*l2*p21 - ca1*ca1;
1778 :
1779 82 : if(sa<0) sa=0;
1780 82 : if(sa1<0)sa1=0;
1781 :
1782 74 : if( TMath::Abs(bq)>1.e-8){
1783 74 : s = TMath::ATan2( bq*( bq1*(dx*px +dy*py) + a ) , cS )/bq;
1784 74 : ds = TMath::ATan2(TMath::Sqrt(sa),ca)/bq;
1785 74 : } else {
1786 0 : s = ( (dx*px + dy*py) + (py*px1-px*py1)/bq1)/p2;
1787 0 : ds = s*s - (d2-2*(px1*dy-py1*dx)/bq1)/p2;
1788 0 : if( ds<0 ) ds = 0;
1789 0 : ds = TMath::Sqrt(ds);
1790 : }
1791 :
1792 74 : if( TMath::Abs(bq1)>1.e-8){
1793 74 : s1 = TMath::ATan2( -bq1*( bq*(dx*px1+dy*py1) + a), cS1 )/bq1;
1794 74 : ds1 = TMath::ATan2(TMath::Sqrt(sa1),ca1)/bq1;
1795 74 : } else {
1796 0 : s1 = (-(dx*px1 + dy*py1) + (py*px1-px*py1)/bq)/p21;
1797 0 : ds1 = s1*s1 - (d2+2*(px*dy-py*dx)/bq)/p21;
1798 0 : if( ds1<0 ) ds1 = 0;
1799 0 : ds1 = TMath::Sqrt(ds1);
1800 : }
1801 74 : }
1802 :
1803 74 : Double_t ss[2], ss1[2], g[2][5],g1[2][5];
1804 :
1805 74 : ss[0] = s + ds;
1806 74 : ss[1] = s - ds;
1807 74 : ss1[0] = s1 + ds1;
1808 74 : ss1[1] = s1 - ds1;
1809 444 : for( Int_t i=0; i<2; i++){
1810 148 : Double_t bs = bq*ss[i];
1811 148 : Double_t c = TMath::Cos(bs), sss = TMath::Sin(bs);
1812 : Double_t cB,sB;
1813 148 : if( TMath::Abs(bq)>1.e-8){
1814 148 : cB= (1-c)/bq;
1815 148 : sB= sss/bq;
1816 148 : }else{
1817 0 : const Double_t kOvSqr6 = 1./TMath::Sqrt(6.);
1818 0 : sB = (1.-bs*kOvSqr6)*(1.+bs*kOvSqr6)*ss[i];
1819 0 : cB = .5*sB*bs;
1820 : }
1821 148 : g[i][0] = fP[0] + sB*px + cB*py;
1822 148 : g[i][1] = fP[1] - cB*px + sB*py;
1823 148 : g[i][2] = fP[2] + ss[i]*pz;
1824 148 : g[i][3] = + c*px + sss*py;
1825 148 : g[i][4] = - sss*px + c*py;
1826 :
1827 148 : bs = bq1*ss1[i];
1828 148 : c = TMath::Cos(bs); sss = TMath::Sin(bs);
1829 148 : if( TMath::Abs(bq1)>1.e-8){
1830 148 : cB= (1-c)/bq1;
1831 148 : sB= sss/bq1;
1832 148 : }else{
1833 0 : const Double_t kOvSqr6 = 1./TMath::Sqrt(6.);
1834 0 : sB = (1.-bs*kOvSqr6)*(1.+bs*kOvSqr6)*ss1[i];
1835 0 : cB = .5*sB*bs;
1836 : }
1837 :
1838 148 : g1[i][0] = p.fP[0] + sB*px1 + cB*py1;
1839 148 : g1[i][1] = p.fP[1] - cB*px1 + sB*py1;
1840 148 : g1[i][2] = p.fP[2] + ss[i]*pz1;
1841 148 : g1[i][3] = + c*px1 + sss*py1;
1842 148 : g1[i][4] = - sss*px1 + c*py1;
1843 : }
1844 :
1845 : Int_t i=0, i1=0;
1846 :
1847 : Double_t dMin = 1.e10;
1848 444 : for( Int_t j=0; j<2; j++){
1849 888 : for( Int_t j1=0; j1<2; j1++){
1850 296 : Double_t xx = g[j][0]-g1[j1][0];
1851 296 : Double_t yy = g[j][1]-g1[j1][1];
1852 296 : Double_t zz = g[j][2]-g1[j1][2];
1853 296 : Double_t d = xx*xx + yy*yy + zz*zz;
1854 296 : if( d<dMin ){
1855 : dMin = d;
1856 : i = j;
1857 : i1 = j1;
1858 90 : }
1859 : }
1860 : }
1861 :
1862 74 : DS = ss[i];
1863 74 : DS1 = ss1[i1];
1864 : if(0){
1865 : Double_t x= g[i][0], y= g[i][1], z= g[i][2], ppx= g[i][3], ppy= g[i][4];
1866 : Double_t x1=g1[i1][0], y1= g1[i1][1], z1= g1[i1][2], ppx1= g1[i1][3], ppy1= g1[i1][4];
1867 : Double_t dx = x1-x;
1868 : Double_t dy = y1-y;
1869 : Double_t dz = z1-z;
1870 : Double_t a = ppx*ppx1 + ppy*ppy1 + pz*pz1;
1871 : Double_t b = dx*ppx1 + dy*ppy1 + dz*pz1;
1872 : Double_t c = dx*ppx + dy*ppy + dz*pz ;
1873 : Double_t pp2 = ppx*ppx + ppy*ppy + pz*pz;
1874 : Double_t pp21= ppx1*ppx1 + ppy1*ppy1 + pz1*pz1;
1875 : Double_t det = pp2*pp21 - a*a;
1876 : if( TMath::Abs(det)>1.e-8 ){
1877 : DS+=(a*b-pp21*c)/det;
1878 : DS1+=(a*c-pp2*b)/det;
1879 : }
1880 : }
1881 148 : }
1882 :
1883 :
1884 :
1885 : void AliKFParticleBase::TransportCBM( Double_t dS,
1886 : Double_t P[], Double_t C[] ) const
1887 : {
1888 : //* Transport the particle on dS, output to P[],C[], for CBM field
1889 :
1890 0 : if( fQ==0 ){
1891 0 : TransportLine( dS, P, C );
1892 0 : return;
1893 : }
1894 :
1895 : const Double_t kCLight = 0.000299792458;
1896 :
1897 0 : Double_t c = fQ*kCLight;
1898 :
1899 : // construct coefficients
1900 :
1901 : Double_t
1902 0 : px = fP[3],
1903 0 : py = fP[4],
1904 0 : pz = fP[5];
1905 :
1906 : Double_t sx=0, sy=0, sz=0, syy=0, syz=0, syyy=0, ssx=0, ssy=0, ssz=0, ssyy=0, ssyz=0, ssyyy=0;
1907 :
1908 : { // get field integrals
1909 :
1910 0 : Double_t fld[3][3];
1911 0 : Double_t p0[3], p1[3], p2[3];
1912 :
1913 : // line track approximation
1914 :
1915 0 : p0[0] = fP[0];
1916 0 : p0[1] = fP[1];
1917 0 : p0[2] = fP[2];
1918 :
1919 0 : p2[0] = fP[0] + px*dS;
1920 0 : p2[1] = fP[1] + py*dS;
1921 0 : p2[2] = fP[2] + pz*dS;
1922 :
1923 0 : p1[0] = 0.5*(p0[0]+p2[0]);
1924 0 : p1[1] = 0.5*(p0[1]+p2[1]);
1925 0 : p1[2] = 0.5*(p0[2]+p2[2]);
1926 :
1927 : // first order track approximation
1928 : {
1929 0 : GetFieldValue( p0, fld[0] );
1930 0 : GetFieldValue( p1, fld[1] );
1931 0 : GetFieldValue( p2, fld[2] );
1932 :
1933 0 : Double_t ssy1 = ( 7*fld[0][1] + 6*fld[1][1]-fld[2][1] )*c*dS*dS/96.;
1934 0 : Double_t ssy2 = ( fld[0][1] + 2*fld[1][1] )*c*dS*dS/6.;
1935 :
1936 0 : p1[0] -= ssy1*pz;
1937 0 : p1[2] += ssy1*px;
1938 0 : p2[0] -= ssy2*pz;
1939 0 : p2[2] += ssy2*px;
1940 : }
1941 :
1942 0 : GetFieldValue( p0, fld[0] );
1943 0 : GetFieldValue( p1, fld[1] );
1944 0 : GetFieldValue( p2, fld[2] );
1945 :
1946 0 : sx = c*( fld[0][0] + 4*fld[1][0] + fld[2][0] )*dS/6.;
1947 0 : sy = c*( fld[0][1] + 4*fld[1][1] + fld[2][1] )*dS/6.;
1948 0 : sz = c*( fld[0][2] + 4*fld[1][2] + fld[2][2] )*dS/6.;
1949 :
1950 0 : ssx = c*( fld[0][0] + 2*fld[1][0])*dS*dS/6.;
1951 0 : ssy = c*( fld[0][1] + 2*fld[1][1])*dS*dS/6.;
1952 0 : ssz = c*( fld[0][2] + 2*fld[1][2])*dS*dS/6.;
1953 :
1954 0 : Double_t c2[3][3] = { { 5, -4, -1},{ 44, 80, -4},{ 11, 44, 5} }; // /=360.
1955 0 : Double_t cc2[3][3] = { { 38, 8, -4},{ 148, 208, -20},{ 3, 36, 3} }; // /=2520.
1956 0 : for(Int_t n=0; n<3; n++)
1957 0 : for(Int_t m=0; m<3; m++)
1958 : {
1959 0 : syz += c2[n][m]*fld[n][1]*fld[m][2];
1960 0 : ssyz += cc2[n][m]*fld[n][1]*fld[m][2];
1961 : }
1962 :
1963 0 : syz *= c*c*dS*dS/360.;
1964 0 : ssyz *= c*c*dS*dS*dS/2520.;
1965 :
1966 0 : syy = c*( fld[0][1] + 4*fld[1][1] + fld[2][1] )*dS;
1967 0 : syyy = syy*syy*syy / 1296;
1968 0 : syy = syy*syy/72;
1969 :
1970 0 : ssyy = ( fld[0][1]*( 38*fld[0][1] + 156*fld[1][1] - fld[2][1] )+
1971 0 : fld[1][1]*( 208*fld[1][1] +16*fld[2][1] )+
1972 0 : fld[2][1]*( 3*fld[2][1] )
1973 0 : )*dS*dS*dS*c*c/2520.;
1974 : ssyyy =
1975 : (
1976 0 : fld[0][1]*( fld[0][1]*( 85*fld[0][1] + 526*fld[1][1] - 7*fld[2][1] )+
1977 0 : fld[1][1]*( 1376*fld[1][1] +84*fld[2][1] )+
1978 0 : fld[2][1]*( 19*fld[2][1] ) )+
1979 0 : fld[1][1]*( fld[1][1]*( 1376*fld[1][1] +256*fld[2][1] )+
1980 0 : fld[2][1]*( 62*fld[2][1] ) )+
1981 0 : fld[2][1]*fld[2][1] *( 3*fld[2][1] )
1982 0 : )*dS*dS*dS*dS*c*c*c/90720.;
1983 :
1984 0 : }
1985 :
1986 0 : Double_t mJ[8][8];
1987 0 : for( Int_t i=0; i<8; i++ ) for( Int_t j=0; j<8; j++) mJ[i][j]=0;
1988 :
1989 0 : mJ[0][0]=1; mJ[0][1]=0; mJ[0][2]=0; mJ[0][3]=dS-ssyy; mJ[0][4]=ssx; mJ[0][5]=ssyyy-ssy;
1990 0 : mJ[1][0]=0; mJ[1][1]=1; mJ[1][2]=0; mJ[1][3]=-ssz; mJ[1][4]=dS; mJ[1][5]=ssx+ssyz;
1991 0 : mJ[2][0]=0; mJ[2][1]=0; mJ[2][2]=1; mJ[2][3]=ssy-ssyyy; mJ[2][4]=-ssx; mJ[2][5]=dS-ssyy;
1992 :
1993 0 : mJ[3][0]=0; mJ[3][1]=0; mJ[3][2]=0; mJ[3][3]=1-syy; mJ[3][4]=sx; mJ[3][5]=syyy-sy;
1994 0 : mJ[4][0]=0; mJ[4][1]=0; mJ[4][2]=0; mJ[4][3]=-sz; mJ[4][4]=1; mJ[4][5]=sx+syz;
1995 0 : mJ[5][0]=0; mJ[5][1]=0; mJ[5][2]=0; mJ[5][3]=sy-syyy; mJ[5][4]=-sx; mJ[5][5]=1-syy;
1996 0 : mJ[6][6] = mJ[7][7] = 1;
1997 :
1998 0 : P[0] = fP[0] + mJ[0][3]*px + mJ[0][4]*py + mJ[0][5]*pz;
1999 0 : P[1] = fP[1] + mJ[1][3]*px + mJ[1][4]*py + mJ[1][5]*pz;
2000 0 : P[2] = fP[2] + mJ[2][3]*px + mJ[2][4]*py + mJ[2][5]*pz;
2001 0 : P[3] = mJ[3][3]*px + mJ[3][4]*py + mJ[3][5]*pz;
2002 0 : P[4] = mJ[4][3]*px + mJ[4][4]*py + mJ[4][5]*pz;
2003 0 : P[5] = mJ[5][3]*px + mJ[5][4]*py + mJ[5][5]*pz;
2004 0 : P[6] = fP[6];
2005 0 : P[7] = fP[7];
2006 :
2007 0 : MultQSQt( mJ[0], fC, C);
2008 :
2009 0 : }
2010 :
2011 :
2012 : void AliKFParticleBase::TransportBz( Double_t b, Double_t t,
2013 : Double_t p[], Double_t e[] ) const
2014 : {
2015 : //* Transport the particle on dS, output to P[],C[], for Bz field
2016 :
2017 : const Double_t kCLight = 0.000299792458;
2018 1480 : b = b*fQ*kCLight;
2019 740 : Double_t bs= b*t;
2020 740 : Double_t s = TMath::Sin(bs), c = TMath::Cos(bs);
2021 : Double_t sB, cB;
2022 740 : if( TMath::Abs(bs)>1.e-10){
2023 378 : sB= s/b;
2024 378 : cB= (1-c)/b;
2025 378 : }else{
2026 362 : const Double_t kOvSqr6 = 1./TMath::Sqrt(6.);
2027 362 : sB = (1.-bs*kOvSqr6)*(1.+bs*kOvSqr6)*t;
2028 362 : cB = .5*sB*bs;
2029 : }
2030 :
2031 740 : Double_t px = fP[3];
2032 740 : Double_t py = fP[4];
2033 740 : Double_t pz = fP[5];
2034 :
2035 740 : p[0] = fP[0] + sB*px + cB*py;
2036 740 : p[1] = fP[1] - cB*px + sB*py;
2037 740 : p[2] = fP[2] + t*pz;
2038 740 : p[3] = c*px + s*py;
2039 740 : p[4] = -s*px + c*py;
2040 740 : p[5] = fP[5];
2041 740 : p[6] = fP[6];
2042 740 : p[7] = fP[7];
2043 :
2044 : /*
2045 : Double_t mJ[8][8] = { {1,0,0, sB, cB, 0, 0, 0 },
2046 : {0,1,0, -cB, sB, 0, 0, 0 },
2047 : {0,0,1, 0, 0, t, 0, 0 },
2048 : {0,0,0, c, s, 0, 0, 0 },
2049 : {0,0,0, -s, c, 0, 0, 0 },
2050 : {0,0,0, 0, 0, 1, 0, 0 },
2051 : {0,0,0, 0, 0, 0, 1, 0 },
2052 : {0,0,0, 0, 0, 0, 0, 1 } };
2053 : Double_t mA[8][8];
2054 : for( Int_t k=0,i=0; i<8; i++)
2055 : for( Int_t j=0; j<=i; j++, k++ ) mA[i][j] = mA[j][i] = fC[k];
2056 :
2057 : Double_t mJC[8][8];
2058 : for( Int_t i=0; i<8; i++ )
2059 : for( Int_t j=0; j<8; j++ ){
2060 : mJC[i][j]=0;
2061 : for( Int_t k=0; k<8; k++ ) mJC[i][j]+=mJ[i][k]*mA[k][j];
2062 : }
2063 :
2064 : for( Int_t k=0,i=0; i<8; i++)
2065 : for( Int_t j=0; j<=i; j++, k++ ){
2066 : e[k] = 0;
2067 : for( Int_t l=0; l<8; l++ ) e[k]+=mJC[i][l]*mJ[j][l];
2068 : }
2069 :
2070 : return;
2071 : */
2072 :
2073 : Double_t
2074 740 : c6=fC[6], c7=fC[7], c8=fC[8], c17=fC[17], c18=fC[18],
2075 740 : c24 = fC[24], c31 = fC[31];
2076 :
2077 : Double_t
2078 740 : cBC13 = cB*fC[13],
2079 740 : mJC13 = c7 - cB*fC[9] + sB*fC[13],
2080 740 : mJC14 = fC[11] - cBC13 + sB*fC[14],
2081 740 : mJC23 = c8 + t*c18,
2082 740 : mJC24 = fC[12] + t*fC[19],
2083 740 : mJC33 = c*fC[9] + s*fC[13],
2084 740 : mJC34 = c*fC[13] + s*fC[14],
2085 740 : mJC43 = -s*fC[9] + c*fC[13],
2086 740 : mJC44 = -s*fC[13] + c*fC[14];
2087 :
2088 :
2089 740 : e[0]= fC[0] + 2*(sB*c6 + cB*fC[10]) + (sB*fC[9] + 2*cBC13)*sB + cB*cB*fC[14];
2090 740 : e[1]= fC[1] - cB*c6 + sB*fC[10] + mJC13*sB + mJC14*cB;
2091 740 : e[2]= fC[2] - cB*c7 + sB*fC[11] - mJC13*cB + mJC14*sB;
2092 740 : e[3]= fC[3] + t*fC[15] + mJC23*sB + mJC24*cB;
2093 740 : e[4]= fC[4] + t*fC[16] - mJC23*cB + mJC24*sB;
2094 :
2095 740 : e[15]= fC[15] + c18*sB + fC[19]*cB;
2096 740 : e[16]= fC[16] - c18*cB + fC[19]*sB;
2097 740 : e[17]= c17 + fC[20]*t;
2098 740 : e[18]= c18*c + fC[19]*s;
2099 740 : e[19]= -c18*s + fC[19]*c;
2100 :
2101 740 : e[5]= fC[5] + (c17 + e[17] )*t;
2102 :
2103 740 : e[6]= c*c6 + s*fC[10] + mJC33*sB + mJC34*cB;
2104 740 : e[7]= c*c7 + s*fC[11] - mJC33*cB + mJC34*sB;
2105 740 : e[8]= c*c8 + s*fC[12] + e[18]*t;
2106 740 : e[9]= mJC33*c + mJC34*s;
2107 740 : e[10]= -s*c6 + c*fC[10] + mJC43*sB + mJC44*cB;
2108 :
2109 :
2110 740 : e[11]= -s*c7 + c*fC[11] - mJC43*cB + mJC44*sB;
2111 740 : e[12]= -s*c8 + c*fC[12] + e[19]*t;
2112 740 : e[13]= mJC43*c + mJC44*s;
2113 740 : e[14]= -mJC43*s + mJC44*c;
2114 740 : e[20]= fC[20];
2115 740 : e[21]= fC[21] + fC[25]*cB + c24*sB;
2116 740 : e[22]= fC[22] - c24*cB + fC[25]*sB;
2117 740 : e[23]= fC[23] + fC[26]*t;
2118 740 : e[24]= c*c24 + s*fC[25];
2119 740 : e[25]= c*fC[25] - c24*s;
2120 740 : e[26]= fC[26];
2121 740 : e[27]= fC[27];
2122 740 : e[28]= fC[28] + fC[32]*cB + c31*sB;
2123 740 : e[29]= fC[29] - c31*cB + fC[32]*sB;
2124 740 : e[30]= fC[30] + fC[33]*t;
2125 740 : e[31]= c*c31 + s*fC[32];
2126 740 : e[32]= c*fC[32] - s*c31;
2127 740 : e[33]= fC[33];
2128 740 : e[34]= fC[34];
2129 740 : e[35]= fC[35];
2130 740 : }
2131 :
2132 :
2133 : Double_t AliKFParticleBase::GetDistanceFromVertex( const AliKFParticleBase &Vtx ) const
2134 : {
2135 : //* Calculate distance from vertex [cm]
2136 :
2137 0 : return GetDistanceFromVertex( Vtx.fP );
2138 : }
2139 :
2140 : Double_t AliKFParticleBase::GetDistanceFromVertex( const Double_t vtx[] ) const
2141 : {
2142 : //* Calculate distance from vertex [cm]
2143 :
2144 0 : Double_t mP[8], mC[36];
2145 0 : Transport( GetDStoPoint(vtx), mP, mC );
2146 0 : Double_t d[3]={ vtx[0]-mP[0], vtx[1]-mP[1], vtx[2]-mP[2]};
2147 0 : return TMath::Sqrt( d[0]*d[0]+d[1]*d[1]+d[2]*d[2] );
2148 0 : }
2149 :
2150 : Double_t AliKFParticleBase::GetDistanceFromParticle( const AliKFParticleBase &p )
2151 : const
2152 : {
2153 : //* Calculate distance to other particle [cm]
2154 :
2155 0 : Double_t dS, dS1;
2156 0 : GetDStoParticle( p, dS, dS1 );
2157 0 : Double_t mP[8], mC[36], mP1[8], mC1[36];
2158 0 : Transport( dS, mP, mC );
2159 0 : p.Transport( dS1, mP1, mC1 );
2160 0 : Double_t dx = mP[0]-mP1[0];
2161 0 : Double_t dy = mP[1]-mP1[1];
2162 0 : Double_t dz = mP[2]-mP1[2];
2163 : dz = 0;
2164 0 : return TMath::Sqrt(dx*dx+dy*dy+dz*dz);
2165 0 : }
2166 :
2167 : Double_t AliKFParticleBase::GetDeviationFromVertex( const AliKFParticleBase &Vtx ) const
2168 : {
2169 : //* Calculate sqrt(Chi2/ndf) deviation from vertex
2170 :
2171 0 : return GetDeviationFromVertex( Vtx.fP, Vtx.fC );
2172 : }
2173 :
2174 :
2175 : Double_t AliKFParticleBase::GetDeviationFromVertex( const Double_t v[], const Double_t Cv[] ) const
2176 : {
2177 : //* Calculate sqrt(Chi2/ndf) deviation from vertex
2178 : //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
2179 :
2180 0 : Double_t mP[8];
2181 0 : Double_t mC[36];
2182 :
2183 0 : Transport( GetDStoPoint(v), mP, mC );
2184 :
2185 0 : Double_t d[3]={ v[0]-mP[0], v[1]-mP[1], v[2]-mP[2]};
2186 :
2187 0 : Double_t sigmaS = .1+10.*TMath::Sqrt( (d[0]*d[0]+d[1]*d[1]+d[2]*d[2])/
2188 0 : (mP[3]*mP[3]+mP[4]*mP[4]+mP[5]*mP[5]) );
2189 :
2190 :
2191 0 : Double_t h[3] = { mP[3]*sigmaS, mP[4]*sigmaS, mP[5]*sigmaS };
2192 :
2193 : Double_t mSi[6] =
2194 0 : { mC[0] +h[0]*h[0],
2195 0 : mC[1] +h[1]*h[0], mC[2] +h[1]*h[1],
2196 0 : mC[3] +h[2]*h[0], mC[4] +h[2]*h[1], mC[5] +h[2]*h[2] };
2197 :
2198 0 : if( Cv ){
2199 0 : mSi[0]+=Cv[0];
2200 0 : mSi[1]+=Cv[1];
2201 0 : mSi[2]+=Cv[2];
2202 0 : mSi[3]+=Cv[3];
2203 0 : mSi[4]+=Cv[4];
2204 0 : mSi[5]+=Cv[5];
2205 0 : }
2206 :
2207 : Double_t mS[6];
2208 :
2209 0 : mS[0] = mSi[2]*mSi[5] - mSi[4]*mSi[4];
2210 0 : mS[1] = mSi[3]*mSi[4] - mSi[1]*mSi[5];
2211 0 : mS[2] = mSi[0]*mSi[5] - mSi[3]*mSi[3];
2212 0 : mS[3] = mSi[1]*mSi[4] - mSi[2]*mSi[3];
2213 0 : mS[4] = mSi[1]*mSi[3] - mSi[0]*mSi[4];
2214 0 : mS[5] = mSi[0]*mSi[2] - mSi[1]*mSi[1];
2215 :
2216 0 : Double_t s = ( mSi[0]*mS[0] + mSi[1]*mS[1] + mSi[3]*mS[3] );
2217 0 : s = ( s > 1.E-20 ) ?1./s :0;
2218 :
2219 0 : return TMath::Sqrt( TMath::Abs(s*( ( mS[0]*d[0] + mS[1]*d[1] + mS[3]*d[2])*d[0]
2220 0 : +(mS[1]*d[0] + mS[2]*d[1] + mS[4]*d[2])*d[1]
2221 0 : +(mS[3]*d[0] + mS[4]*d[1] + mS[5]*d[2])*d[2] ))/2);
2222 0 : }
2223 :
2224 :
2225 : Double_t AliKFParticleBase::GetDeviationFromParticle( const AliKFParticleBase &p )
2226 : const
2227 : {
2228 : //* Calculate sqrt(Chi2/ndf) deviation from other particle
2229 :
2230 0 : Double_t dS, dS1;
2231 0 : GetDStoParticle( p, dS, dS1 );
2232 0 : Double_t mP1[8], mC1[36];
2233 0 : p.Transport( dS1, mP1, mC1 );
2234 :
2235 0 : Double_t d[3]={ fP[0]-mP1[0], fP[1]-mP1[1], fP[2]-mP1[2]};
2236 :
2237 0 : Double_t sigmaS = .1+10.*TMath::Sqrt( (d[0]*d[0]+d[1]*d[1]+d[2]*d[2])/
2238 0 : (mP1[3]*mP1[3]+mP1[4]*mP1[4]+mP1[5]*mP1[5]) );
2239 :
2240 0 : Double_t h[3] = { mP1[3]*sigmaS, mP1[4]*sigmaS, mP1[5]*sigmaS };
2241 :
2242 0 : mC1[0] +=h[0]*h[0];
2243 0 : mC1[1] +=h[1]*h[0];
2244 0 : mC1[2] +=h[1]*h[1];
2245 0 : mC1[3] +=h[2]*h[0];
2246 0 : mC1[4] +=h[2]*h[1];
2247 0 : mC1[5] +=h[2]*h[2];
2248 :
2249 0 : return GetDeviationFromVertex( mP1, mC1 )*TMath::Sqrt(2./1.);
2250 0 : }
2251 :
2252 :
2253 :
2254 : void AliKFParticleBase::SubtractFromVertex( AliKFParticleBase &Vtx ) const
2255 : {
2256 : //* Subtract the particle from the vertex
2257 :
2258 0 : Double_t fld[3];
2259 : {
2260 0 : GetFieldValue( Vtx.fP, fld );
2261 : const Double_t kCLight = 0.000299792458;
2262 0 : fld[0]*=kCLight; fld[1]*=kCLight; fld[2]*=kCLight;
2263 : }
2264 :
2265 0 : Double_t m[8];
2266 0 : Double_t mCm[36];
2267 :
2268 0 : if( Vtx.fIsLinearized ){
2269 0 : GetMeasurement( Vtx.fVtxGuess, m, mCm );
2270 0 : } else {
2271 0 : GetMeasurement( Vtx.fP, m, mCm );
2272 : }
2273 :
2274 : Double_t mV[6];
2275 :
2276 0 : mV[ 0] = mCm[ 0];
2277 0 : mV[ 1] = mCm[ 1];
2278 0 : mV[ 2] = mCm[ 2];
2279 0 : mV[ 3] = mCm[ 3];
2280 0 : mV[ 4] = mCm[ 4];
2281 0 : mV[ 5] = mCm[ 5];
2282 :
2283 : //*
2284 :
2285 : Double_t mS[6];
2286 : {
2287 0 : Double_t mSi[6] = { mV[0]-Vtx.fC[0],
2288 0 : mV[1]-Vtx.fC[1], mV[2]-Vtx.fC[2],
2289 0 : mV[3]-Vtx.fC[3], mV[4]-Vtx.fC[4], mV[5]-Vtx.fC[5] };
2290 :
2291 0 : mS[0] = mSi[2]*mSi[5] - mSi[4]*mSi[4];
2292 0 : mS[1] = mSi[3]*mSi[4] - mSi[1]*mSi[5];
2293 0 : mS[2] = mSi[0]*mSi[5] - mSi[3]*mSi[3];
2294 0 : mS[3] = mSi[1]*mSi[4] - mSi[2]*mSi[3];
2295 0 : mS[4] = mSi[1]*mSi[3] - mSi[0]*mSi[4];
2296 0 : mS[5] = mSi[0]*mSi[2] - mSi[1]*mSi[1];
2297 :
2298 0 : Double_t s = ( mSi[0]*mS[0] + mSi[1]*mS[1] + mSi[3]*mS[3] );
2299 0 : s = ( s > 1.E-20 ) ?1./s :0;
2300 0 : mS[0]*=s;
2301 0 : mS[1]*=s;
2302 0 : mS[2]*=s;
2303 0 : mS[3]*=s;
2304 0 : mS[4]*=s;
2305 0 : mS[5]*=s;
2306 : }
2307 :
2308 : //* Residual (measured - estimated)
2309 :
2310 0 : Double_t zeta[3] = { m[0]-Vtx.fP[0], m[1]-Vtx.fP[1], m[2]-Vtx.fP[2] };
2311 :
2312 : //* mCHt = mCH' - D'
2313 :
2314 0 : Double_t mCHt0[3], mCHt1[3], mCHt2[3];
2315 :
2316 0 : mCHt0[0]=Vtx.fC[ 0] ; mCHt1[0]=Vtx.fC[ 1] ; mCHt2[0]=Vtx.fC[ 3] ;
2317 0 : mCHt0[1]=Vtx.fC[ 1] ; mCHt1[1]=Vtx.fC[ 2] ; mCHt2[1]=Vtx.fC[ 4] ;
2318 0 : mCHt0[2]=Vtx.fC[ 3] ; mCHt1[2]=Vtx.fC[ 4] ; mCHt2[2]=Vtx.fC[ 5] ;
2319 :
2320 : //* Kalman gain K = mCH'*S
2321 :
2322 0 : Double_t k0[3], k1[3], k2[3];
2323 :
2324 0 : for(Int_t i=0;i<3;++i){
2325 0 : k0[i] = mCHt0[i]*mS[0] + mCHt1[i]*mS[1] + mCHt2[i]*mS[3];
2326 0 : k1[i] = mCHt0[i]*mS[1] + mCHt1[i]*mS[2] + mCHt2[i]*mS[4];
2327 0 : k2[i] = mCHt0[i]*mS[3] + mCHt1[i]*mS[4] + mCHt2[i]*mS[5];
2328 : }
2329 :
2330 : //* New estimation of the vertex position r += K*zeta
2331 :
2332 0 : Double_t dChi2 = -(mS[0]*zeta[0] + mS[1]*zeta[1] + mS[3]*zeta[2])*zeta[0]
2333 0 : + (mS[1]*zeta[0] + mS[2]*zeta[1] + mS[4]*zeta[2])*zeta[1]
2334 0 : + (mS[3]*zeta[0] + mS[4]*zeta[1] + mS[5]*zeta[2])*zeta[2];
2335 :
2336 0 : if( Vtx.fChi2 - dChi2 < 0 ) return;
2337 :
2338 0 : for(Int_t i=0;i<3;++i)
2339 0 : Vtx.fP[i] -= k0[i]*zeta[0] + k1[i]*zeta[1] + k2[i]*zeta[2];
2340 :
2341 : //* New covariance matrix C -= K*(mCH')'
2342 :
2343 0 : for(Int_t i=0, k=0;i<3;++i){
2344 0 : for(Int_t j=0;j<=i;++j,++k)
2345 0 : Vtx.fC[k] += k0[i]*mCHt0[j] + k1[i]*mCHt1[j] + k2[i]*mCHt2[j];
2346 : }
2347 :
2348 : //* Calculate Chi^2
2349 :
2350 0 : Vtx.fNDF -= 2;
2351 0 : Vtx.fChi2 -= dChi2;
2352 0 : }
2353 :
2354 :
2355 :
2356 : void AliKFParticleBase::TransportLine( Double_t dS,
2357 : Double_t P[], Double_t C[] ) const
2358 : {
2359 : //* Transport the particle as a straight line
2360 :
2361 0 : P[0] = fP[0] + dS*fP[3];
2362 0 : P[1] = fP[1] + dS*fP[4];
2363 0 : P[2] = fP[2] + dS*fP[5];
2364 0 : P[3] = fP[3];
2365 0 : P[4] = fP[4];
2366 0 : P[5] = fP[5];
2367 0 : P[6] = fP[6];
2368 0 : P[7] = fP[7];
2369 :
2370 0 : Double_t c6 = fC[ 6] + dS*fC[ 9];
2371 0 : Double_t c11 = fC[11] + dS*fC[14];
2372 0 : Double_t c17 = fC[17] + dS*fC[20];
2373 0 : Double_t sc13 = dS*fC[13];
2374 0 : Double_t sc18 = dS*fC[18];
2375 0 : Double_t sc19 = dS*fC[19];
2376 :
2377 0 : C[ 0] = fC[ 0] + dS*( fC[ 6] + c6 );
2378 0 : C[ 2] = fC[ 2] + dS*( fC[11] + c11 );
2379 0 : C[ 5] = fC[ 5] + dS*( fC[17] + c17 );
2380 :
2381 0 : C[ 7] = fC[ 7] + sc13;
2382 0 : C[ 8] = fC[ 8] + sc18;
2383 0 : C[ 9] = fC[ 9];
2384 :
2385 0 : C[12] = fC[12] + sc19;
2386 :
2387 0 : C[ 1] = fC[ 1] + dS*( fC[10] + C[ 7] );
2388 0 : C[ 3] = fC[ 3] + dS*( fC[15] + C[ 8] );
2389 0 : C[ 4] = fC[ 4] + dS*( fC[16] + C[12] );
2390 0 : C[ 6] = c6;
2391 :
2392 0 : C[10] = fC[10] + sc13;
2393 0 : C[11] = c11;
2394 :
2395 0 : C[13] = fC[13];
2396 0 : C[14] = fC[14];
2397 0 : C[15] = fC[15] + sc18;
2398 0 : C[16] = fC[16] + sc19;
2399 0 : C[17] = c17;
2400 :
2401 0 : C[18] = fC[18];
2402 0 : C[19] = fC[19];
2403 0 : C[20] = fC[20];
2404 0 : C[21] = fC[21] + dS*fC[24];
2405 0 : C[22] = fC[22] + dS*fC[25];
2406 0 : C[23] = fC[23] + dS*fC[26];
2407 :
2408 0 : C[24] = fC[24];
2409 0 : C[25] = fC[25];
2410 0 : C[26] = fC[26];
2411 0 : C[27] = fC[27];
2412 0 : C[28] = fC[28] + dS*fC[31];
2413 0 : C[29] = fC[29] + dS*fC[32];
2414 0 : C[30] = fC[30] + dS*fC[33];
2415 :
2416 0 : C[31] = fC[31];
2417 0 : C[32] = fC[32];
2418 0 : C[33] = fC[33];
2419 0 : C[34] = fC[34];
2420 0 : C[35] = fC[35];
2421 0 : }
2422 :
2423 :
2424 : void AliKFParticleBase::ConstructGammaBz( const AliKFParticleBase &daughter1,
2425 : const AliKFParticleBase &daughter2, double Bz )
2426 : {
2427 : //* Create gamma
2428 :
2429 0 : const AliKFParticleBase *daughters[2] = { &daughter1, &daughter2};
2430 :
2431 0 : double v0[3];
2432 :
2433 0 : if( !fIsLinearized ){
2434 0 : Double_t ds, ds1;
2435 0 : Double_t m[8];
2436 0 : Double_t mCd[36];
2437 0 : daughter1.GetDStoParticle(daughter2, ds, ds1);
2438 0 : daughter1.Transport( ds, m, mCd );
2439 0 : fP[0] = m[0];
2440 0 : fP[1] = m[1];
2441 0 : fP[2] = m[2];
2442 0 : daughter2.Transport( ds1, m, mCd );
2443 0 : fP[0] = .5*( fP[0] + m[0] );
2444 0 : fP[1] = .5*( fP[1] + m[1] );
2445 0 : fP[2] = .5*( fP[2] + m[2] );
2446 0 : } else {
2447 0 : fP[0] = fVtxGuess[0];
2448 0 : fP[1] = fVtxGuess[1];
2449 0 : fP[2] = fVtxGuess[2];
2450 : }
2451 :
2452 0 : double daughterP[2][8], daughterC[2][36];
2453 0 : double vtxMom[2][3];
2454 :
2455 0 : int nIter = fIsLinearized ?1 :2;
2456 :
2457 0 : for( int iter=0; iter<nIter; iter++){
2458 :
2459 0 : v0[0] = fP[0];
2460 0 : v0[1] = fP[1];
2461 0 : v0[2] = fP[2];
2462 :
2463 0 : fAtProductionVertex = 0;
2464 0 : fSFromDecay = 0;
2465 0 : fP[0] = v0[0];
2466 0 : fP[1] = v0[1];
2467 0 : fP[2] = v0[2];
2468 0 : fP[3] = 0;
2469 0 : fP[4] = 0;
2470 0 : fP[5] = 0;
2471 0 : fP[6] = 0;
2472 0 : fP[7] = 0;
2473 :
2474 :
2475 : // fit daughters to the vertex guess
2476 :
2477 : {
2478 0 : for( int id=0; id<2; id++ ){
2479 :
2480 0 : double *p = daughterP[id];
2481 0 : double *mC = daughterC[id];
2482 :
2483 0 : daughters[id]->GetMeasurement( v0, p, mC );
2484 :
2485 0 : Double_t mAi[6];
2486 0 : InvertSym3(mC, mAi );
2487 :
2488 : Double_t mB[3][3];
2489 :
2490 0 : mB[0][0] = mC[ 6]*mAi[0] + mC[ 7]*mAi[1] + mC[ 8]*mAi[3];
2491 0 : mB[0][1] = mC[ 6]*mAi[1] + mC[ 7]*mAi[2] + mC[ 8]*mAi[4];
2492 0 : mB[0][2] = mC[ 6]*mAi[3] + mC[ 7]*mAi[4] + mC[ 8]*mAi[5];
2493 :
2494 0 : mB[1][0] = mC[10]*mAi[0] + mC[11]*mAi[1] + mC[12]*mAi[3];
2495 0 : mB[1][1] = mC[10]*mAi[1] + mC[11]*mAi[2] + mC[12]*mAi[4];
2496 0 : mB[1][2] = mC[10]*mAi[3] + mC[11]*mAi[4] + mC[12]*mAi[5];
2497 :
2498 0 : mB[2][0] = mC[15]*mAi[0] + mC[16]*mAi[1] + mC[17]*mAi[3];
2499 0 : mB[2][1] = mC[15]*mAi[1] + mC[16]*mAi[2] + mC[17]*mAi[4];
2500 0 : mB[2][2] = mC[15]*mAi[3] + mC[16]*mAi[4] + mC[17]*mAi[5];
2501 :
2502 0 : Double_t z[3] = { v0[0]-p[0], v0[1]-p[1], v0[2]-p[2] };
2503 :
2504 0 : vtxMom[id][0] = p[3] + mB[0][0]*z[0] + mB[0][1]*z[1] + mB[0][2]*z[2];
2505 0 : vtxMom[id][1] = p[4] + mB[1][0]*z[0] + mB[1][1]*z[1] + mB[1][2]*z[2];
2506 0 : vtxMom[id][2] = p[5] + mB[2][0]*z[0] + mB[2][1]*z[1] + mB[2][2]*z[2];
2507 :
2508 0 : daughters[id]->Transport( daughters[id]->GetDStoPoint(v0), p, mC );
2509 :
2510 0 : }
2511 :
2512 : } // fit daughters to guess
2513 :
2514 :
2515 : // fit new vertex
2516 : {
2517 :
2518 0 : double mpx0 = vtxMom[0][0]+vtxMom[1][0];
2519 0 : double mpy0 = vtxMom[0][1]+vtxMom[1][1];
2520 0 : double mpt0 = TMath::Sqrt(mpx0*mpx0 + mpy0*mpy0);
2521 : // double a0 = TMath::ATan2(mpy0,mpx0);
2522 :
2523 0 : double ca0 = mpx0/mpt0;
2524 0 : double sa0 = mpy0/mpt0;
2525 0 : double r[3] = { v0[0], v0[1], v0[2] };
2526 0 : double mC[3][3] = {{1000., 0 , 0 },
2527 : {0, 1000., 0 },
2528 : {0, 0, 1000. } };
2529 : double chi2=0;
2530 :
2531 0 : for( int id=0; id<2; id++ ){
2532 : const Double_t kCLight = 0.000299792458;
2533 0 : Double_t q = Bz*daughters[id]->GetQ()*kCLight;
2534 0 : Double_t px0 = vtxMom[id][0];
2535 0 : Double_t py0 = vtxMom[id][1];
2536 0 : Double_t pz0 = vtxMom[id][2];
2537 0 : Double_t pt0 = TMath::Sqrt(px0*px0+py0*py0);
2538 0 : Double_t mG[3][6], mB[3], mH[3][3];
2539 : // r = {vx,vy,vz};
2540 : // m = {x,y,z,Px,Py,Pz};
2541 : // V = daughter.C
2542 : // G*m + B = H*r;
2543 : // q*x + Py - q*vx - sin(a)*Pt = 0
2544 : // q*y - Px - q*vy + cos(a)*Pt = 0
2545 : // (Px*cos(a) + Py*sin(a) ) (vz -z) - Pz( cos(a)*(vx-x) + sin(a)*(vy-y)) = 0
2546 :
2547 0 : mG[0][0] = q;
2548 0 : mG[0][1] = 0;
2549 0 : mG[0][2] = 0;
2550 0 : mG[0][3] = -sa0*px0/pt0;
2551 0 : mG[0][4] = 1 -sa0*py0/pt0;
2552 0 : mG[0][5] = 0;
2553 0 : mH[0][0] = q;
2554 0 : mH[0][1] = 0;
2555 0 : mH[0][2] = 0;
2556 0 : mB[0] = py0 - sa0*pt0 - mG[0][3]*px0 - mG[0][4]*py0 ;
2557 :
2558 : // q*y - Px - q*vy + cos(a)*Pt = 0
2559 :
2560 0 : mG[1][0] = 0;
2561 0 : mG[1][1] = q;
2562 0 : mG[1][2] = 0;
2563 0 : mG[1][3] = -1 + ca0*px0/pt0;
2564 0 : mG[1][4] = + ca0*py0/pt0;
2565 0 : mG[1][5] = 0;
2566 0 : mH[1][0] = 0;
2567 0 : mH[1][1] = q;
2568 0 : mH[1][2] = 0;
2569 0 : mB[1] = -px0 + ca0*pt0 - mG[1][3]*px0 - mG[1][4]*py0 ;
2570 :
2571 : // (Px*cos(a) + Py*sin(a) ) (z -vz) - Pz( cos(a)*(x-vx) + sin(a)*(y-vy)) = 0
2572 :
2573 0 : mG[2][0] = -pz0*ca0;
2574 0 : mG[2][1] = -pz0*sa0;
2575 0 : mG[2][2] = px0*ca0 + py0*sa0;
2576 0 : mG[2][3] = 0;
2577 0 : mG[2][4] = 0;
2578 0 : mG[2][5] = 0;
2579 :
2580 0 : mH[2][0] = mG[2][0];
2581 0 : mH[2][1] = mG[2][1];
2582 0 : mH[2][2] = mG[2][2];
2583 :
2584 0 : mB[2] = 0;
2585 :
2586 : // fit the vertex
2587 :
2588 : // V = GVGt
2589 :
2590 0 : double mGV[3][6];
2591 0 : double mV[6];
2592 0 : double m[3];
2593 0 : for( int i=0; i<3; i++ ){
2594 0 : m[i] = mB[i];
2595 0 : for( int k=0; k<6; k++ ) m[i]+=mG[i][k]*daughterP[id][k];
2596 : }
2597 0 : for( int i=0; i<3; i++ ){
2598 0 : for( int j=0; j<6; j++ ){
2599 0 : mGV[i][j] = 0;
2600 0 : for( int k=0; k<6; k++ ) mGV[i][j]+=mG[i][k]*daughterC[id][ IJ(k,j) ];
2601 : }
2602 : }
2603 0 : for( int i=0, k=0; i<3; i++ ){
2604 0 : for( int j=0; j<=i; j++,k++ ){
2605 0 : mV[k] = 0;
2606 0 : for( int l=0; l<6; l++ ) mV[k]+=mGV[i][l]*mG[j][l];
2607 : }
2608 : }
2609 :
2610 :
2611 : //* CHt
2612 :
2613 0 : Double_t mCHt[3][3];
2614 0 : Double_t mHCHt[6];
2615 0 : Double_t mHr[3];
2616 0 : for( int i=0; i<3; i++ ){
2617 0 : mHr[i] = 0;
2618 0 : for( int k=0; k<3; k++ ) mHr[i]+= mH[i][k]*r[k];
2619 : }
2620 :
2621 0 : for( int i=0; i<3; i++ ){
2622 0 : for( int j=0; j<3; j++){
2623 0 : mCHt[i][j] = 0;
2624 0 : for( int k=0; k<3; k++ ) mCHt[i][j]+= mC[i][k]*mH[j][k];
2625 : }
2626 : }
2627 :
2628 0 : for( int i=0, k=0; i<3; i++ ){
2629 0 : for( int j=0; j<=i; j++, k++ ){
2630 0 : mHCHt[k] = 0;
2631 0 : for( int l=0; l<3; l++ ) mHCHt[k]+= mH[i][l]*mCHt[l][j];
2632 : }
2633 : }
2634 :
2635 0 : Double_t mS[6] = { mHCHt[0]+mV[0],
2636 0 : mHCHt[1]+mV[1], mHCHt[2]+mV[2],
2637 0 : mHCHt[3]+mV[3], mHCHt[4]+mV[4], mHCHt[5]+mV[5] };
2638 :
2639 :
2640 0 : InvertSym3(mS,mS);
2641 :
2642 : //* Residual (measured - estimated)
2643 :
2644 0 : Double_t zeta[3] = { m[0]-mHr[0], m[1]-mHr[1], m[2]-mHr[2] };
2645 :
2646 : //* Kalman gain K = mCH'*S
2647 :
2648 0 : Double_t k[3][3];
2649 :
2650 0 : for(Int_t i=0;i<3;++i){
2651 0 : k[i][0] = mCHt[i][0]*mS[0] + mCHt[i][1]*mS[1] + mCHt[i][2]*mS[3];
2652 0 : k[i][1] = mCHt[i][0]*mS[1] + mCHt[i][1]*mS[2] + mCHt[i][2]*mS[4];
2653 0 : k[i][2] = mCHt[i][0]*mS[3] + mCHt[i][1]*mS[4] + mCHt[i][2]*mS[5];
2654 : }
2655 :
2656 : //* New estimation of the vertex position r += K*zeta
2657 :
2658 0 : for(Int_t i=0;i<3;++i)
2659 0 : r[i] = r[i] + k[i][0]*zeta[0] + k[i][1]*zeta[1] + k[i][2]*zeta[2];
2660 :
2661 : //* New covariance matrix C -= K*(mCH')'
2662 :
2663 0 : for(Int_t i=0;i<3;++i){
2664 0 : for(Int_t j=0;j<=i;++j){
2665 0 : mC[i][j] = mC[i][j] - (k[i][0]*mCHt[j][0] + k[i][1]*mCHt[j][1] + k[i][2]*mCHt[j][2]);
2666 0 : mC[j][i] = mC[i][j];
2667 : }
2668 : }
2669 :
2670 : //* Calculate Chi^2
2671 :
2672 0 : chi2 += ( ( mS[0]*zeta[0] + mS[1]*zeta[1] + mS[3]*zeta[2] )*zeta[0]
2673 0 : +(mS[1]*zeta[0] + mS[2]*zeta[1] + mS[4]*zeta[2] )*zeta[1]
2674 0 : +(mS[3]*zeta[0] + mS[4]*zeta[1] + mS[5]*zeta[2] )*zeta[2] );
2675 0 : }
2676 :
2677 : // store vertex
2678 :
2679 0 : fNDF = 2;
2680 0 : fChi2 = chi2;
2681 0 : for( int i=0; i<3; i++ ) fP[i] = r[i];
2682 0 : for( int i=0,k=0; i<3; i++ ){
2683 0 : for( int j=0; j<=i; j++,k++ ){
2684 0 : fC[k] = mC[i][j];
2685 : }
2686 : }
2687 0 : }
2688 :
2689 : } // iterations
2690 :
2691 : // now fit daughters to the vertex
2692 :
2693 0 : fQ = 0;
2694 0 : fSFromDecay = 0;
2695 :
2696 0 : for(Int_t i=3;i<8;++i) fP[i]=0.;
2697 0 : for(Int_t i=6;i<35;++i) fC[i]=0.;
2698 0 : fC[35] = 100.;
2699 :
2700 0 : for( int id=0; id<2; id++ ){
2701 :
2702 0 : double *p = daughterP[id];
2703 0 : double *mC = daughterC[id];
2704 0 : daughters[id]->GetMeasurement( v0, p, mC );
2705 :
2706 0 : const Double_t *m = fP, *mV = fC;
2707 :
2708 0 : Double_t mAi[6];
2709 0 : InvertSym3(mC, mAi );
2710 :
2711 : Double_t mB[4][3];
2712 :
2713 0 : mB[0][0] = mC[ 6]*mAi[0] + mC[ 7]*mAi[1] + mC[ 8]*mAi[3];
2714 0 : mB[0][1] = mC[ 6]*mAi[1] + mC[ 7]*mAi[2] + mC[ 8]*mAi[4];
2715 0 : mB[0][2] = mC[ 6]*mAi[3] + mC[ 7]*mAi[4] + mC[ 8]*mAi[5];
2716 :
2717 0 : mB[1][0] = mC[10]*mAi[0] + mC[11]*mAi[1] + mC[12]*mAi[3];
2718 0 : mB[1][1] = mC[10]*mAi[1] + mC[11]*mAi[2] + mC[12]*mAi[4];
2719 0 : mB[1][2] = mC[10]*mAi[3] + mC[11]*mAi[4] + mC[12]*mAi[5];
2720 :
2721 0 : mB[2][0] = mC[15]*mAi[0] + mC[16]*mAi[1] + mC[17]*mAi[3];
2722 0 : mB[2][1] = mC[15]*mAi[1] + mC[16]*mAi[2] + mC[17]*mAi[4];
2723 0 : mB[2][2] = mC[15]*mAi[3] + mC[16]*mAi[4] + mC[17]*mAi[5];
2724 :
2725 0 : mB[3][0] = mC[21]*mAi[0] + mC[22]*mAi[1] + mC[23]*mAi[3];
2726 0 : mB[3][1] = mC[21]*mAi[1] + mC[22]*mAi[2] + mC[23]*mAi[4];
2727 0 : mB[3][2] = mC[21]*mAi[3] + mC[22]*mAi[4] + mC[23]*mAi[5];
2728 :
2729 :
2730 0 : Double_t z[3] = { m[0]-p[0], m[1]-p[1], m[2]-p[2] };
2731 :
2732 : // {
2733 : // Double_t mAV[6] = { mC[0]-mV[0], mC[1]-mV[1], mC[2]-mV[2],
2734 : // mC[3]-mV[3], mC[4]-mV[4], mC[5]-mV[5] };
2735 : //
2736 : // Double_t mAVi[6];
2737 : // if( !InvertSym3(mAV, mAVi) ){
2738 : // Double_t dChi2 = ( +(mAVi[0]*z[0] + mAVi[1]*z[1] + mAVi[3]*z[2])*z[0]
2739 : // +(mAVi[1]*z[0] + mAVi[2]*z[1] + mAVi[4]*z[2])*z[1]
2740 : // +(mAVi[3]*z[0] + mAVi[4]*z[1] + mAVi[5]*z[2])*z[2] );
2741 : // fChi2+= TMath::Abs( dChi2 );
2742 : // }
2743 : // fNDF += 2;
2744 : // }
2745 :
2746 : //* Add the daughter momentum to the particle momentum
2747 :
2748 0 : fP[3]+= p[3] + mB[0][0]*z[0] + mB[0][1]*z[1] + mB[0][2]*z[2];
2749 0 : fP[4]+= p[4] + mB[1][0]*z[0] + mB[1][1]*z[1] + mB[1][2]*z[2];
2750 0 : fP[5]+= p[5] + mB[2][0]*z[0] + mB[2][1]*z[1] + mB[2][2]*z[2];
2751 0 : fP[6]+= p[6] + mB[3][0]*z[0] + mB[3][1]*z[1] + mB[3][2]*z[2];
2752 :
2753 : Double_t d0, d1, d2;
2754 :
2755 0 : d0= mB[0][0]*mV[0] + mB[0][1]*mV[1] + mB[0][2]*mV[3] - mC[ 6];
2756 0 : d1= mB[0][0]*mV[1] + mB[0][1]*mV[2] + mB[0][2]*mV[4] - mC[ 7];
2757 0 : d2= mB[0][0]*mV[3] + mB[0][1]*mV[4] + mB[0][2]*mV[5] - mC[ 8];
2758 :
2759 : //fC[6]+= mC[ 6] + d0;
2760 : //fC[7]+= mC[ 7] + d1;
2761 : //fC[8]+= mC[ 8] + d2;
2762 0 : fC[9]+= mC[ 9] + d0*mB[0][0] + d1*mB[0][1] + d2*mB[0][2];
2763 :
2764 0 : d0= mB[1][0]*mV[0] + mB[1][1]*mV[1] + mB[1][2]*mV[3] - mC[10];
2765 0 : d1= mB[1][0]*mV[1] + mB[1][1]*mV[2] + mB[1][2]*mV[4] - mC[11];
2766 0 : d2= mB[1][0]*mV[3] + mB[1][1]*mV[4] + mB[1][2]*mV[5] - mC[12];
2767 :
2768 : //fC[10]+= mC[10]+ d0;
2769 : //fC[11]+= mC[11]+ d1;
2770 : //fC[12]+= mC[12]+ d2;
2771 0 : fC[13]+= mC[13]+ d0*mB[0][0] + d1*mB[0][1] + d2*mB[0][2];
2772 0 : fC[14]+= mC[14]+ d0*mB[1][0] + d1*mB[1][1] + d2*mB[1][2];
2773 :
2774 0 : d0= mB[2][0]*mV[0] + mB[2][1]*mV[1] + mB[2][2]*mV[3] - mC[15];
2775 0 : d1= mB[2][0]*mV[1] + mB[2][1]*mV[2] + mB[2][2]*mV[4] - mC[16];
2776 0 : d2= mB[2][0]*mV[3] + mB[2][1]*mV[4] + mB[2][2]*mV[5] - mC[17];
2777 :
2778 : //fC[15]+= mC[15]+ d0;
2779 : //fC[16]+= mC[16]+ d1;
2780 : //fC[17]+= mC[17]+ d2;
2781 0 : fC[18]+= mC[18]+ d0*mB[0][0] + d1*mB[0][1] + d2*mB[0][2];
2782 0 : fC[19]+= mC[19]+ d0*mB[1][0] + d1*mB[1][1] + d2*mB[1][2];
2783 0 : fC[20]+= mC[20]+ d0*mB[2][0] + d1*mB[2][1] + d2*mB[2][2];
2784 :
2785 0 : d0= mB[3][0]*mV[0] + mB[3][1]*mV[1] + mB[3][2]*mV[3] - mC[21];
2786 0 : d1= mB[3][0]*mV[1] + mB[3][1]*mV[2] + mB[3][2]*mV[4] - mC[22];
2787 0 : d2= mB[3][0]*mV[3] + mB[3][1]*mV[4] + mB[3][2]*mV[5] - mC[23];
2788 :
2789 : //fC[21]+= mC[21] + d0;
2790 : //fC[22]+= mC[22] + d1;
2791 : //fC[23]+= mC[23] + d2;
2792 0 : fC[24]+= mC[24] + d0*mB[0][0] + d1*mB[0][1] + d2*mB[0][2];
2793 0 : fC[25]+= mC[25] + d0*mB[1][0] + d1*mB[1][1] + d2*mB[1][2];
2794 0 : fC[26]+= mC[26] + d0*mB[2][0] + d1*mB[2][1] + d2*mB[2][2];
2795 0 : fC[27]+= mC[27] + d0*mB[3][0] + d1*mB[3][1] + d2*mB[3][2];
2796 0 : }
2797 :
2798 : // SetMassConstraint(0,0);
2799 0 : SetNonlinearMassConstraint(0);
2800 0 : }
2801 :
2802 : void AliKFParticleBase::GetArmenterosPodolanski(AliKFParticleBase& positive, AliKFParticleBase& negative, Double_t QtAlfa[2] )
2803 : {
2804 : // example:
2805 : // AliKFParticle PosParticle(...)
2806 : // AliKFParticle NegParticle(...)
2807 : // Gamma.ConstructGamma(PosParticle, NegParticle);
2808 : // Double_t VertexGamma[3] = {Gamma.GetX(), Gamma.GetY(), Gamma.GetZ()};
2809 : // PosParticle.TransportToPoint(VertexGamma);
2810 : // NegParticle.TransportToPoint(VertexGamma);
2811 : // Double_t armenterosQtAlfa[2] = {0.};
2812 : // AliKFParticle::GetArmenterosPodolanski(PosParticle, NegParticle, armenterosQtAlfa );
2813 :
2814 : Double_t alpha = 0., qt = 0.;
2815 0 : Double_t spx = positive.GetPx() + negative.GetPx();
2816 0 : Double_t spy = positive.GetPy() + negative.GetPy();
2817 0 : Double_t spz = positive.GetPz() + negative.GetPz();
2818 0 : Double_t sp = sqrt(spx*spx + spy*spy + spz*spz);
2819 0 : if( sp == 0.0) return;
2820 : Double_t pn, pln, plp; // ,pp;
2821 :
2822 0 : pn = TMath::Sqrt(negative.GetPx()*negative.GetPx() + negative.GetPy()*negative.GetPy() + negative.GetPz()*negative.GetPz());
2823 : // pp = TMath::Sqrt(positive.GetPx()*positive.GetPx() + positive.GetPy()*positive.GetPy() + positive.GetPz()*positive.GetPz());
2824 0 : pln = (negative.GetPx()*spx+negative.GetPy()*spy+negative.GetPz()*spz)/sp;
2825 0 : plp = (positive.GetPx()*spx+positive.GetPy()*spy+positive.GetPz()*spz)/sp;
2826 :
2827 0 : if( pn == 0.0) return;
2828 0 : Double_t ptm = (1.-((pln/pn)*(pln/pn)));
2829 0 : qt= (ptm>=0.)? pn*sqrt(ptm) :0;
2830 0 : alpha = (plp-pln)/(plp+pln);
2831 :
2832 0 : QtAlfa[0] = qt;
2833 0 : QtAlfa[1] = alpha;
2834 0 : }
2835 :
2836 : void AliKFParticleBase::RotateXY(Double_t angle, Double_t Vtx[3])
2837 : {
2838 : // Rotates the KFParticle object around OZ axis, OZ axis is set by the vertex position
2839 : // Double_t angle - angle of rotation in XY plane in [rad]
2840 : // Double_t Vtx[3] - position of the vertex in [cm]
2841 :
2842 : // Before rotation the center of the coordinat system should be moved to the vertex position; move back after rotation
2843 0 : X() = X() - Vtx[0];
2844 0 : Y() = Y() - Vtx[1];
2845 0 : Z() = Z() - Vtx[2];
2846 :
2847 : // Rotate the kf particle
2848 0 : Double_t c = TMath::Cos(angle);
2849 0 : Double_t s = TMath::Sin(angle);
2850 :
2851 0 : Double_t mA[8][ 8];
2852 0 : for( Int_t i=0; i<8; i++ ){
2853 0 : for( Int_t j=0; j<8; j++){
2854 0 : mA[i][j] = 0;
2855 : }
2856 : }
2857 0 : for( int i=0; i<8; i++ ){
2858 0 : mA[i][i] = 1;
2859 : }
2860 0 : mA[0][0] = c; mA[0][1] = s;
2861 0 : mA[1][0] = -s; mA[1][1] = c;
2862 0 : mA[3][3] = c; mA[3][4] = s;
2863 0 : mA[4][3] = -s; mA[4][4] = c;
2864 :
2865 0 : Double_t mAC[8][8];
2866 0 : Double_t mAp[8];
2867 :
2868 0 : for( Int_t i=0; i<8; i++ ){
2869 0 : mAp[i] = 0;
2870 0 : for( Int_t k=0; k<8; k++){
2871 0 : mAp[i]+=mA[i][k] * fP[k];
2872 : }
2873 : }
2874 :
2875 0 : for( Int_t i=0; i<8; i++){
2876 0 : fP[i] = mAp[i];
2877 : }
2878 :
2879 0 : for( Int_t i=0; i<8; i++ ){
2880 0 : for( Int_t j=0; j<8; j++ ){
2881 0 : mAC[i][j] = 0;
2882 0 : for( Int_t k=0; k<8; k++ ){
2883 0 : mAC[i][j]+= mA[i][k] * GetCovariance(k,j);
2884 : }
2885 : }
2886 : }
2887 :
2888 0 : for( Int_t i=0; i<8; i++ ){
2889 0 : for( Int_t j=0; j<=i; j++ ){
2890 : Double_t xx = 0;
2891 0 : for( Int_t k=0; k<8; k++){
2892 0 : xx+= mAC[i][k]*mA[j][k];
2893 : }
2894 0 : Covariance(i,j) = xx;
2895 : }
2896 : }
2897 :
2898 0 : X() = GetX() + Vtx[0];
2899 0 : Y() = GetY() + Vtx[1];
2900 0 : Z() = GetZ() + Vtx[2];
2901 0 : }
2902 :
2903 : Bool_t AliKFParticleBase::InvertSym3( const Double_t A[], Double_t Ai[] )
2904 : {
2905 : //* Invert symmetric matric stored in low-triagonal form
2906 :
2907 : bool ret = 0;
2908 296 : double a0 = A[0], a1 = A[1], a2 = A[2], a3 = A[3];
2909 :
2910 148 : Ai[0] = a2*A[5] - A[4]*A[4];
2911 148 : Ai[1] = a3*A[4] - a1*A[5];
2912 148 : Ai[3] = a1*A[4] - a2*a3;
2913 148 : Double_t det = (a0*Ai[0] + a1*Ai[1] + a3*Ai[3]);
2914 296 : if( TMath::Abs(det)>1.e-20 ) det = 1./det;
2915 : else{
2916 : det = 0;
2917 : ret = 1;
2918 : }
2919 148 : Ai[0] *= det;
2920 148 : Ai[1] *= det;
2921 148 : Ai[3] *= det;
2922 148 : Ai[2] = ( a0*A[5] - a3*a3 )*det;
2923 148 : Ai[4] = ( a1*a3 - a0*A[4] )*det;
2924 148 : Ai[5] = ( a0*a2 - a1*a1 )*det;
2925 148 : return ret;
2926 : }
2927 :
2928 : void AliKFParticleBase::MultQSQt( const Double_t Q[], const Double_t S[], Double_t SOut[] )
2929 : {
2930 : //* Matrix multiplication Q*S*Q^T, Q - square matrix, S - symmetric
2931 :
2932 : const Int_t kN= 8;
2933 0 : Double_t mA[kN*kN];
2934 :
2935 0 : for( Int_t i=0, ij=0; i<kN; i++ ){
2936 0 : for( Int_t j=0; j<kN; j++, ++ij ){
2937 0 : mA[ij] = 0 ;
2938 0 : for( Int_t k=0; k<kN; ++k ) mA[ij]+= S[( k<=i ) ? i*(i+1)/2+k :k*(k+1)/2+i] * Q[ j*kN+k];
2939 : }
2940 : }
2941 :
2942 0 : for( Int_t i=0; i<kN; i++ ){
2943 0 : for( Int_t j=0; j<=i; j++ ){
2944 0 : Int_t ij = ( j<=i ) ? i*(i+1)/2+j :j*(j+1)/2+i;
2945 0 : SOut[ij] = 0 ;
2946 0 : for( Int_t k=0; k<kN; k++ ) SOut[ij] += Q[ i*kN+k ] * mA[ k*kN+j ];
2947 : }
2948 : }
2949 0 : }
2950 :
2951 :
2952 : // 72-charachters line to define the printer border
2953 : //3456789012345678901234567890123456789012345678901234567890123456789012
2954 :
|
