Implemented Copy() function for all esd objects to allow for assignment of AliESDEven...
[u/mrichter/AliRoot.git] / STEER / 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 <Riostream.h>
29 #include <TMath.h>
30 #include <TDatabasePDG.h>
31 #include <TPDGCode.h>
32 #include <TParticlePDG.h>
33
34 #include "AliLog.h"
35 #include "AliESDv0.h"
36 #include "AliExternalTrackParam.h"
37
38 ClassImp(AliESDv0)
39
40 const AliESDV0Params  AliESDv0::fgkParams;
41
42 AliESDv0::AliESDv0() :
43   TObject(),
44   fParamN(),
45   fParamP(),
46   fEffMass(TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass()),
47   fDcaV0Daughters(0),
48   fChi2V0(0.),
49   fRr(0),
50   fDistSigma(0),
51   fChi2Before(0),
52   fChi2After(0),
53   fPointAngleFi(0),
54   fPointAngleTh(0),
55   fPointAngle(0),
56   fPdgCode(kK0Short),
57   fNidx(0),
58   fPidx(0),
59   fStatus(0),
60   fNBefore(0),
61   fNAfter(0),
62   fOnFlyStatus(kFALSE)
63 {
64   //--------------------------------------------------------------------
65   // Default constructor  (K0s)
66   //--------------------------------------------------------------------
67
68   for (Int_t i=0; i<3; i++) {
69     fPos[i] = 0.;
70     fNmom[i] = 0.;
71     fPmom[i] = 0.;
72   }
73
74   for (Int_t i=0; i<6; i++) {
75     fPosCov[i]= 0.;
76   }
77
78   for (Int_t i=0;i<6;i++){fClusters[0][i]=0; fClusters[1][i]=0;}
79   fNormDCAPrim[0]=fNormDCAPrim[1]=0;
80   for (Int_t i=0;i<3;i++){fAngle[i]=0;}
81   for (Int_t i=0;i<4;i++){fCausality[i]=0;}
82 }
83
84 AliESDv0::AliESDv0(const AliESDv0& v0) :
85   TObject(v0),
86   fParamN(v0.fParamN),
87   fParamP(v0.fParamP),
88   fEffMass(v0.fEffMass),
89   fDcaV0Daughters(v0.fDcaV0Daughters),
90   fChi2V0(v0.fChi2V0),
91   fRr(v0.fRr),
92   fDistSigma(v0.fDistSigma),
93   fChi2Before(v0.fChi2Before),
94   fChi2After(v0.fChi2After),
95   fPointAngleFi(v0.fPointAngleFi),
96   fPointAngleTh(v0.fPointAngleTh),
97   fPointAngle(v0.fPointAngle),
98   fPdgCode(v0.fPdgCode),
99   fNidx(v0.fNidx),
100   fPidx(v0.fPidx),
101   fStatus(v0.fStatus),
102   fNBefore(v0.fNBefore),
103   fNAfter(v0.fNAfter),
104   fOnFlyStatus(v0.fOnFlyStatus)
105 {
106   //--------------------------------------------------------------------
107   // The copy constructor
108   //--------------------------------------------------------------------
109
110   for (int i=0; i<3; i++) {
111     fPos[i]  = v0.fPos[i];
112     fNmom[i] = v0.fNmom[i];
113     fPmom[i] = v0.fPmom[i];
114   }
115   for (int i=0; i<6; i++) {
116     fPosCov[i]  = v0.fPosCov[i];
117   }
118
119   for (Int_t i=0; i<2; i++) {
120     fNormDCAPrim[i]=v0.fNormDCAPrim[i];
121   }
122   for (Int_t i=0;i<6;i++){
123       fClusters[0][i]=v0.fClusters[0][i]; 
124       fClusters[1][i]=v0.fClusters[1][i];
125   }
126   for (Int_t i=0;i<3;i++){
127       fAngle[i]=v0.fAngle[i];
128   }
129   for (Int_t i=0;i<4;i++){fCausality[i]=v0.fCausality[i];}
130 }
131
132
133 AliESDv0::AliESDv0(const AliExternalTrackParam &t1, Int_t i1,
134                    const AliExternalTrackParam &t2, Int_t i2) :
135   TObject(),
136   fParamN(t1),
137   fParamP(t2),
138   fEffMass(TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass()),
139   fDcaV0Daughters(0),
140   fChi2V0(0.),
141   fRr(0),
142   fDistSigma(0),
143   fChi2Before(0),
144   fChi2After(0),
145   fPointAngleFi(0),
146   fPointAngleTh(0),
147   fPointAngle(0),
148   fPdgCode(kK0Short),
149   fNidx(i1),
150   fPidx(i2),
151   fStatus(0),
152   fNBefore(0),
153   fNAfter(0),
154   fOnFlyStatus(kFALSE)
155 {
156   //--------------------------------------------------------------------
157   // Main constructor  (K0s)
158   //--------------------------------------------------------------------
159
160   for (Int_t i=0; i<6; i++) {
161     fPosCov[i]= 0.;
162   }
163
164   //Trivial estimation of the vertex parameters
165   Double_t alpha=t1.GetAlpha(), cs=TMath::Cos(alpha), sn=TMath::Sin(alpha);
166   Double_t tmp[3];
167   t1.GetPxPyPz(tmp);
168   Double_t px1=tmp[0], py1=tmp[1], pz1=tmp[2];
169   t1.GetXYZ(tmp);
170   Double_t  x1=tmp[0],  y1=tmp[1],  z1=tmp[2];
171   const Double_t ss=0.0005*0.0005;//a kind of a residual misalignment precision
172   Double_t sx1=sn*sn*t1.GetSigmaY2()+ss, sy1=cs*cs*t1.GetSigmaY2()+ss; 
173
174
175   alpha=t2.GetAlpha(); cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
176   t2.GetPxPyPz(tmp);
177   Double_t px2=tmp[0], py2=tmp[1], pz2=tmp[2];
178   t2.GetXYZ(tmp);
179   Double_t  x2=tmp[0],  y2=tmp[1],  z2=tmp[2];
180   Double_t sx2=sn*sn*t2.GetSigmaY2()+ss, sy2=cs*cs*t2.GetSigmaY2()+ss; 
181     
182   Double_t sz1=t1.GetSigmaZ2(), sz2=t2.GetSigmaZ2();
183   Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
184   Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
185   Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
186   fPos[0]=wx1*x1 + wx2*x2; fPos[1]=wy1*y1 + wy2*y2; fPos[2]=wz1*z1 + wz2*z2;
187
188   //fPos[0]=0.5*(x1+x2); fPos[1]=0.5*(y1+y2); fPos[2]=0.5*(z1+z2);
189   fNmom[0]=px1; fNmom[1]=py1; fNmom[2]=pz1; 
190   fPmom[0]=px2; fPmom[1]=py2; fPmom[2]=pz2;
191
192   for (Int_t i=0;i<6;i++){fClusters[0][i]=0; fClusters[1][i]=0;}
193   fNormDCAPrim[0]=fNormDCAPrim[1]=0;
194   for (Int_t i=0;i<3;i++){fAngle[i]=0;}
195   for (Int_t i=0;i<4;i++){fCausality[i]=0;}
196 }
197
198 AliESDv0& AliESDv0::operator=(const AliESDv0 &v0){
199
200
201   //--------------------------------------------------------------------
202   // The assingment operator
203   //--------------------------------------------------------------------
204
205   if(this==&v0)return *this;
206   TObject::operator=(v0);
207   fParamN  = v0.fParamN;
208   fParamP  = v0.fParamP;
209   fEffMass = v0.fEffMass;
210   fDcaV0Daughters = v0.fDcaV0Daughters;
211   fChi2V0 = v0.fChi2V0;
212   fRr = v0.fRr;
213   fDistSigma    = v0.fDistSigma;
214   fChi2Before   = v0.fChi2Before;
215   fChi2After    = v0.fChi2After;
216   fPointAngleFi = v0.fPointAngleFi;
217   fPointAngleTh = v0.fPointAngleTh;
218   fPointAngle   = v0.fPointAngle;
219   fPdgCode      = v0.fPdgCode;
220   fNidx         = v0.fNidx;
221   fPidx         = v0.fPidx;
222   fStatus       = v0.fStatus;
223   fNBefore      = v0.fNBefore;
224   fNAfter       = v0.fNAfter;
225   fOnFlyStatus  = v0.fOnFlyStatus;
226
227   for (int i=0; i<3; i++) {
228     fPos[i]  = v0.fPos[i];
229     fNmom[i] = v0.fNmom[i];
230     fPmom[i] = v0.fPmom[i];
231   }
232   for (int i=0; i<6; i++) {
233     fPosCov[i]  = v0.fPosCov[i];
234   }
235   for (Int_t i=0; i<2; i++) {
236     fNormDCAPrim[i]=v0.fNormDCAPrim[i];
237   }
238   for (Int_t i=0;i<6;i++){
239       fClusters[0][i]=v0.fClusters[0][i]; 
240       fClusters[1][i]=v0.fClusters[1][i];
241   }
242   for (Int_t i=0;i<3;i++){
243       fAngle[i]=v0.fAngle[i];
244   }
245   for (Int_t i=0;i<4;i++){fCausality[i]=v0.fCausality[i];}
246
247   return *this;
248 }
249
250 void AliESDv0::Copy(TObject& obj) const {
251
252    // this overwrites the virtual TOBject::Copy()
253   // to allow run time copying without casting
254   // in AliESDEvent
255
256   if(this==&obj)return;
257   AliESDv0 *robj = dynamic_cast<AliESDv0*>(&obj);
258   if(!robj)return; // not an aliesesv0
259   *robj = *this;
260 }
261
262
263 AliESDv0::~AliESDv0(){
264   //--------------------------------------------------------------------
265   // Empty destructor
266   //--------------------------------------------------------------------
267 }
268
269
270
271 Double_t AliESDv0::ChangeMassHypothesis(Int_t code) {
272   //--------------------------------------------------------------------
273   // This function changes the mass hypothesis for this V0
274   // and returns the "kinematical quality" of this hypothesis 
275   //--------------------------------------------------------------------
276   static
277   Double_t piMass=TDatabasePDG::Instance()->GetParticle(kPiPlus)->Mass();
278   static
279   Double_t prMass=TDatabasePDG::Instance()->GetParticle(kProton)->Mass();
280   static
281   Double_t k0Mass=TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
282   static
283   Double_t l0Mass=TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
284
285   Double_t nmass=piMass, pmass=piMass, mass=k0Mass, ps=0.206;
286
287   fPdgCode=code;
288
289   switch (code) {
290   case kLambda0:
291     nmass=piMass; pmass=prMass; mass=l0Mass; ps=0.101; break;
292   case kLambda0Bar:
293     pmass=piMass; nmass=prMass; mass=l0Mass; ps=0.101; break;
294   case kK0Short: 
295     break;
296   default:
297     AliError("invalide PDG code ! Assuming K0s...");
298     fPdgCode=kK0Short;
299     break;
300   }
301
302   Double_t pxn=fNmom[0], pyn=fNmom[1], pzn=fNmom[2]; 
303   Double_t pxp=fPmom[0], pyp=fPmom[1], pzp=fPmom[2];
304
305   Double_t en=TMath::Sqrt(nmass*nmass + pxn*pxn + pyn*pyn + pzn*pzn);
306   Double_t ep=TMath::Sqrt(pmass*pmass + pxp*pxp + pyp*pyp + pzp*pzp);
307   Double_t pxl=pxn+pxp, pyl=pyn+pyp, pzl=pzn+pzp;
308   Double_t pl=TMath::Sqrt(pxl*pxl + pyl*pyl + pzl*pzl);
309
310   fEffMass=TMath::Sqrt((en+ep)*(en+ep)-pl*pl);
311
312   Double_t beta=pl/(en+ep);
313   Double_t pln=(pxn*pxl + pyn*pyl + pzn*pzl)/pl;
314   Double_t plp=(pxp*pxl + pyp*pyl + pzp*pzl)/pl;
315
316   Double_t pt2=pxp*pxp + pyp*pyp + pzp*pzp - plp*plp;
317
318   Double_t a=(plp-pln)/(plp+pln);
319   a -= (pmass*pmass-nmass*nmass)/(mass*mass);
320   a = 0.25*beta*beta*mass*mass*a*a + pt2;
321
322   return (a - ps*ps);
323   
324 }
325
326 void AliESDv0::GetPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const {
327   //--------------------------------------------------------------------
328   // This function returns V0's momentum (global)
329   //--------------------------------------------------------------------
330   px=fNmom[0]+fPmom[0]; 
331   py=fNmom[1]+fPmom[1]; 
332   pz=fNmom[2]+fPmom[2]; 
333 }
334
335 void AliESDv0::GetXYZ(Double_t &x, Double_t &y, Double_t &z) const {
336   //--------------------------------------------------------------------
337   // This function returns V0's position (global)
338   //--------------------------------------------------------------------
339   x=fPos[0]; 
340   y=fPos[1]; 
341   z=fPos[2]; 
342 }
343
344 Float_t AliESDv0::GetD(Double_t x0, Double_t y0, Double_t z0) const {
345   //--------------------------------------------------------------------
346   // This function returns V0's impact parameter
347   //--------------------------------------------------------------------
348   Double_t x=fPos[0],y=fPos[1],z=fPos[2];
349   Double_t px=fNmom[0]+fPmom[0];
350   Double_t py=fNmom[1]+fPmom[1];
351   Double_t pz=fNmom[2]+fPmom[2];
352
353   Double_t dx=(y0-y)*pz - (z0-z)*py; 
354   Double_t dy=(x0-x)*pz - (z0-z)*px;
355   Double_t dz=(x0-x)*py - (y0-y)*px;
356   Double_t d=TMath::Sqrt((dx*dx+dy*dy+dz*dz)/(px*px+py*py+pz*pz));
357   return d;
358 }
359
360
361 Float_t AliESDv0::GetV0CosineOfPointingAngle(Double_t& refPointX, Double_t& refPointY, Double_t& refPointZ) const {
362   // calculates the pointing angle of the V0 wrt a reference point
363
364   Double_t momV0[3]; //momentum of the V0
365   GetPxPyPz(momV0[0],momV0[1],momV0[2]);
366
367   Double_t deltaPos[3]; //vector between the reference point and the V0 vertex
368   deltaPos[0] = fPos[0] - refPointX;
369   deltaPos[1] = fPos[1] - refPointY;
370   deltaPos[2] = fPos[2] - refPointZ;
371
372   Double_t momV02    = momV0[0]*momV0[0] + momV0[1]*momV0[1] + momV0[2]*momV0[2];
373   Double_t deltaPos2 = deltaPos[0]*deltaPos[0] + deltaPos[1]*deltaPos[1] + deltaPos[2]*deltaPos[2];
374
375   Double_t cosinePointingAngle = (deltaPos[0]*momV0[0] +
376                                   deltaPos[1]*momV0[1] +
377                                   deltaPos[2]*momV0[2] ) /
378     TMath::Sqrt(momV02 * deltaPos2);
379   
380   return cosinePointingAngle;
381 }
382
383
384 // **** The following functions need to be revised
385
386 void AliESDv0::GetPosCov(Double_t cov[6]) const {
387
388   for (Int_t i=0; i<6; ++i) cov[i] = fPosCov[i];
389
390 }
391
392 Double_t AliESDv0::GetSigmaY(){
393   //
394   // return sigmay in y  at vertex position  using covariance matrix 
395   //
396   const Double_t * cp  = fParamP.GetCovariance();
397   const Double_t * cm  = fParamN.GetCovariance();
398   Double_t sigmay = cp[0]+cm[0]+ cp[5]*(fParamP.GetX()-fRr)*(fParamP.GetX()-fRr)+ cm[5]*(fParamN.GetX()-fRr)*(fParamN.GetX()-fRr);
399   return (sigmay>0) ? TMath::Sqrt(sigmay):100;
400 }
401
402 Double_t AliESDv0::GetSigmaZ(){
403   //
404   // return sigmay in y  at vertex position  using covariance matrix 
405   //
406   const Double_t * cp  = fParamP.GetCovariance();
407   const Double_t * cm  = fParamN.GetCovariance();
408   Double_t sigmaz = cp[2]+cm[2]+ cp[9]*(fParamP.GetX()-fRr)*(fParamP.GetX()-fRr)+ cm[9]*(fParamN.GetX()-fRr)*(fParamN.GetX()-fRr);
409   return (sigmaz>0) ? TMath::Sqrt(sigmaz):100;
410 }
411
412 Double_t AliESDv0::GetSigmaD0(){
413   //
414   // Sigma parameterization using covariance matrix
415   //
416   // sigma of distance between two tracks in vertex position 
417   // sigma of DCA is proportianal to sigmaD0
418   // factor 2 difference is explained by the fact that the DCA is calculated at the position 
419   // where the tracks as closest together ( not exact position of the vertex)
420   //
421   const Double_t * cp      = fParamP.GetCovariance();
422   const Double_t * cm      = fParamN.GetCovariance();
423   Double_t sigmaD0   = cp[0]+cm[0]+cp[2]+cm[2]+fgkParams.fPSigmaOffsetD0*fgkParams.fPSigmaOffsetD0;
424   sigmaD0           += ((fParamP.GetX()-fRr)*(fParamP.GetX()-fRr))*(cp[5]+cp[9]);
425   sigmaD0           += ((fParamN.GetX()-fRr)*(fParamN.GetX()-fRr))*(cm[5]+cm[9]);
426   return (sigmaD0>0)? TMath::Sqrt(sigmaD0):100;
427 }
428
429
430 Double_t AliESDv0::GetSigmaAP0(){
431   //
432   //Sigma parameterization using covariance matrices
433   //
434   Double_t prec  = TMath::Sqrt((fNmom[0]+fPmom[0])*(fNmom[0]+fPmom[0])
435                               +(fNmom[1]+fPmom[1])*(fNmom[1]+fPmom[1])
436                               +(fNmom[2]+fPmom[2])*(fNmom[2]+fPmom[2]));
437   Double_t normp = TMath::Sqrt(fPmom[0]*fPmom[0]+fPmom[1]*fPmom[1]+fPmom[2]*fPmom[2])/prec;  // fraction of the momenta
438   Double_t normm = TMath::Sqrt(fNmom[0]*fNmom[0]+fNmom[1]*fNmom[1]+fNmom[2]*fNmom[2])/prec;  
439   const Double_t * cp      = fParamP.GetCovariance();
440   const Double_t * cm      = fParamN.GetCovariance();
441   Double_t sigmaAP0 = fgkParams.fPSigmaOffsetAP0*fgkParams.fPSigmaOffsetAP0;                           // minimal part
442   sigmaAP0 +=  (cp[5]+cp[9])*(normp*normp)+(cm[5]+cm[9])*(normm*normm);          // angular resolution part
443   Double_t sigmaAP1 = GetSigmaD0()/(TMath::Abs(fRr)+0.01);                       // vertex position part
444   sigmaAP0 +=  0.5*sigmaAP1*sigmaAP1;                              
445   return (sigmaAP0>0)? TMath::Sqrt(sigmaAP0):100;
446 }
447
448 Double_t AliESDv0::GetEffectiveSigmaD0(){
449   //
450   // minimax - effective Sigma parameterization 
451   // p12 effective curvature and v0 radius postion used as parameters  
452   //  
453   Double_t p12 = TMath::Sqrt(fParamP.GetParameter()[4]*fParamP.GetParameter()[4]+
454                              fParamN.GetParameter()[4]*fParamN.GetParameter()[4]);
455   Double_t sigmaED0= TMath::Max(TMath::Sqrt(fRr)-fgkParams.fPSigmaRminDE,0.0)*fgkParams.fPSigmaCoefDE*p12*p12;
456   sigmaED0*= sigmaED0;
457   sigmaED0*= sigmaED0;
458   sigmaED0 = TMath::Sqrt(sigmaED0+fgkParams.fPSigmaOffsetDE*fgkParams.fPSigmaOffsetDE);
459   return (sigmaED0<fgkParams.fPSigmaMaxDE) ? sigmaED0: fgkParams.fPSigmaMaxDE;
460 }
461
462
463 Double_t AliESDv0::GetEffectiveSigmaAP0(){
464   //
465   // effective Sigma parameterization of point angle resolution 
466   //
467   Double_t p12 = TMath::Sqrt(fParamP.GetParameter()[4]*fParamP.GetParameter()[4]+
468                              fParamN.GetParameter()[4]*fParamN.GetParameter()[4]);
469   Double_t sigmaAPE= fgkParams.fPSigmaBase0APE;
470   sigmaAPE+= fgkParams.fPSigmaR0APE/(fgkParams.fPSigmaR1APE+fRr);
471   sigmaAPE*= (fgkParams.fPSigmaP0APE+fgkParams.fPSigmaP1APE*p12);
472   sigmaAPE = TMath::Min(sigmaAPE,fgkParams.fPSigmaMaxAPE);
473   return sigmaAPE;
474 }
475
476
477 Double_t  AliESDv0::GetMinimaxSigmaAP0(){
478   //
479   // calculate mini-max effective sigma of point angle resolution
480   //
481   //compv0->fTree->SetAlias("SigmaAP2","max(min((SigmaAP0+SigmaAPE0)*0.5,1.5*SigmaAPE0),0.5*SigmaAPE0+0.003)");
482   Double_t    effectiveSigma = GetEffectiveSigmaAP0();
483   Double_t    sigmaMMAP = 0.5*(GetSigmaAP0()+effectiveSigma);
484   sigmaMMAP  = TMath::Min(sigmaMMAP, fgkParams.fPMaxFractionAP0*effectiveSigma);
485   sigmaMMAP  = TMath::Max(sigmaMMAP, fgkParams.fPMinFractionAP0*effectiveSigma+fgkParams.fPMinAP0);
486   return sigmaMMAP;
487 }
488 Double_t  AliESDv0::GetMinimaxSigmaD0(){
489   //
490   // calculate mini-max sigma of dca resolution
491   // 
492   //compv0->fTree->SetAlias("SigmaD2","max(min((SigmaD0+SigmaDE0)*0.5,1.5*SigmaDE0),0.5*SigmaDE0)");
493   Double_t    effectiveSigma = GetEffectiveSigmaD0();
494   Double_t    sigmaMMD0 = 0.5*(GetSigmaD0()+effectiveSigma);
495   sigmaMMD0  = TMath::Min(sigmaMMD0, fgkParams.fPMaxFractionD0*effectiveSigma);
496   sigmaMMD0  = TMath::Max(sigmaMMD0, fgkParams.fPMinFractionD0*effectiveSigma+fgkParams.fPMinD0);
497   return sigmaMMD0;
498 }
499
500
501 Double_t AliESDv0::GetLikelihoodAP(Int_t mode0, Int_t mode1){
502   //
503   // get likelihood for point angle
504   //
505   Double_t sigmaAP = 0.007;            //default sigma
506   switch (mode0){
507   case 0:
508     sigmaAP = GetSigmaAP0();           // mode 0  - covariance matrix estimates used 
509     break;
510   case 1:
511     sigmaAP = GetEffectiveSigmaAP0();  // mode 1 - effective sigma used
512     break;
513   case 2:
514     sigmaAP = GetMinimaxSigmaAP0();    // mode 2 - minimax sigma
515     break;
516   }
517   Double_t apNorm = TMath::Min(TMath::ACos(fPointAngle)/sigmaAP,50.);  
518   //normalized point angle, restricted - because of overflow problems in Exp
519   Double_t likelihood = 0;
520   switch(mode1){
521   case 0:
522     likelihood = TMath::Exp(-0.5*apNorm*apNorm);   
523     // one component
524     break;
525   case 1:
526     likelihood = (TMath::Exp(-0.5*apNorm*apNorm)+0.5* TMath::Exp(-0.25*apNorm*apNorm))/1.5;
527     // two components
528     break;
529   case 2:
530     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;
531     // three components
532     break;
533   }
534   return likelihood;
535 }
536
537 Double_t AliESDv0::GetLikelihoodD(Int_t mode0, Int_t mode1){
538   //
539   // get likelihood for DCA
540   //
541   Double_t sigmaD = 0.03;            //default sigma
542   switch (mode0){
543   case 0:
544     sigmaD = GetSigmaD0();           // mode 0  - covariance matrix estimates used 
545     break;
546   case 1:
547     sigmaD = GetEffectiveSigmaD0();  // mode 1 - effective sigma used
548     break;
549   case 2:
550     sigmaD = GetMinimaxSigmaD0();    // mode 2 - minimax sigma
551     break;
552   }
553
554   //Bo:  Double_t dNorm = TMath::Min(fDist2/sigmaD,50.);
555   Double_t dNorm = TMath::Min(fDcaV0Daughters/sigmaD,50.);//Bo:
556   //normalized point angle, restricted - because of overflow problems in Exp
557   Double_t likelihood = 0;
558   switch(mode1){
559   case 0:
560     likelihood = TMath::Exp(-2.*dNorm);   
561     // one component
562     break;
563   case 1:
564     likelihood = (TMath::Exp(-2.*dNorm)+0.5* TMath::Exp(-dNorm))/1.5;
565     // two components
566     break;
567   case 2:
568     likelihood = (TMath::Exp(-2.*dNorm)+0.5* TMath::Exp(-dNorm)+0.25*TMath::Exp(-0.5*dNorm))/1.75;
569     // three components
570     break;
571   }
572   return likelihood;
573
574 }
575
576 Double_t AliESDv0::GetLikelihoodC(Int_t mode0, Int_t /*mode1*/){
577   //
578   // get likelihood for Causality
579   // !!!  Causality variables defined in AliITStrackerMI !!! 
580   //      when more information was available
581   //  
582   Double_t likelihood = 0.5;
583   Double_t minCausal  = TMath::Min(fCausality[0],fCausality[1]);
584   Double_t maxCausal  = TMath::Max(fCausality[0],fCausality[1]);
585   //  minCausal           = TMath::Max(minCausal,0.5*maxCausal);
586   //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");
587   
588   switch(mode0){
589   case 0:
590     //normalization 
591     likelihood = TMath::Power((1.05-2*(0.8-TMath::Exp(-maxCausal))),4.);
592     break;
593   case 1:
594     likelihood = TMath::Power(1.05-(2*(0.8-TMath::Exp(-maxCausal))+(2*(0.8-TMath::Exp(-minCausal))))*0.5,4.);
595     break;
596   }
597   return likelihood;
598   
599 }
600
601 void AliESDv0::SetCausality(Float_t pb0, Float_t pb1, Float_t pa0, Float_t pa1)
602 {
603   //
604   // set probabilities
605   //
606   fCausality[0] = pb0;     // probability - track 0 exist before vertex
607   fCausality[1] = pb1;     // probability - track 1 exist before vertex
608   fCausality[2] = pa0;     // probability - track 0 exist close after vertex
609   fCausality[3] = pa1;     // probability - track 1 exist close after vertex
610 }
611 void  AliESDv0::SetClusters(Int_t *clp, Int_t *clm)
612 {
613   //
614   // Set its clusters indexes
615   //
616   for (Int_t i=0;i<6;i++) fClusters[0][i] = clp[i]; 
617   for (Int_t i=0;i<6;i++) fClusters[1][i] = clm[i]; 
618 }
619
620 Float_t AliESDv0::GetEffMass(UInt_t p1, UInt_t p2){
621   //
622   // calculate effective mass
623   //
624   const Float_t kpmass[5] = {5.10000000000000037e-04,1.05660000000000004e-01,1.39570000000000000e-01,
625                       4.93599999999999983e-01, 9.38270000000000048e-01};
626   if (p1>4) return -1;
627   if (p2>4) return -1;
628   Float_t mass1 = kpmass[p1]; 
629   Float_t mass2 = kpmass[p2];   
630   Double_t *m1 = fPmom;
631   Double_t *m2 = fNmom;
632   //
633   //if (fRP[p1]+fRM[p2]<fRP[p2]+fRM[p1]){
634   //  m1 = fPM;
635   //  m2 = fPP;
636   //}
637   //
638   Float_t e1    = TMath::Sqrt(mass1*mass1+
639                               m1[0]*m1[0]+
640                               m1[1]*m1[1]+
641                               m1[2]*m1[2]);
642   Float_t e2    = TMath::Sqrt(mass2*mass2+
643                               m2[0]*m2[0]+
644                               m2[1]*m2[1]+
645                               m2[2]*m2[2]);  
646   Float_t mass =  
647     (m2[0]+m1[0])*(m2[0]+m1[0])+
648     (m2[1]+m1[1])*(m2[1]+m1[1])+
649     (m2[2]+m1[2])*(m2[2]+m1[2]);
650   
651   mass = TMath::Sqrt((e1+e2)*(e1+e2)-mass);
652   return mass;
653 }