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