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