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