STEER/ESD/AliESDv0.cxx

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