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