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   //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, Double_t z0) const {
 421   //--------------------------------------------------------------------
 422   // This function returns V0's impact parameter
 423   //--------------------------------------------------------------------
 424   Double_t x=fPos[0],y=fPos[1],z=fPos[2];
 425   Double_t px=fNmom[0]+fPmom[0];
 426   Double_t py=fNmom[1]+fPmom[1];
 427   Double_t pz=fNmom[2]+fPmom[2];
 428
 429   Double_t dx=(y0-y)*pz - (z0-z)*py;
 430   Double_t dy=(x0-x)*pz - (z0-z)*px;
 431   Double_t dz=(x0-x)*py - (y0-y)*px;
 432   Double_t d=TMath::Sqrt((dx*dx+dy*dy+dz*dz)/(px*px+py*py+pz*pz));
 433   return d;
 434 }
 435
 436 Float_t AliESDv0::GetV0CosineOfPointingAngle(Double_t refPointX, Double_t refPointY, Double_t refPointZ) const {
 437   // calculates the pointing angle of the V0 wrt a reference point
 438
 439   Double_t momV0[3]; //momentum of the V0
 440   GetPxPyPz(momV0[0],momV0[1],momV0[2]);
 441
 442   Double_t deltaPos[3]; //vector between the reference point and the V0 vertex
 443   deltaPos[0] = fPos[0] - refPointX;
 444   deltaPos[1] = fPos[1] - refPointY;
 445   deltaPos[2] = fPos[2] - refPointZ;
 446
 447   Double_t momV02    = momV0[0]*momV0[0] + momV0[1]*momV0[1] + momV0[2]*momV0[2];
 448   Double_t deltaPos2 = deltaPos[0]*deltaPos[0] + deltaPos[1]*deltaPos[1] + deltaPos[2]*deltaPos[2];
 449
 450   Double_t cosinePointingAngle = (deltaPos[0]*momV0[0] +
 451                                   deltaPos[1]*momV0[1] +
 452                                   deltaPos[2]*momV0[2] ) /
 453     TMath::Sqrt(momV02 * deltaPos2);
 454
 455   return cosinePointingAngle;
 456 }
 457
 458
 459 // **** The following functions need to be revised
 460
 461 void AliESDv0::GetPosCov(Double_t cov[6]) const {
 462
 463   for (Int_t i=0; i<6; ++i) cov[i] = fPosCov[i];
 464
 465 }
 466
 467 Double_t AliESDv0::GetSigmaY(){
 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 sigmay = cp[0]+cm[0]+ cp[5]*(fParamP.GetX()-fRr)*(fParamP.GetX()-fRr)+ cm[5]*(fParamN.GetX()-fRr)*(fParamN.GetX()-fRr);
 474   return (sigmay>0) ? TMath::Sqrt(sigmay):100;
 475 }
 476
 477 Double_t AliESDv0::GetSigmaZ(){
 478   //
 479   // return sigmay in y  at vertex position  using covariance matrix
 480   //
 481   const Double_t * cp  = fParamP.GetCovariance();
 482   const Double_t * cm  = fParamN.GetCovariance();
 483   Double_t sigmaz = cp[2]+cm[2]+ cp[9]*(fParamP.GetX()-fRr)*(fParamP.GetX()-fRr)+ cm[9]*(fParamN.GetX()-fRr)*(fParamN.GetX()-fRr);
 484   return (sigmaz>0) ? TMath::Sqrt(sigmaz):100;
 485 }
 486
 487 Double_t AliESDv0::GetSigmaD0(){
 488   //
 489   // Sigma parameterization using covariance matrix
 490   //
 491   // sigma of distance between two tracks in vertex position
 492   // sigma of DCA is proportianal to sigmaD0
 493   // factor 2 difference is explained by the fact that the DCA is calculated at the position
 494   // where the tracks as closest together ( not exact position of the vertex)
 495   //
 496   const Double_t * cp      = fParamP.GetCovariance();
 497   const Double_t * cm      = fParamN.GetCovariance();
 498   Double_t sigmaD0   = cp[0]+cm[0]+cp[2]+cm[2]+fgkParams.fPSigmaOffsetD0*fgkParams.fPSigmaOffsetD0;
 499   sigmaD0           += ((fParamP.GetX()-fRr)*(fParamP.GetX()-fRr))*(cp[5]+cp[9]);
 500   sigmaD0           += ((fParamN.GetX()-fRr)*(fParamN.GetX()-fRr))*(cm[5]+cm[9]);
 501   return (sigmaD0>0)? TMath::Sqrt(sigmaD0):100;
 502 }
 503
 504
 505 Double_t AliESDv0::GetSigmaAP0(){
 506   //
 507   //Sigma parameterization using covariance matrices
 508   //
 509   Double_t prec  = TMath::Sqrt((fNmom[0]+fPmom[0])*(fNmom[0]+fPmom[0])
 510                               +(fNmom[1]+fPmom[1])*(fNmom[1]+fPmom[1])
 511                               +(fNmom[2]+fPmom[2])*(fNmom[2]+fPmom[2]));
 512   Double_t normp = TMath::Sqrt(fPmom[0]*fPmom[0]+fPmom[1]*fPmom[1]+fPmom[2]*fPmom[2])/prec;  // fraction of the momenta
 513   Double_t normm = TMath::Sqrt(fNmom[0]*fNmom[0]+fNmom[1]*fNmom[1]+fNmom[2]*fNmom[2])/prec;
 514   const Double_t * cp      = fParamP.GetCovariance();
 515   const Double_t * cm      = fParamN.GetCovariance();
 516   Double_t sigmaAP0 = fgkParams.fPSigmaOffsetAP0*fgkParams.fPSigmaOffsetAP0;                           // minimal part
 517   sigmaAP0 +=  (cp[5]+cp[9])*(normp*normp)+(cm[5]+cm[9])*(normm*normm);          // angular resolution part
 518   Double_t sigmaAP1 = GetSigmaD0()/(TMath::Abs(fRr)+0.01);                       // vertex position part
 519   sigmaAP0 +=  0.5*sigmaAP1*sigmaAP1;
 520   return (sigmaAP0>0)? TMath::Sqrt(sigmaAP0):100;
 521 }
 522
 523 Double_t AliESDv0::GetEffectiveSigmaD0(){
 524   //
 525   // minimax - effective Sigma parameterization
 526   // p12 effective curvature and v0 radius postion used as parameters
 527   //
 528   Double_t p12 = TMath::Sqrt(fParamP.GetParameter()[4]*fParamP.GetParameter()[4]+
 529                              fParamN.GetParameter()[4]*fParamN.GetParameter()[4]);
 530   Double_t sigmaED0= TMath::Max(TMath::Sqrt(fRr)-fgkParams.fPSigmaRminDE,0.0)*fgkParams.fPSigmaCoefDE*p12*p12;
 531   sigmaED0*= sigmaED0;
 532   sigmaED0*= sigmaED0;
 533   sigmaED0 = TMath::Sqrt(sigmaED0+fgkParams.fPSigmaOffsetDE*fgkParams.fPSigmaOffsetDE);
 534   return (sigmaED0<fgkParams.fPSigmaMaxDE) ? sigmaED0: fgkParams.fPSigmaMaxDE;
 535 }
 536
 537
 538 Double_t AliESDv0::GetEffectiveSigmaAP0(){
 539   //
 540   // effective Sigma parameterization of point angle resolution
 541   //
 542   Double_t p12 = TMath::Sqrt(fParamP.GetParameter()[4]*fParamP.GetParameter()[4]+
 543                              fParamN.GetParameter()[4]*fParamN.GetParameter()[4]);
 544   Double_t sigmaAPE= fgkParams.fPSigmaBase0APE;
 545   sigmaAPE+= fgkParams.fPSigmaR0APE/(fgkParams.fPSigmaR1APE+fRr);
 546   sigmaAPE*= (fgkParams.fPSigmaP0APE+fgkParams.fPSigmaP1APE*p12);
 547   sigmaAPE = TMath::Min(sigmaAPE,fgkParams.fPSigmaMaxAPE);
 548   return sigmaAPE;
 549 }
 550
 551
 552 Double_t  AliESDv0::GetMinimaxSigmaAP0(){
 553   //
 554   // calculate mini-max effective sigma of point angle resolution
 555   //
 556   //compv0->fTree->SetAlias("SigmaAP2","max(min((SigmaAP0+SigmaAPE0)*0.5,1.5*SigmaAPE0),0.5*SigmaAPE0+0.003)");
 557   Double_t    effectiveSigma = GetEffectiveSigmaAP0();
 558   Double_t    sigmaMMAP = 0.5*(GetSigmaAP0()+effectiveSigma);
 559   sigmaMMAP  = TMath::Min(sigmaMMAP, fgkParams.fPMaxFractionAP0*effectiveSigma);
 560   sigmaMMAP  = TMath::Max(sigmaMMAP, fgkParams.fPMinFractionAP0*effectiveSigma+fgkParams.fPMinAP0);
 561   return sigmaMMAP;
 562 }
 563 Double_t  AliESDv0::GetMinimaxSigmaD0(){
 564   //
 565   // calculate mini-max sigma of dca resolution
 566   //
 567   //compv0->fTree->SetAlias("SigmaD2","max(min((SigmaD0+SigmaDE0)*0.5,1.5*SigmaDE0),0.5*SigmaDE0)");
 568   Double_t    effectiveSigma = GetEffectiveSigmaD0();
 569   Double_t    sigmaMMD0 = 0.5*(GetSigmaD0()+effectiveSigma);
 570   sigmaMMD0  = TMath::Min(sigmaMMD0, fgkParams.fPMaxFractionD0*effectiveSigma);
 571   sigmaMMD0  = TMath::Max(sigmaMMD0, fgkParams.fPMinFractionD0*effectiveSigma+fgkParams.fPMinD0);
 572   return sigmaMMD0;
 573 }
 574
 575
 576 Double_t AliESDv0::GetLikelihoodAP(Int_t mode0, Int_t mode1){
 577   //
 578   // get likelihood for point angle
 579   //
 580   Double_t sigmaAP = 0.007;            //default sigma
 581   switch (mode0){
 582   case 0:
 583     sigmaAP = GetSigmaAP0();           // mode 0  - covariance matrix estimates used
 584     break;
 585   case 1:
 586     sigmaAP = GetEffectiveSigmaAP0();  // mode 1 - effective sigma used
 587     break;
 588   case 2:
 589     sigmaAP = GetMinimaxSigmaAP0();    // mode 2 - minimax sigma
 590     break;
 591   }
 592   Double_t apNorm = TMath::Min(TMath::ACos(fPointAngle)/sigmaAP,50.);
 593   //normalized point angle, restricted - because of overflow problems in Exp
 594   Double_t likelihood = 0;
 595   switch(mode1){
 596   case 0:
 597     likelihood = TMath::Exp(-0.5*apNorm*apNorm);
 598     // one component
 599     break;
 600   case 1:
 601     likelihood = (TMath::Exp(-0.5*apNorm*apNorm)+0.5* TMath::Exp(-0.25*apNorm*apNorm))/1.5;
 602     // two components
 603     break;
 604   case 2:
 605     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;
 606     // three components
 607     break;
 608   }
 609   return likelihood;
 610 }
 611
 612 Double_t AliESDv0::GetLikelihoodD(Int_t mode0, Int_t mode1){
 613   //
 614   // get likelihood for DCA
 615   //
 616   Double_t sigmaD = 0.03;            //default sigma
 617   switch (mode0){
 618   case 0:
 619     sigmaD = GetSigmaD0();           // mode 0  - covariance matrix estimates used
 620     break;
 621   case 1:
 622     sigmaD = GetEffectiveSigmaD0();  // mode 1 - effective sigma used
 623     break;
 624   case 2:
 625     sigmaD = GetMinimaxSigmaD0();    // mode 2 - minimax sigma
 626     break;
 627   }
 628
 629   //Bo:  Double_t dNorm = TMath::Min(fDist2/sigmaD,50.);
 630   Double_t dNorm = TMath::Min(fDcaV0Daughters/sigmaD,50.);//Bo:
 631   //normalized point angle, restricted - because of overflow problems in Exp
 632   Double_t likelihood = 0;
 633   switch(mode1){
 634   case 0:
 635     likelihood = TMath::Exp(-2.*dNorm);
 636     // one component
 637     break;
 638   case 1:
 639     likelihood = (TMath::Exp(-2.*dNorm)+0.5* TMath::Exp(-dNorm))/1.5;
 640     // two components
 641     break;
 642   case 2:
 643     likelihood = (TMath::Exp(-2.*dNorm)+0.5* TMath::Exp(-dNorm)+0.25*TMath::Exp(-0.5*dNorm))/1.75;
 644     // three components
 645     break;
 646   }
 647   return likelihood;
 648
 649 }
 650
 651 Double_t AliESDv0::GetLikelihoodC(Int_t mode0, Int_t /*mode1*/) const {
 652   //
 653   // get likelihood for Causality
 654   // !!!  Causality variables defined in AliITStrackerMI !!!
 655   //      when more information was available
 656   //
 657   Double_t likelihood = 0.5;
 658   Double_t minCausal  = TMath::Min(fCausality[0],fCausality[1]);
 659   Double_t maxCausal  = TMath::Max(fCausality[0],fCausality[1]);
 660   //  minCausal           = TMath::Max(minCausal,0.5*maxCausal);
 661   //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");
 662
 663   switch(mode0){
 664   case 0:
 665     //normalization
 666     likelihood = TMath::Power((1.05-2*(0.8-TMath::Exp(-maxCausal))),4.);
 667     break;
 668   case 1:
 669     likelihood = TMath::Power(1.05-(2*(0.8-TMath::Exp(-maxCausal))+(2*(0.8-TMath::Exp(-minCausal))))*0.5,4.);
 670     break;
 671   }
 672   return likelihood;
 673
 674 }
 675
 676 void AliESDv0::SetCausality(Float_t pb0, Float_t pb1, Float_t pa0, Float_t pa1)
 677 {
 678   //
 679   // set probabilities
 680   //
 681   fCausality[0] = pb0;     // probability - track 0 exist before vertex
 682   fCausality[1] = pb1;     // probability - track 1 exist before vertex
 683   fCausality[2] = pa0;     // probability - track 0 exist close after vertex
 684   fCausality[3] = pa1;     // probability - track 1 exist close after vertex
 685 }
 686 void  AliESDv0::SetClusters(const Int_t *clp, const Int_t *clm)
 687 {
 688   //
 689   // Set its clusters indexes
 690   //
 691   for (Int_t i=0;i<6;i++) fClusters[0][i] = clp[i];
 692   for (Int_t i=0;i<6;i++) fClusters[1][i] = clm[i];
 693 }
 694
 695 Double_t AliESDv0::GetEffMass(UInt_t p1, UInt_t p2) const{
 696   //
 697   // calculate effective mass
 698   //
 699   const Double_t kpmass[5] = {TDatabasePDG::Instance()->GetParticle(kElectron)->Mass(),
 700                              TDatabasePDG::Instance()->GetParticle(kMuonMinus)->Mass(),
 701                              TDatabasePDG::Instance()->GetParticle(kPiPlus)->Mass(),
 702                              TDatabasePDG::Instance()->GetParticle(kKPlus)->Mass(),
 703                              TDatabasePDG::Instance()->GetParticle(kProton)->Mass()};
 704   /*
 705   if (p1>4) return -1;
 706   if (p2>4) return -1;
 707   Double_t mass1 = kpmass[p1];
 708   Double_t mass2 = kpmass[p2];
 709   const Double_t *m1 = fPmom;
 710   const Double_t *m2 = fNmom;
 711   //
 712   //if (fRP[p1]+fRM[p2]<fRP[p2]+fRM[p1]){
 713   //  m1 = fPM;
 714   //  m2 = fPP;
 715   //}
 716   //
 717   Double_t e1    = TMath::Sqrt(mass1*mass1+
 718                               m1[0]*m1[0]+
 719                               m1[1]*m1[1]+
 720                               m1[2]*m1[2]);
 721   Double_t e2    = TMath::Sqrt(mass2*mass2+
 722                               m2[0]*m2[0]+
 723                               m2[1]*m2[1]+
 724                               m2[2]*m2[2]);
 725   Double_t mass =
 726     (m2[0]+m1[0])*(m2[0]+m1[0])+
 727     (m2[1]+m1[1])*(m2[1]+m1[1])+
 728     (m2[2]+m1[2])*(m2[2]+m1[2]);
 729
 730   mass = (e1+e2)*(e1+e2)-mass;
 731   if (mass < 0.) mass = 0.;
 732   return (TMath::Sqrt(mass));
 733   */
 734   if(p1>4 || p2>4) return -1;
 735   Double_t e12   = kpmass[p1]*kpmass[p1]+fPmom[0]*fPmom[0]+fPmom[1]*fPmom[1]+fPmom[2]*fPmom[2];
 736   Double_t e22   = kpmass[p2]*kpmass[p2]+fNmom[0]*fNmom[0]+fNmom[1]*fNmom[1]+fNmom[2]*fNmom[2];
 737   Double_t cmass = TMath::Sqrt(TMath::Max(kpmass[p1]*kpmass[p1]+kpmass[p2]*kpmass[p2]
 738                                           +2.*(TMath::Sqrt(e12*e22)-fPmom[0]*fNmom[0]-fPmom[1]*fNmom[1]-fPmom[2]*fNmom[2]),0.));
 739   return cmass;
 740
 741 }