Coding conventions fixed
[u/mrichter/AliRoot.git] / STEER / ESD / AliESDv0.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 //-------------------------------------------------------------------------
19 //               Implementation of the ESD V0 vertex class
20 //            This class is part of the Event Data Summary
21 //            set of classes and contains information about
22 //            V0 kind vertexes generated by a neutral particle
23 //     Origin: Iouri Belikov, IReS, Strasbourg, Jouri.Belikov@cern.ch
24 //     Modified by: Marian Ivanov,  CERN, Marian.Ivanov@cern.ch
25 //            and  Boris Hippolyte,IPHC, hippolyt@in2p3.fr 
26 //-------------------------------------------------------------------------
27
28 #include <TMath.h>
29 #include <TDatabasePDG.h>
30 #include <TParticlePDG.h>
31 #include <TVector3.h>
32
33 #include "AliLog.h"
34 #include "AliESDv0.h"
35 #include "AliESDV0Params.h"
36 #include "AliKFParticle.h"
37 #include "AliKFVertex.h"
38
39 ClassImp(AliESDv0)
40
41 const AliESDV0Params  AliESDv0::fgkParams;
42
43 AliESDv0::AliESDv0() :
44   AliVParticle(),
45   fParamN(),
46   fParamP(),
47   fEffMass(TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass()),
48   fDcaV0Daughters(0),
49   fChi2V0(0.),
50   fRr(0),
51   fDistSigma(0),
52   fChi2Before(0),
53   fChi2After(0),
54   fPointAngleFi(0),
55   fPointAngleTh(0),
56   fPointAngle(0),
57   fPdgCode(kK0Short),
58   fNidx(0),
59   fPidx(0),
60   fStatus(0),
61   fNBefore(0),
62   fNAfter(0),
63   fOnFlyStatus(kFALSE)
64 {
65   //--------------------------------------------------------------------
66   // Default constructor  (K0s)
67   //--------------------------------------------------------------------
68
69   for (Int_t i=0; i<3; i++) {
70     fPos[i] = 0.;
71     fNmom[i] = 0.;
72     fPmom[i] = 0.;
73   }
74
75   for (Int_t i=0; i<6; i++) {
76     fPosCov[i]= 0.;
77   }
78
79   for (Int_t i=0;i<6;i++){fClusters[0][i]=0; fClusters[1][i]=0;}
80   fNormDCAPrim[0]=fNormDCAPrim[1]=0;
81   for (Int_t i=0;i<3;i++){fAngle[i]=0;}
82   for (Int_t i=0;i<4;i++){fCausality[i]=0;}
83 }
84
85 AliESDv0::AliESDv0(const AliESDv0& v0) :
86   AliVParticle(v0),
87   fParamN(v0.fParamN),
88   fParamP(v0.fParamP),
89   fEffMass(v0.fEffMass),
90   fDcaV0Daughters(v0.fDcaV0Daughters),
91   fChi2V0(v0.fChi2V0),
92   fRr(v0.fRr),
93   fDistSigma(v0.fDistSigma),
94   fChi2Before(v0.fChi2Before),
95   fChi2After(v0.fChi2After),
96   fPointAngleFi(v0.fPointAngleFi),
97   fPointAngleTh(v0.fPointAngleTh),
98   fPointAngle(v0.fPointAngle),
99   fPdgCode(v0.fPdgCode),
100   fNidx(v0.fNidx),
101   fPidx(v0.fPidx),
102   fStatus(v0.fStatus),
103   fNBefore(v0.fNBefore),
104   fNAfter(v0.fNAfter),
105   fOnFlyStatus(v0.fOnFlyStatus)
106 {
107   //--------------------------------------------------------------------
108   // The copy constructor
109   //--------------------------------------------------------------------
110
111   for (int i=0; i<3; i++) {
112     fPos[i]  = v0.fPos[i];
113     fNmom[i] = v0.fNmom[i];
114     fPmom[i] = v0.fPmom[i];
115   }
116   for (int i=0; i<6; i++) {
117     fPosCov[i]  = v0.fPosCov[i];
118   }
119
120   for (Int_t i=0; i<2; i++) {
121     fNormDCAPrim[i]=v0.fNormDCAPrim[i];
122   }
123   for (Int_t i=0;i<6;i++){
124       fClusters[0][i]=v0.fClusters[0][i]; 
125       fClusters[1][i]=v0.fClusters[1][i];
126   }
127   for (Int_t i=0;i<3;i++){
128       fAngle[i]=v0.fAngle[i];
129   }
130   for (Int_t i=0;i<4;i++){fCausality[i]=v0.fCausality[i];}
131 }
132
133
134 AliESDv0::AliESDv0(const AliExternalTrackParam &t1, Int_t i1,
135                    const AliExternalTrackParam &t2, Int_t i2) :
136   AliVParticle(),
137   fParamN(t1),
138   fParamP(t2),
139   fEffMass(TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass()),
140   fDcaV0Daughters(0),
141   fChi2V0(0.),
142   fRr(0),
143   fDistSigma(0),
144   fChi2Before(0),
145   fChi2After(0),
146   fPointAngleFi(0),
147   fPointAngleTh(0),
148   fPointAngle(0),
149   fPdgCode(kK0Short),
150   fNidx(i1),
151   fPidx(i2),
152   fStatus(0),
153   fNBefore(0),
154   fNAfter(0),
155   fOnFlyStatus(kFALSE)
156 {
157   //--------------------------------------------------------------------
158   // Main constructor  (K0s)
159   //--------------------------------------------------------------------
160
161   //Make sure the daughters are ordered (needed for the on-the-fly V0s)
162   Short_t cN=t1.Charge(), cP=t2.Charge();
163   if ((cN>0) && (cN != cP)) {
164      fParamN.~AliExternalTrackParam();
165      new (&fParamN) AliExternalTrackParam(t2);
166      fParamP.~AliExternalTrackParam();
167      new (&fParamP) AliExternalTrackParam(t1);
168
169      Int_t index=fNidx;
170      fNidx=fPidx;
171      fPidx=index;
172   }
173
174   for (Int_t i=0; i<6; i++) {
175     fPosCov[i]= 0.;
176   }
177
178   //Trivial estimation of the vertex parameters
179   Double_t alpha=t1.GetAlpha(), cs=TMath::Cos(alpha), sn=TMath::Sin(alpha);
180   Double_t tmp[3];
181   t1.GetPxPyPz(tmp);
182   Double_t px1=tmp[0], py1=tmp[1], pz1=tmp[2];
183   t1.GetXYZ(tmp);
184   Double_t  x1=tmp[0],  y1=tmp[1],  z1=tmp[2];
185   const Double_t ss=0.0005*0.0005;//a kind of a residual misalignment precision
186   Double_t sx1=sn*sn*t1.GetSigmaY2()+ss, sy1=cs*cs*t1.GetSigmaY2()+ss; 
187
188
189   alpha=t2.GetAlpha(); cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
190   t2.GetPxPyPz(tmp);
191   Double_t px2=tmp[0], py2=tmp[1], pz2=tmp[2];
192   t2.GetXYZ(tmp);
193   Double_t  x2=tmp[0],  y2=tmp[1],  z2=tmp[2];
194   Double_t sx2=sn*sn*t2.GetSigmaY2()+ss, sy2=cs*cs*t2.GetSigmaY2()+ss; 
195     
196   Double_t sz1=t1.GetSigmaZ2(), sz2=t2.GetSigmaZ2();
197   Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
198   Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
199   Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
200   fPos[0]=wx1*x1 + wx2*x2; fPos[1]=wy1*y1 + wy2*y2; fPos[2]=wz1*z1 + wz2*z2;
201
202   //fPos[0]=0.5*(x1+x2); fPos[1]=0.5*(y1+y2); fPos[2]=0.5*(z1+z2);
203   fNmom[0]=px1; fNmom[1]=py1; fNmom[2]=pz1; 
204   fPmom[0]=px2; fPmom[1]=py2; fPmom[2]=pz2;
205
206   for (Int_t i=0;i<6;i++){fClusters[0][i]=0; fClusters[1][i]=0;}
207   fNormDCAPrim[0]=fNormDCAPrim[1]=0;
208   for (Int_t i=0;i<3;i++){fAngle[i]=0;}
209   for (Int_t i=0;i<4;i++){fCausality[i]=0;}
210 }
211
212 AliESDv0& AliESDv0::operator=(const AliESDv0 &v0)
213 {
214   //--------------------------------------------------------------------
215   // The assignment operator
216   //--------------------------------------------------------------------
217
218   if(this==&v0)return *this;
219   AliVParticle::operator=(v0);
220   fParamN  = v0.fParamN;
221   fParamP  = v0.fParamP;
222   fEffMass = v0.fEffMass;
223   fDcaV0Daughters = v0.fDcaV0Daughters;
224   fChi2V0 = v0.fChi2V0;
225   fRr = v0.fRr;
226   fDistSigma    = v0.fDistSigma;
227   fChi2Before   = v0.fChi2Before;
228   fChi2After    = v0.fChi2After;
229   fPointAngleFi = v0.fPointAngleFi;
230   fPointAngleTh = v0.fPointAngleTh;
231   fPointAngle   = v0.fPointAngle;
232   fPdgCode      = v0.fPdgCode;
233   fNidx         = v0.fNidx;
234   fPidx         = v0.fPidx;
235   fStatus       = v0.fStatus;
236   fNBefore      = v0.fNBefore;
237   fNAfter       = v0.fNAfter;
238   fOnFlyStatus  = v0.fOnFlyStatus;
239
240   for (int i=0; i<3; i++) {
241     fPos[i]  = v0.fPos[i];
242     fNmom[i] = v0.fNmom[i];
243     fPmom[i] = v0.fPmom[i];
244   }
245   for (int i=0; i<6; i++) {
246     fPosCov[i]  = v0.fPosCov[i];
247   }
248   for (Int_t i=0; i<2; i++) {
249     fNormDCAPrim[i]=v0.fNormDCAPrim[i];
250   }
251   for (Int_t i=0;i<6;i++){
252       fClusters[0][i]=v0.fClusters[0][i]; 
253       fClusters[1][i]=v0.fClusters[1][i];
254   }
255   for (Int_t i=0;i<3;i++){
256       fAngle[i]=v0.fAngle[i];
257   }
258   for (Int_t i=0;i<4;i++){fCausality[i]=v0.fCausality[i];}
259
260   return *this;
261 }
262
263 void AliESDv0::Copy(TObject& obj) const {
264
265   // this overwrites the virtual TOBject::Copy()
266   // to allow run time copying without casting
267   // in AliESDEvent
268
269   if(this==&obj)return;
270   AliESDv0 *robj = dynamic_cast<AliESDv0*>(&obj);
271   if(!robj)return; // not an aliesesv0
272   *robj = *this;
273 }
274
275 AliESDv0::~AliESDv0(){
276   //--------------------------------------------------------------------
277   // Empty destructor
278   //--------------------------------------------------------------------
279 }
280
281 // Start with AliVParticle functions
282 Double_t AliESDv0::E() const {
283   //--------------------------------------------------------------------
284   // This gives the energy assuming the ChangeMassHypothesis was called
285   //--------------------------------------------------------------------
286   return E(fPdgCode);
287 }
288
289 Double_t AliESDv0::Y() const {
290   //--------------------------------------------------------------------
291   // This gives the energy assuming the ChangeMassHypothesis was called
292   //--------------------------------------------------------------------
293   return Y(fPdgCode);
294 }
295
296 // Then extend AliVParticle functions
297 Double_t AliESDv0::E(Int_t pdg) const {
298   //--------------------------------------------------------------------
299   // This gives the energy with the particle hypothesis as argument 
300   //--------------------------------------------------------------------
301   Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
302   return TMath::Sqrt(mass*mass+P()*P());
303 }
304
305 Double_t AliESDv0::Y(Int_t pdg) const {
306   //--------------------------------------------------------------------
307   // This gives the rapidity with the particle hypothesis as argument 
308   //--------------------------------------------------------------------
309   return 0.5*TMath::Log((E(pdg)+Pz())/(E(pdg)-Pz()+1.e-13));
310 }
311
312 // Now the functions for analysis consistency
313 Double_t AliESDv0::RapK0Short() const {
314   //--------------------------------------------------------------------
315   // This gives the pseudorapidity assuming a K0s particle
316   //--------------------------------------------------------------------
317   return Y(kK0Short);
318 }
319
320 Double_t AliESDv0::RapLambda() const {
321   //--------------------------------------------------------------------
322   // This gives the pseudorapidity assuming a (Anti) Lambda particle
323   //--------------------------------------------------------------------
324   return Y(kLambda0);
325 }
326
327 Double_t AliESDv0::AlphaV0() const {
328   //--------------------------------------------------------------------
329   // This gives the Armenteros-Podolanski alpha
330   //--------------------------------------------------------------------
331   TVector3 momNeg(fNmom[0],fNmom[1],fNmom[2]);
332   TVector3 momPos(fPmom[0],fPmom[1],fPmom[2]);
333   TVector3 momTot(Px(),Py(),Pz());
334
335   Double_t lQlNeg = momNeg.Dot(momTot)/momTot.Mag();
336   Double_t lQlPos = momPos.Dot(momTot)/momTot.Mag();
337
338   //return 1.-2./(1.+lQlNeg/lQlPos);
339   return (lQlPos - lQlNeg)/(lQlPos + lQlNeg);
340 }
341
342 Double_t AliESDv0::PtArmV0() const {
343   //--------------------------------------------------------------------
344   // This gives the Armenteros-Podolanski ptarm
345   //--------------------------------------------------------------------
346   TVector3 momNeg(fNmom[0],fNmom[1],fNmom[2]);
347   TVector3 momTot(Px(),Py(),Pz());
348
349   return momNeg.Perp(momTot);
350 }
351
352 // Eventually the older functions
353 Double_t AliESDv0::ChangeMassHypothesis(Int_t code) {
354   //--------------------------------------------------------------------
355   // This function changes the mass hypothesis for this V0
356   // and returns the "kinematical quality" of this hypothesis 
357   //--------------------------------------------------------------------
358   static
359   Double_t piMass=TDatabasePDG::Instance()->GetParticle(kPiPlus)->Mass();
360   static
361   Double_t prMass=TDatabasePDG::Instance()->GetParticle(kProton)->Mass();
362   static
363   Double_t k0Mass=TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
364   static
365   Double_t l0Mass=TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
366
367   Double_t nmass=piMass, pmass=piMass, mass=k0Mass, ps=0.206;
368
369   fPdgCode=code;
370
371   switch (code) {
372   case kLambda0:
373     nmass=piMass; pmass=prMass; mass=l0Mass; ps=0.101; break;
374   case kLambda0Bar:
375     pmass=piMass; nmass=prMass; mass=l0Mass; ps=0.101; break;
376   case kK0Short: 
377     break;
378   default:
379     AliError("invalide PDG code ! Assuming K0s...");
380     fPdgCode=kK0Short;
381     break;
382   }
383
384   Double_t pxn=fNmom[0], pyn=fNmom[1], pzn=fNmom[2]; 
385   Double_t pxp=fPmom[0], pyp=fPmom[1], pzp=fPmom[2];
386
387   Double_t en=TMath::Sqrt(nmass*nmass + pxn*pxn + pyn*pyn + pzn*pzn);
388   Double_t ep=TMath::Sqrt(pmass*pmass + pxp*pxp + pyp*pyp + pzp*pzp);
389   Double_t pxl=pxn+pxp, pyl=pyn+pyp, pzl=pzn+pzp;
390   Double_t pl=TMath::Sqrt(pxl*pxl + pyl*pyl + pzl*pzl);
391
392   fEffMass=TMath::Sqrt((en+ep)*(en+ep)-pl*pl);
393
394   Double_t beta=pl/(en+ep);
395   Double_t pln=(pxn*pxl + pyn*pyl + pzn*pzl)/pl;
396   Double_t plp=(pxp*pxl + pyp*pyl + pzp*pzl)/pl;
397
398   Double_t pt2=pxp*pxp + pyp*pyp + pzp*pzp - plp*plp;
399
400   Double_t a=(plp-pln)/(plp+pln);
401   a -= (pmass*pmass-nmass*nmass)/(mass*mass);
402   a = 0.25*beta*beta*mass*mass*a*a + pt2;
403
404   return (a - ps*ps);
405   
406 }
407
408 void AliESDv0::GetPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const {
409   //--------------------------------------------------------------------
410   // This function returns V0's momentum (global)
411   //--------------------------------------------------------------------
412   px=fNmom[0]+fPmom[0]; 
413   py=fNmom[1]+fPmom[1]; 
414   pz=fNmom[2]+fPmom[2]; 
415 }
416
417 void AliESDv0::GetXYZ(Double_t &x, Double_t &y, Double_t &z) const {
418   //--------------------------------------------------------------------
419   // This function returns V0's position (global)
420   //--------------------------------------------------------------------
421   x=fPos[0]; 
422   y=fPos[1]; 
423   z=fPos[2]; 
424 }
425
426 Float_t AliESDv0::GetD(Double_t x0, Double_t y0) const {
427   //--------------------------------------------------------------------
428   // This function returns V0's impact parameter calculated in 2D in XY plane
429   //--------------------------------------------------------------------
430   Double_t x=fPos[0],y=fPos[1];
431   Double_t px=fNmom[0]+fPmom[0];
432   Double_t py=fNmom[1]+fPmom[1];
433
434   Double_t dz=(x0-x)*py - (y0-y)*px;
435   Double_t d=TMath::Sqrt(dz*dz/(px*px+py*py));
436   return d;
437 }
438
439 Float_t AliESDv0::GetD(Double_t x0, Double_t y0, Double_t z0) const {
440   //--------------------------------------------------------------------
441   // This function returns V0's impact parameter calculated in 3D
442   //--------------------------------------------------------------------
443   Double_t x=fPos[0],y=fPos[1],z=fPos[2];
444   Double_t px=fNmom[0]+fPmom[0];
445   Double_t py=fNmom[1]+fPmom[1];
446   Double_t pz=fNmom[2]+fPmom[2];
447
448   Double_t dx=(y0-y)*pz - (z0-z)*py; 
449   Double_t dy=(x0-x)*pz - (z0-z)*px;
450   Double_t dz=(x0-x)*py - (y0-y)*px;
451   Double_t d=TMath::Sqrt((dx*dx+dy*dy+dz*dz)/(px*px+py*py+pz*pz));
452   return d;
453 }
454
455 Float_t AliESDv0::GetV0CosineOfPointingAngle(Double_t refPointX, Double_t refPointY, Double_t refPointZ) const {
456   // calculates the pointing angle of the V0 wrt a reference point
457
458   Double_t momV0[3]; //momentum of the V0
459   GetPxPyPz(momV0[0],momV0[1],momV0[2]);
460
461   Double_t deltaPos[3]; //vector between the reference point and the V0 vertex
462   deltaPos[0] = fPos[0] - refPointX;
463   deltaPos[1] = fPos[1] - refPointY;
464   deltaPos[2] = fPos[2] - refPointZ;
465
466   Double_t momV02    = momV0[0]*momV0[0] + momV0[1]*momV0[1] + momV0[2]*momV0[2];
467   Double_t deltaPos2 = deltaPos[0]*deltaPos[0] + deltaPos[1]*deltaPos[1] + deltaPos[2]*deltaPos[2];
468
469   Double_t cosinePointingAngle = (deltaPos[0]*momV0[0] +
470                                   deltaPos[1]*momV0[1] +
471                                   deltaPos[2]*momV0[2] ) /
472     TMath::Sqrt(momV02 * deltaPos2);
473   
474   return cosinePointingAngle;
475 }
476
477
478 // **** The following functions need to be revised
479
480 void AliESDv0::GetPosCov(Double_t cov[6]) const {
481
482   for (Int_t i=0; i<6; ++i) cov[i] = fPosCov[i];
483
484 }
485
486 Double_t AliESDv0::GetSigmaY(){
487   //
488   // return sigmay in y  at vertex position  using covariance matrix 
489   //
490   const Double_t * cp  = fParamP.GetCovariance();
491   const Double_t * cm  = fParamN.GetCovariance();
492   Double_t sigmay = cp[0]+cm[0]+ cp[5]*(fParamP.GetX()-fRr)*(fParamP.GetX()-fRr)+ cm[5]*(fParamN.GetX()-fRr)*(fParamN.GetX()-fRr);
493   return (sigmay>0) ? TMath::Sqrt(sigmay):100;
494 }
495
496 Double_t AliESDv0::GetSigmaZ(){
497   //
498   // return sigmay in y  at vertex position  using covariance matrix 
499   //
500   const Double_t * cp  = fParamP.GetCovariance();
501   const Double_t * cm  = fParamN.GetCovariance();
502   Double_t sigmaz = cp[2]+cm[2]+ cp[9]*(fParamP.GetX()-fRr)*(fParamP.GetX()-fRr)+ cm[9]*(fParamN.GetX()-fRr)*(fParamN.GetX()-fRr);
503   return (sigmaz>0) ? TMath::Sqrt(sigmaz):100;
504 }
505
506 Double_t AliESDv0::GetSigmaD0(){
507   //
508   // Sigma parameterization using covariance matrix
509   //
510   // sigma of distance between two tracks in vertex position 
511   // sigma of DCA is proportianal to sigmaD0
512   // factor 2 difference is explained by the fact that the DCA is calculated at the position 
513   // where the tracks as closest together ( not exact position of the vertex)
514   //
515   const Double_t * cp      = fParamP.GetCovariance();
516   const Double_t * cm      = fParamN.GetCovariance();
517   Double_t sigmaD0   = cp[0]+cm[0]+cp[2]+cm[2]+fgkParams.fPSigmaOffsetD0*fgkParams.fPSigmaOffsetD0;
518   sigmaD0           += ((fParamP.GetX()-fRr)*(fParamP.GetX()-fRr))*(cp[5]+cp[9]);
519   sigmaD0           += ((fParamN.GetX()-fRr)*(fParamN.GetX()-fRr))*(cm[5]+cm[9]);
520   return (sigmaD0>0)? TMath::Sqrt(sigmaD0):100;
521 }
522
523
524 Double_t AliESDv0::GetSigmaAP0(){
525   //
526   //Sigma parameterization using covariance matrices
527   //
528   Double_t prec  = TMath::Sqrt((fNmom[0]+fPmom[0])*(fNmom[0]+fPmom[0])
529                               +(fNmom[1]+fPmom[1])*(fNmom[1]+fPmom[1])
530                               +(fNmom[2]+fPmom[2])*(fNmom[2]+fPmom[2]));
531   Double_t normp = TMath::Sqrt(fPmom[0]*fPmom[0]+fPmom[1]*fPmom[1]+fPmom[2]*fPmom[2])/prec;  // fraction of the momenta
532   Double_t normm = TMath::Sqrt(fNmom[0]*fNmom[0]+fNmom[1]*fNmom[1]+fNmom[2]*fNmom[2])/prec;  
533   const Double_t * cp      = fParamP.GetCovariance();
534   const Double_t * cm      = fParamN.GetCovariance();
535   Double_t sigmaAP0 = fgkParams.fPSigmaOffsetAP0*fgkParams.fPSigmaOffsetAP0;                           // minimal part
536   sigmaAP0 +=  (cp[5]+cp[9])*(normp*normp)+(cm[5]+cm[9])*(normm*normm);          // angular resolution part
537   Double_t sigmaAP1 = GetSigmaD0()/(TMath::Abs(fRr)+0.01);                       // vertex position part
538   sigmaAP0 +=  0.5*sigmaAP1*sigmaAP1;                              
539   return (sigmaAP0>0)? TMath::Sqrt(sigmaAP0):100;
540 }
541
542 Double_t AliESDv0::GetEffectiveSigmaD0(){
543   //
544   // minimax - effective Sigma parameterization 
545   // p12 effective curvature and v0 radius postion used as parameters  
546   //  
547   Double_t p12 = TMath::Sqrt(fParamP.GetParameter()[4]*fParamP.GetParameter()[4]+
548                              fParamN.GetParameter()[4]*fParamN.GetParameter()[4]);
549   Double_t sigmaED0= TMath::Max(TMath::Sqrt(fRr)-fgkParams.fPSigmaRminDE,0.0)*fgkParams.fPSigmaCoefDE*p12*p12;
550   sigmaED0*= sigmaED0;
551   sigmaED0*= sigmaED0;
552   sigmaED0 = TMath::Sqrt(sigmaED0+fgkParams.fPSigmaOffsetDE*fgkParams.fPSigmaOffsetDE);
553   return (sigmaED0<fgkParams.fPSigmaMaxDE) ? sigmaED0: fgkParams.fPSigmaMaxDE;
554 }
555
556
557 Double_t AliESDv0::GetEffectiveSigmaAP0(){
558   //
559   // effective Sigma parameterization of point angle resolution 
560   //
561   Double_t p12 = TMath::Sqrt(fParamP.GetParameter()[4]*fParamP.GetParameter()[4]+
562                              fParamN.GetParameter()[4]*fParamN.GetParameter()[4]);
563   Double_t sigmaAPE= fgkParams.fPSigmaBase0APE;
564   sigmaAPE+= fgkParams.fPSigmaR0APE/(fgkParams.fPSigmaR1APE+fRr);
565   sigmaAPE*= (fgkParams.fPSigmaP0APE+fgkParams.fPSigmaP1APE*p12);
566   sigmaAPE = TMath::Min(sigmaAPE,fgkParams.fPSigmaMaxAPE);
567   return sigmaAPE;
568 }
569
570
571 Double_t  AliESDv0::GetMinimaxSigmaAP0(){
572   //
573   // calculate mini-max effective sigma of point angle resolution
574   //
575   //compv0->fTree->SetAlias("SigmaAP2","max(min((SigmaAP0+SigmaAPE0)*0.5,1.5*SigmaAPE0),0.5*SigmaAPE0+0.003)");
576   Double_t    effectiveSigma = GetEffectiveSigmaAP0();
577   Double_t    sigmaMMAP = 0.5*(GetSigmaAP0()+effectiveSigma);
578   sigmaMMAP  = TMath::Min(sigmaMMAP, fgkParams.fPMaxFractionAP0*effectiveSigma);
579   sigmaMMAP  = TMath::Max(sigmaMMAP, fgkParams.fPMinFractionAP0*effectiveSigma+fgkParams.fPMinAP0);
580   return sigmaMMAP;
581 }
582 Double_t  AliESDv0::GetMinimaxSigmaD0(){
583   //
584   // calculate mini-max sigma of dca resolution
585   // 
586   //compv0->fTree->SetAlias("SigmaD2","max(min((SigmaD0+SigmaDE0)*0.5,1.5*SigmaDE0),0.5*SigmaDE0)");
587   Double_t    effectiveSigma = GetEffectiveSigmaD0();
588   Double_t    sigmaMMD0 = 0.5*(GetSigmaD0()+effectiveSigma);
589   sigmaMMD0  = TMath::Min(sigmaMMD0, fgkParams.fPMaxFractionD0*effectiveSigma);
590   sigmaMMD0  = TMath::Max(sigmaMMD0, fgkParams.fPMinFractionD0*effectiveSigma+fgkParams.fPMinD0);
591   return sigmaMMD0;
592 }
593
594
595 Double_t AliESDv0::GetLikelihoodAP(Int_t mode0, Int_t mode1){
596   //
597   // get likelihood for point angle
598   //
599   Double_t sigmaAP = 0.007;            //default sigma
600   switch (mode0){
601   case 0:
602     sigmaAP = GetSigmaAP0();           // mode 0  - covariance matrix estimates used 
603     break;
604   case 1:
605     sigmaAP = GetEffectiveSigmaAP0();  // mode 1 - effective sigma used
606     break;
607   case 2:
608     sigmaAP = GetMinimaxSigmaAP0();    // mode 2 - minimax sigma
609     break;
610   }
611   Double_t apNorm = TMath::Min(TMath::ACos(fPointAngle)/sigmaAP,50.);  
612   //normalized point angle, restricted - because of overflow problems in Exp
613   Double_t likelihood = 0;
614   switch(mode1){
615   case 0:
616     likelihood = TMath::Exp(-0.5*apNorm*apNorm);   
617     // one component
618     break;
619   case 1:
620     likelihood = (TMath::Exp(-0.5*apNorm*apNorm)+0.5* TMath::Exp(-0.25*apNorm*apNorm))/1.5;
621     // two components
622     break;
623   case 2:
624     likelihood = (TMath::Exp(-0.5*apNorm*apNorm)+0.5* TMath::Exp(-0.25*apNorm*apNorm)+0.25*TMath::Exp(-0.125*apNorm*apNorm))/1.75;
625     // three components
626     break;
627   }
628   return likelihood;
629 }
630
631 Double_t AliESDv0::GetLikelihoodD(Int_t mode0, Int_t mode1){
632   //
633   // get likelihood for DCA
634   //
635   Double_t sigmaD = 0.03;            //default sigma
636   switch (mode0){
637   case 0:
638     sigmaD = GetSigmaD0();           // mode 0  - covariance matrix estimates used 
639     break;
640   case 1:
641     sigmaD = GetEffectiveSigmaD0();  // mode 1 - effective sigma used
642     break;
643   case 2:
644     sigmaD = GetMinimaxSigmaD0();    // mode 2 - minimax sigma
645     break;
646   }
647
648   //Bo:  Double_t dNorm = TMath::Min(fDist2/sigmaD,50.);
649   Double_t dNorm = TMath::Min(fDcaV0Daughters/sigmaD,50.);//Bo:
650   //normalized point angle, restricted - because of overflow problems in Exp
651   Double_t likelihood = 0;
652   switch(mode1){
653   case 0:
654     likelihood = TMath::Exp(-2.*dNorm);   
655     // one component
656     break;
657   case 1:
658     likelihood = (TMath::Exp(-2.*dNorm)+0.5* TMath::Exp(-dNorm))/1.5;
659     // two components
660     break;
661   case 2:
662     likelihood = (TMath::Exp(-2.*dNorm)+0.5* TMath::Exp(-dNorm)+0.25*TMath::Exp(-0.5*dNorm))/1.75;
663     // three components
664     break;
665   }
666   return likelihood;
667
668 }
669
670 Double_t AliESDv0::GetLikelihoodC(Int_t mode0, Int_t /*mode1*/) const {
671   //
672   // get likelihood for Causality
673   // !!!  Causality variables defined in AliITStrackerMI !!! 
674   //      when more information was available
675   //  
676   Double_t likelihood = 0.5;
677   Double_t minCausal  = TMath::Min(fCausality[0],fCausality[1]);
678   Double_t maxCausal  = TMath::Max(fCausality[0],fCausality[1]);
679   //  minCausal           = TMath::Max(minCausal,0.5*maxCausal);
680   //compv0->fTree->SetAlias("LCausal","(1.05-(2*(0.8-exp(-max(RC.fV0rec.fCausality[0],RC.fV0rec.fCausality[1])))+2*(0.8-exp(-min(RC.fV0rec.fCausality[0],RC.fV0rec.fCausality[1]))))/2)**4");
681   
682   switch(mode0){
683   case 0:
684     //normalization 
685     likelihood = TMath::Power((1.05-2*(0.8-TMath::Exp(-maxCausal))),4.);
686     break;
687   case 1:
688     likelihood = TMath::Power(1.05-(2*(0.8-TMath::Exp(-maxCausal))+(2*(0.8-TMath::Exp(-minCausal))))*0.5,4.);
689     break;
690   }
691   return likelihood;
692   
693 }
694
695 void AliESDv0::SetCausality(Float_t pb0, Float_t pb1, Float_t pa0, Float_t pa1)
696 {
697   //
698   // set probabilities
699   //
700   fCausality[0] = pb0;     // probability - track 0 exist before vertex
701   fCausality[1] = pb1;     // probability - track 1 exist before vertex
702   fCausality[2] = pa0;     // probability - track 0 exist close after vertex
703   fCausality[3] = pa1;     // probability - track 1 exist close after vertex
704 }
705 void  AliESDv0::SetClusters(const Int_t *clp, const Int_t *clm)
706 {
707   //
708   // Set its clusters indexes
709   //
710   for (Int_t i=0;i<6;i++) fClusters[0][i] = clp[i]; 
711   for (Int_t i=0;i<6;i++) fClusters[1][i] = clm[i]; 
712 }
713
714 Double_t AliESDv0::GetEffMass(UInt_t p1, UInt_t p2) const{
715   //
716   // calculate effective mass
717   //
718   const Double_t kpmass[5] = {TDatabasePDG::Instance()->GetParticle(kElectron)->Mass(),
719                              TDatabasePDG::Instance()->GetParticle(kMuonMinus)->Mass(),
720                              TDatabasePDG::Instance()->GetParticle(kPiPlus)->Mass(),
721                              TDatabasePDG::Instance()->GetParticle(kKPlus)->Mass(),
722                              TDatabasePDG::Instance()->GetParticle(kProton)->Mass()};
723   /*
724   if (p1>4) return -1;
725   if (p2>4) return -1;
726   Double_t mass1 = kpmass[p1]; 
727   Double_t mass2 = kpmass[p2];   
728   const Double_t *m1 = fPmom;
729   const Double_t *m2 = fNmom;
730   //
731   //if (fRP[p1]+fRM[p2]<fRP[p2]+fRM[p1]){
732   //  m1 = fPM;
733   //  m2 = fPP;
734   //}
735   //
736   Double_t e1    = TMath::Sqrt(mass1*mass1+
737                               m1[0]*m1[0]+
738                               m1[1]*m1[1]+
739                               m1[2]*m1[2]);
740   Double_t e2    = TMath::Sqrt(mass2*mass2+
741                               m2[0]*m2[0]+
742                               m2[1]*m2[1]+
743                               m2[2]*m2[2]);  
744   Double_t mass =  
745     (m2[0]+m1[0])*(m2[0]+m1[0])+
746     (m2[1]+m1[1])*(m2[1]+m1[1])+
747     (m2[2]+m1[2])*(m2[2]+m1[2]);
748   
749   mass = (e1+e2)*(e1+e2)-mass;
750   if (mass < 0.) mass = 0.;
751   return (TMath::Sqrt(mass));
752   */
753   if(p1>4 || p2>4) return -1;
754   const AliExternalTrackParam *paramP = GetParamP();
755   const AliExternalTrackParam *paramN = GetParamN();
756   if (paramP->GetParameter()[4]<0){
757     paramP=GetParamN();
758     paramN=GetParamP();
759   }
760   Double_t pmom[3]={0}, nmom[3]={0};
761   paramP->GetPxPyPz(pmom);
762   paramN->GetPxPyPz(nmom);
763   Double_t e12   = kpmass[p1]*kpmass[p1]+pmom[0]*pmom[0]+pmom[1]*pmom[1]+pmom[2]*pmom[2];
764   Double_t e22   = kpmass[p2]*kpmass[p2]+nmom[0]*nmom[0]+nmom[1]*nmom[1]+nmom[2]*nmom[2];
765   Double_t cmass = TMath::Sqrt(TMath::Max(kpmass[p1]*kpmass[p1]+kpmass[p2]*kpmass[p2]
766                                           +2.*(TMath::Sqrt(e12*e22)-pmom[0]*nmom[0]-pmom[1]*nmom[1]-pmom[2]*nmom[2]),0.));
767   return cmass;
768                                
769 }
770
771
772
773 Double_t AliESDv0::GetKFInfo(UInt_t p1, UInt_t p2, Int_t type) const{
774   //
775   // type:
776   //   0 - return mass
777   //   1 - return err mass
778   //   2 - return chi2
779   // 
780   const Int_t spdg[5]={kPositron,kMuonPlus,kPiPlus, kKPlus, kProton};
781   const AliExternalTrackParam *paramP = GetParamP();
782   const AliExternalTrackParam *paramN = GetParamN();
783   if (paramP->GetSign()<0){
784     paramP=GetParamN();
785     paramN=GetParamP();
786   }
787   AliKFParticle kfp1(  *(paramP), spdg[p1] *TMath::Sign(1,p1) );
788   AliKFParticle kfp2( *(paramN), spdg[p2] *TMath::Sign(1,p2) );
789   AliKFParticle *v0KF = new AliKFParticle;
790   *(v0KF)+=kfp1;
791   *(v0KF)+=kfp2;
792   if (type==0) return v0KF->GetMass();
793   if (type==1) return v0KF->GetErrMass();
794   if (type==2) return v0KF->GetChi2();
795   return 0;
796 }
797
798
799 Double_t AliESDv0::GetKFInfoScale(UInt_t p1, UInt_t p2, Int_t type, Double_t d1pt, Double_t s1pt) const{
800   //
801   // type
802   //   0 - return mass
803   //   1 - return err mass
804   //   2 - return chi2
805   //   d1pt - 1/pt shift
806   //   s1pt - scaling of 1/pt
807   // Important function to benchmark the pt resolution, and to find out systematic distortion
808   //
809   const Int_t spdg[5]={kPositron,kMuonPlus,kPiPlus, kKPlus, kProton};
810   const AliExternalTrackParam *paramP = GetParamP();
811   const AliExternalTrackParam *paramN = GetParamN();
812   if (paramP->GetSign()<0){
813     paramP=GetParamP();
814     paramN=GetParamN();
815   }
816   Double_t *pparam1 = (Double_t*)paramP->GetParameter();
817   Double_t *pparam2 = (Double_t*)paramN->GetParameter();
818   pparam1[4]+=d1pt;
819   pparam2[4]+=d1pt;
820   pparam1[4]*=(1+s1pt);
821   pparam2[4]*=(1+s1pt);
822   //
823   AliKFParticle kfp1( *paramP, spdg[p1] *TMath::Sign(1,p1) );
824   AliKFParticle kfp2( *paramN, spdg[p2] *TMath::Sign(1,p2) );
825   AliKFParticle *v0KF = new AliKFParticle;
826   *(v0KF)+=kfp1;
827   *(v0KF)+=kfp2;
828   if (type==0) return v0KF->GetMass();
829   if (type==1) return v0KF->GetErrMass();
830   if (type==2) return v0KF->GetChi2();
831   return 0;
832 }