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