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