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 }