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