Change from Int_t/SHort_t to smaller data types, the change from Float_t/Double_t...
[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 AliESDv0::AliESDv0(const AliExternalTrackParam &t1, Int_t i1,
133                    const AliExternalTrackParam &t2, Int_t i2) :
134   TObject(),
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(){
198   //--------------------------------------------------------------------
199   // Empty destructor
200   //--------------------------------------------------------------------
201 }
202
203
204
205 Double_t AliESDv0::ChangeMassHypothesis(Int_t code) {
206   //--------------------------------------------------------------------
207   // This function changes the mass hypothesis for this V0
208   // and returns the "kinematical quality" of this hypothesis 
209   //--------------------------------------------------------------------
210   static
211   Double_t piMass=TDatabasePDG::Instance()->GetParticle(kPiPlus)->Mass();
212   static
213   Double_t prMass=TDatabasePDG::Instance()->GetParticle(kProton)->Mass();
214   static
215   Double_t k0Mass=TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
216   static
217   Double_t l0Mass=TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
218
219   Double_t nmass=piMass, pmass=piMass, mass=k0Mass, ps=0.206;
220
221   fPdgCode=code;
222
223   switch (code) {
224   case kLambda0:
225     nmass=piMass; pmass=prMass; mass=l0Mass; ps=0.101; break;
226   case kLambda0Bar:
227     pmass=piMass; nmass=prMass; mass=l0Mass; ps=0.101; break;
228   case kK0Short: 
229     break;
230   default:
231     AliError("invalide PDG code ! Assuming K0s...");
232     fPdgCode=kK0Short;
233     break;
234   }
235
236   Double_t pxn=fNmom[0], pyn=fNmom[1], pzn=fNmom[2]; 
237   Double_t pxp=fPmom[0], pyp=fPmom[1], pzp=fPmom[2];
238
239   Double_t en=TMath::Sqrt(nmass*nmass + pxn*pxn + pyn*pyn + pzn*pzn);
240   Double_t ep=TMath::Sqrt(pmass*pmass + pxp*pxp + pyp*pyp + pzp*pzp);
241   Double_t pxl=pxn+pxp, pyl=pyn+pyp, pzl=pzn+pzp;
242   Double_t pl=TMath::Sqrt(pxl*pxl + pyl*pyl + pzl*pzl);
243
244   fEffMass=TMath::Sqrt((en+ep)*(en+ep)-pl*pl);
245
246   Double_t beta=pl/(en+ep);
247   Double_t pln=(pxn*pxl + pyn*pyl + pzn*pzl)/pl;
248   Double_t plp=(pxp*pxl + pyp*pyl + pzp*pzl)/pl;
249
250   Double_t pt2=pxp*pxp + pyp*pyp + pzp*pzp - plp*plp;
251
252   Double_t a=(plp-pln)/(plp+pln);
253   a -= (pmass*pmass-nmass*nmass)/(mass*mass);
254   a = 0.25*beta*beta*mass*mass*a*a + pt2;
255
256   return (a - ps*ps);
257   
258 }
259
260 void AliESDv0::GetPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const {
261   //--------------------------------------------------------------------
262   // This function returns V0's momentum (global)
263   //--------------------------------------------------------------------
264   px=fNmom[0]+fPmom[0]; 
265   py=fNmom[1]+fPmom[1]; 
266   pz=fNmom[2]+fPmom[2]; 
267 }
268
269 void AliESDv0::GetXYZ(Double_t &x, Double_t &y, Double_t &z) const {
270   //--------------------------------------------------------------------
271   // This function returns V0's position (global)
272   //--------------------------------------------------------------------
273   x=fPos[0]; 
274   y=fPos[1]; 
275   z=fPos[2]; 
276 }
277
278 Float_t AliESDv0::GetD(Double_t x0, Double_t y0, Double_t z0) const {
279   //--------------------------------------------------------------------
280   // This function returns V0's impact parameter
281   //--------------------------------------------------------------------
282   Double_t x=fPos[0],y=fPos[1],z=fPos[2];
283   Double_t px=fNmom[0]+fPmom[0];
284   Double_t py=fNmom[1]+fPmom[1];
285   Double_t pz=fNmom[2]+fPmom[2];
286
287   Double_t dx=(y0-y)*pz - (z0-z)*py; 
288   Double_t dy=(x0-x)*pz - (z0-z)*px;
289   Double_t dz=(x0-x)*py - (y0-y)*px;
290   Double_t d=TMath::Sqrt((dx*dx+dy*dy+dz*dz)/(px*px+py*py+pz*pz));
291   return d;
292 }
293
294
295 Float_t AliESDv0::GetV0CosineOfPointingAngle(Double_t& refPointX, Double_t& refPointY, Double_t& refPointZ) const {
296   // calculates the pointing angle of the V0 wrt a reference point
297
298   Double_t momV0[3]; //momentum of the V0
299   GetPxPyPz(momV0[0],momV0[1],momV0[2]);
300
301   Double_t deltaPos[3]; //vector between the reference point and the V0 vertex
302   deltaPos[0] = fPos[0] - refPointX;
303   deltaPos[1] = fPos[1] - refPointY;
304   deltaPos[2] = fPos[2] - refPointZ;
305
306   Double_t momV02    = momV0[0]*momV0[0] + momV0[1]*momV0[1] + momV0[2]*momV0[2];
307   Double_t deltaPos2 = deltaPos[0]*deltaPos[0] + deltaPos[1]*deltaPos[1] + deltaPos[2]*deltaPos[2];
308
309   Double_t cosinePointingAngle = (deltaPos[0]*momV0[0] +
310                                   deltaPos[1]*momV0[1] +
311                                   deltaPos[2]*momV0[2] ) /
312     TMath::Sqrt(momV02 * deltaPos2);
313   
314   return cosinePointingAngle;
315 }
316
317
318 // **** The following functions need to be revised
319
320 void AliESDv0::GetPosCov(Double_t cov[6]) const {
321
322   for (Int_t i=0; i<6; ++i) cov[i] = fPosCov[i];
323
324 }
325
326 Double_t AliESDv0::GetSigmaY(){
327   //
328   // return sigmay in y  at vertex position  using covariance matrix 
329   //
330   const Double_t * cp  = fParamP.GetCovariance();
331   const Double_t * cm  = fParamN.GetCovariance();
332   Double_t sigmay = cp[0]+cm[0]+ cp[5]*(fParamP.GetX()-fRr)*(fParamP.GetX()-fRr)+ cm[5]*(fParamN.GetX()-fRr)*(fParamN.GetX()-fRr);
333   return (sigmay>0) ? TMath::Sqrt(sigmay):100;
334 }
335
336 Double_t AliESDv0::GetSigmaZ(){
337   //
338   // return sigmay in y  at vertex position  using covariance matrix 
339   //
340   const Double_t * cp  = fParamP.GetCovariance();
341   const Double_t * cm  = fParamN.GetCovariance();
342   Double_t sigmaz = cp[2]+cm[2]+ cp[9]*(fParamP.GetX()-fRr)*(fParamP.GetX()-fRr)+ cm[9]*(fParamN.GetX()-fRr)*(fParamN.GetX()-fRr);
343   return (sigmaz>0) ? TMath::Sqrt(sigmaz):100;
344 }
345
346 Double_t AliESDv0::GetSigmaD0(){
347   //
348   // Sigma parameterization using covariance matrix
349   //
350   // sigma of distance between two tracks in vertex position 
351   // sigma of DCA is proportianal to sigmaD0
352   // factor 2 difference is explained by the fact that the DCA is calculated at the position 
353   // where the tracks as closest together ( not exact position of the vertex)
354   //
355   const Double_t * cp      = fParamP.GetCovariance();
356   const Double_t * cm      = fParamN.GetCovariance();
357   Double_t sigmaD0   = cp[0]+cm[0]+cp[2]+cm[2]+fgkParams.fPSigmaOffsetD0*fgkParams.fPSigmaOffsetD0;
358   sigmaD0           += ((fParamP.GetX()-fRr)*(fParamP.GetX()-fRr))*(cp[5]+cp[9]);
359   sigmaD0           += ((fParamN.GetX()-fRr)*(fParamN.GetX()-fRr))*(cm[5]+cm[9]);
360   return (sigmaD0>0)? TMath::Sqrt(sigmaD0):100;
361 }
362
363
364 Double_t AliESDv0::GetSigmaAP0(){
365   //
366   //Sigma parameterization using covariance matrices
367   //
368   Double_t prec  = TMath::Sqrt((fNmom[0]+fPmom[0])*(fNmom[0]+fPmom[0])
369                               +(fNmom[1]+fPmom[1])*(fNmom[1]+fPmom[1])
370                               +(fNmom[2]+fPmom[2])*(fNmom[2]+fPmom[2]));
371   Double_t normp = TMath::Sqrt(fPmom[0]*fPmom[0]+fPmom[1]*fPmom[1]+fPmom[2]*fPmom[2])/prec;  // fraction of the momenta
372   Double_t normm = TMath::Sqrt(fNmom[0]*fNmom[0]+fNmom[1]*fNmom[1]+fNmom[2]*fNmom[2])/prec;  
373   const Double_t * cp      = fParamP.GetCovariance();
374   const Double_t * cm      = fParamN.GetCovariance();
375   Double_t sigmaAP0 = fgkParams.fPSigmaOffsetAP0*fgkParams.fPSigmaOffsetAP0;                           // minimal part
376   sigmaAP0 +=  (cp[5]+cp[9])*(normp*normp)+(cm[5]+cm[9])*(normm*normm);          // angular resolution part
377   Double_t sigmaAP1 = GetSigmaD0()/(TMath::Abs(fRr)+0.01);                       // vertex position part
378   sigmaAP0 +=  0.5*sigmaAP1*sigmaAP1;                              
379   return (sigmaAP0>0)? TMath::Sqrt(sigmaAP0):100;
380 }
381
382 Double_t AliESDv0::GetEffectiveSigmaD0(){
383   //
384   // minimax - effective Sigma parameterization 
385   // p12 effective curvature and v0 radius postion used as parameters  
386   //  
387   Double_t p12 = TMath::Sqrt(fParamP.GetParameter()[4]*fParamP.GetParameter()[4]+
388                              fParamN.GetParameter()[4]*fParamN.GetParameter()[4]);
389   Double_t sigmaED0= TMath::Max(TMath::Sqrt(fRr)-fgkParams.fPSigmaRminDE,0.0)*fgkParams.fPSigmaCoefDE*p12*p12;
390   sigmaED0*= sigmaED0;
391   sigmaED0*= sigmaED0;
392   sigmaED0 = TMath::Sqrt(sigmaED0+fgkParams.fPSigmaOffsetDE*fgkParams.fPSigmaOffsetDE);
393   return (sigmaED0<fgkParams.fPSigmaMaxDE) ? sigmaED0: fgkParams.fPSigmaMaxDE;
394 }
395
396
397 Double_t AliESDv0::GetEffectiveSigmaAP0(){
398   //
399   // effective Sigma parameterization of point angle resolution 
400   //
401   Double_t p12 = TMath::Sqrt(fParamP.GetParameter()[4]*fParamP.GetParameter()[4]+
402                              fParamN.GetParameter()[4]*fParamN.GetParameter()[4]);
403   Double_t sigmaAPE= fgkParams.fPSigmaBase0APE;
404   sigmaAPE+= fgkParams.fPSigmaR0APE/(fgkParams.fPSigmaR1APE+fRr);
405   sigmaAPE*= (fgkParams.fPSigmaP0APE+fgkParams.fPSigmaP1APE*p12);
406   sigmaAPE = TMath::Min(sigmaAPE,fgkParams.fPSigmaMaxAPE);
407   return sigmaAPE;
408 }
409
410
411 Double_t  AliESDv0::GetMinimaxSigmaAP0(){
412   //
413   // calculate mini-max effective sigma of point angle resolution
414   //
415   //compv0->fTree->SetAlias("SigmaAP2","max(min((SigmaAP0+SigmaAPE0)*0.5,1.5*SigmaAPE0),0.5*SigmaAPE0+0.003)");
416   Double_t    effectiveSigma = GetEffectiveSigmaAP0();
417   Double_t    sigmaMMAP = 0.5*(GetSigmaAP0()+effectiveSigma);
418   sigmaMMAP  = TMath::Min(sigmaMMAP, fgkParams.fPMaxFractionAP0*effectiveSigma);
419   sigmaMMAP  = TMath::Max(sigmaMMAP, fgkParams.fPMinFractionAP0*effectiveSigma+fgkParams.fPMinAP0);
420   return sigmaMMAP;
421 }
422 Double_t  AliESDv0::GetMinimaxSigmaD0(){
423   //
424   // calculate mini-max sigma of dca resolution
425   // 
426   //compv0->fTree->SetAlias("SigmaD2","max(min((SigmaD0+SigmaDE0)*0.5,1.5*SigmaDE0),0.5*SigmaDE0)");
427   Double_t    effectiveSigma = GetEffectiveSigmaD0();
428   Double_t    sigmaMMD0 = 0.5*(GetSigmaD0()+effectiveSigma);
429   sigmaMMD0  = TMath::Min(sigmaMMD0, fgkParams.fPMaxFractionD0*effectiveSigma);
430   sigmaMMD0  = TMath::Max(sigmaMMD0, fgkParams.fPMinFractionD0*effectiveSigma+fgkParams.fPMinD0);
431   return sigmaMMD0;
432 }
433
434
435 Double_t AliESDv0::GetLikelihoodAP(Int_t mode0, Int_t mode1){
436   //
437   // get likelihood for point angle
438   //
439   Double_t sigmaAP = 0.007;            //default sigma
440   switch (mode0){
441   case 0:
442     sigmaAP = GetSigmaAP0();           // mode 0  - covariance matrix estimates used 
443     break;
444   case 1:
445     sigmaAP = GetEffectiveSigmaAP0();  // mode 1 - effective sigma used
446     break;
447   case 2:
448     sigmaAP = GetMinimaxSigmaAP0();    // mode 2 - minimax sigma
449     break;
450   }
451   Double_t apNorm = TMath::Min(TMath::ACos(fPointAngle)/sigmaAP,50.);  
452   //normalized point angle, restricted - because of overflow problems in Exp
453   Double_t likelihood = 0;
454   switch(mode1){
455   case 0:
456     likelihood = TMath::Exp(-0.5*apNorm*apNorm);   
457     // one component
458     break;
459   case 1:
460     likelihood = (TMath::Exp(-0.5*apNorm*apNorm)+0.5* TMath::Exp(-0.25*apNorm*apNorm))/1.5;
461     // two components
462     break;
463   case 2:
464     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;
465     // three components
466     break;
467   }
468   return likelihood;
469 }
470
471 Double_t AliESDv0::GetLikelihoodD(Int_t mode0, Int_t mode1){
472   //
473   // get likelihood for DCA
474   //
475   Double_t sigmaD = 0.03;            //default sigma
476   switch (mode0){
477   case 0:
478     sigmaD = GetSigmaD0();           // mode 0  - covariance matrix estimates used 
479     break;
480   case 1:
481     sigmaD = GetEffectiveSigmaD0();  // mode 1 - effective sigma used
482     break;
483   case 2:
484     sigmaD = GetMinimaxSigmaD0();    // mode 2 - minimax sigma
485     break;
486   }
487
488   //Bo:  Double_t dNorm = TMath::Min(fDist2/sigmaD,50.);
489   Double_t dNorm = TMath::Min(fDcaV0Daughters/sigmaD,50.);//Bo:
490   //normalized point angle, restricted - because of overflow problems in Exp
491   Double_t likelihood = 0;
492   switch(mode1){
493   case 0:
494     likelihood = TMath::Exp(-2.*dNorm);   
495     // one component
496     break;
497   case 1:
498     likelihood = (TMath::Exp(-2.*dNorm)+0.5* TMath::Exp(-dNorm))/1.5;
499     // two components
500     break;
501   case 2:
502     likelihood = (TMath::Exp(-2.*dNorm)+0.5* TMath::Exp(-dNorm)+0.25*TMath::Exp(-0.5*dNorm))/1.75;
503     // three components
504     break;
505   }
506   return likelihood;
507
508 }
509
510 Double_t AliESDv0::GetLikelihoodC(Int_t mode0, Int_t /*mode1*/){
511   //
512   // get likelihood for Causality
513   // !!!  Causality variables defined in AliITStrackerMI !!! 
514   //      when more information was available
515   //  
516   Double_t likelihood = 0.5;
517   Double_t minCausal  = TMath::Min(fCausality[0],fCausality[1]);
518   Double_t maxCausal  = TMath::Max(fCausality[0],fCausality[1]);
519   //  minCausal           = TMath::Max(minCausal,0.5*maxCausal);
520   //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");
521   
522   switch(mode0){
523   case 0:
524     //normalization 
525     likelihood = TMath::Power((1.05-2*(0.8-TMath::Exp(-maxCausal))),4.);
526     break;
527   case 1:
528     likelihood = TMath::Power(1.05-(2*(0.8-TMath::Exp(-maxCausal))+(2*(0.8-TMath::Exp(-minCausal))))*0.5,4.);
529     break;
530   }
531   return likelihood;
532   
533 }
534
535 void AliESDv0::SetCausality(Float_t pb0, Float_t pb1, Float_t pa0, Float_t pa1)
536 {
537   //
538   // set probabilities
539   //
540   fCausality[0] = pb0;     // probability - track 0 exist before vertex
541   fCausality[1] = pb1;     // probability - track 1 exist before vertex
542   fCausality[2] = pa0;     // probability - track 0 exist close after vertex
543   fCausality[3] = pa1;     // probability - track 1 exist close after vertex
544 }
545 void  AliESDv0::SetClusters(Int_t *clp, Int_t *clm)
546 {
547   //
548   // Set its clusters indexes
549   //
550   for (Int_t i=0;i<6;i++) fClusters[0][i] = clp[i]; 
551   for (Int_t i=0;i<6;i++) fClusters[1][i] = clm[i]; 
552 }
553
554 Float_t AliESDv0::GetEffMass(UInt_t p1, UInt_t p2){
555   //
556   // calculate effective mass
557   //
558   const Float_t kpmass[5] = {5.10000000000000037e-04,1.05660000000000004e-01,1.39570000000000000e-01,
559                       4.93599999999999983e-01, 9.38270000000000048e-01};
560   if (p1>4) return -1;
561   if (p2>4) return -1;
562   Float_t mass1 = kpmass[p1]; 
563   Float_t mass2 = kpmass[p2];   
564   Double_t *m1 = fPmom;
565   Double_t *m2 = fNmom;
566   //
567   //if (fRP[p1]+fRM[p2]<fRP[p2]+fRM[p1]){
568   //  m1 = fPM;
569   //  m2 = fPP;
570   //}
571   //
572   Float_t e1    = TMath::Sqrt(mass1*mass1+
573                               m1[0]*m1[0]+
574                               m1[1]*m1[1]+
575                               m1[2]*m1[2]);
576   Float_t e2    = TMath::Sqrt(mass2*mass2+
577                               m2[0]*m2[0]+
578                               m2[1]*m2[1]+
579                               m2[2]*m2[2]);  
580   Float_t mass =  
581     (m2[0]+m1[0])*(m2[0]+m1[0])+
582     (m2[1]+m1[1])*(m2[1]+m1[1])+
583     (m2[2]+m1[2])*(m2[2]+m1[2]);
584   
585   mass = TMath::Sqrt((e1+e2)*(e1+e2)-mass);
586   return mass;
587 }