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   return (lQlPos - lQlNeg)/(lQlPos + lQlNeg);
 325 }
 326
 327 Double_t AliESDv0::PtArmV0() const {
 328   //--------------------------------------------------------------------
 329   // This gives the Armenteros-Podolanski ptarm
 330   //--------------------------------------------------------------------
 331   TVector3 momNeg(fNmom[0],fNmom[1],fNmom[2]);
 332   TVector3 momTot(Px(),Py(),Pz());
 333
 334   return momNeg.Perp(momTot);
 335 }
 336
 337 // Eventually the older functions
 338 Double_t AliESDv0::ChangeMassHypothesis(Int_t code) {
 339   //--------------------------------------------------------------------
 340   // This function changes the mass hypothesis for this V0
 341   // and returns the "kinematical quality" of this hypothesis
 342   //--------------------------------------------------------------------
 343   static
 344   Double_t piMass=TDatabasePDG::Instance()->GetParticle(kPiPlus)->Mass();
 345   static
 346   Double_t prMass=TDatabasePDG::Instance()->GetParticle(kProton)->Mass();
 347   static
 348   Double_t k0Mass=TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
 349   static
 350   Double_t l0Mass=TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
 351
 352   Double_t nmass=piMass, pmass=piMass, mass=k0Mass, ps=0.206;
 353
 354   fPdgCode=code;
 355
 356   switch (code) {
 357   case kLambda0:
 358     nmass=piMass; pmass=prMass; mass=l0Mass; ps=0.101; break;
 359   case kLambda0Bar:
 360     pmass=piMass; nmass=prMass; mass=l0Mass; ps=0.101; break;
 361   case kK0Short:
 362     break;
 363   default:
 364     AliError("invalide PDG code ! Assuming K0s...");
 365     fPdgCode=kK0Short;
 366     break;
 367   }
 368
 369   Double_t pxn=fNmom[0], pyn=fNmom[1], pzn=fNmom[2];
 370   Double_t pxp=fPmom[0], pyp=fPmom[1], pzp=fPmom[2];
 371
 372   Double_t en=TMath::Sqrt(nmass*nmass + pxn*pxn + pyn*pyn + pzn*pzn);
 373   Double_t ep=TMath::Sqrt(pmass*pmass + pxp*pxp + pyp*pyp + pzp*pzp);
 374   Double_t pxl=pxn+pxp, pyl=pyn+pyp, pzl=pzn+pzp;
 375   Double_t pl=TMath::Sqrt(pxl*pxl + pyl*pyl + pzl*pzl);
 376
 377   fEffMass=TMath::Sqrt((en+ep)*(en+ep)-pl*pl);
 378
 379   Double_t beta=pl/(en+ep);
 380   Double_t pln=(pxn*pxl + pyn*pyl + pzn*pzl)/pl;
 381   Double_t plp=(pxp*pxl + pyp*pyl + pzp*pzl)/pl;
 382
 383   Double_t pt2=pxp*pxp + pyp*pyp + pzp*pzp - plp*plp;
 384
 385   Double_t a=(plp-pln)/(plp+pln);
 386   a -= (pmass*pmass-nmass*nmass)/(mass*mass);
 387   a = 0.25*beta*beta*mass*mass*a*a + pt2;
 388
 389   return (a - ps*ps);
 390
 391 }
 392
 393 void AliESDv0::GetPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const {
 394   //--------------------------------------------------------------------
 395   // This function returns V0's momentum (global)
 396   //--------------------------------------------------------------------
 397   px=fNmom[0]+fPmom[0];
 398   py=fNmom[1]+fPmom[1];
 399   pz=fNmom[2]+fPmom[2];
 400 }
 401
 402 void AliESDv0::GetXYZ(Double_t &x, Double_t &y, Double_t &z) const {
 403   //--------------------------------------------------------------------
 404   // This function returns V0's position (global)
 405   //--------------------------------------------------------------------
 406   x=fPos[0];
 407   y=fPos[1];
 408   z=fPos[2];
 409 }
 410
 411 Float_t AliESDv0::GetD(Double_t x0, Double_t y0, Double_t z0) const {
 412   //--------------------------------------------------------------------
 413   // This function returns V0's impact parameter
 414   //--------------------------------------------------------------------
 415   Double_t x=fPos[0],y=fPos[1],z=fPos[2];
 416   Double_t px=fNmom[0]+fPmom[0];
 417   Double_t py=fNmom[1]+fPmom[1];
 418   Double_t pz=fNmom[2]+fPmom[2];
 419
 420   Double_t dx=(y0-y)*pz - (z0-z)*py;
 421   Double_t dy=(x0-x)*pz - (z0-z)*px;
 422   Double_t dz=(x0-x)*py - (y0-y)*px;
 423   Double_t d=TMath::Sqrt((dx*dx+dy*dy+dz*dz)/(px*px+py*py+pz*pz));
 424   return d;
 425 }
 426
 427 Float_t AliESDv0::GetV0CosineOfPointingAngle(Double_t refPointX, Double_t refPointY, Double_t refPointZ) const {
 428   // calculates the pointing angle of the V0 wrt a reference point
 429
 430   Double_t momV0[3]; //momentum of the V0
 431   GetPxPyPz(momV0[0],momV0[1],momV0[2]);
 432
 433   Double_t deltaPos[3]; //vector between the reference point and the V0 vertex
 434   deltaPos[0] = fPos[0] - refPointX;
 435   deltaPos[1] = fPos[1] - refPointY;
 436   deltaPos[2] = fPos[2] - refPointZ;
 437
 438   Double_t momV02    = momV0[0]*momV0[0] + momV0[1]*momV0[1] + momV0[2]*momV0[2];
 439   Double_t deltaPos2 = deltaPos[0]*deltaPos[0] + deltaPos[1]*deltaPos[1] + deltaPos[2]*deltaPos[2];
 440
 441   Double_t cosinePointingAngle = (deltaPos[0]*momV0[0] +
 442                                   deltaPos[1]*momV0[1] +
 443                                   deltaPos[2]*momV0[2] ) /
 444     TMath::Sqrt(momV02 * deltaPos2);
 445
 446   return cosinePointingAngle;
 447 }
 448
 449
 450 // **** The following functions need to be revised
 451
 452 void AliESDv0::GetPosCov(Double_t cov[6]) const {
 453
 454   for (Int_t i=0; i<6; ++i) cov[i] = fPosCov[i];
 455
 456 }
 457
 458 Double_t AliESDv0::GetSigmaY(){
 459   //
 460   // return sigmay in y  at vertex position  using covariance matrix
 461   //
 462   const Double_t * cp  = fParamP.GetCovariance();
 463   const Double_t * cm  = fParamN.GetCovariance();
 464   Double_t sigmay = cp[0]+cm[0]+ cp[5]*(fParamP.GetX()-fRr)*(fParamP.GetX()-fRr)+ cm[5]*(fParamN.GetX()-fRr)*(fParamN.GetX()-fRr);
 465   return (sigmay>0) ? TMath::Sqrt(sigmay):100;
 466 }
 467
 468 Double_t AliESDv0::GetSigmaZ(){
 469   //
 470   // return sigmay in y  at vertex position  using covariance matrix
 471   //
 472   const Double_t * cp  = fParamP.GetCovariance();
 473   const Double_t * cm  = fParamN.GetCovariance();
 474   Double_t sigmaz = cp[2]+cm[2]+ cp[9]*(fParamP.GetX()-fRr)*(fParamP.GetX()-fRr)+ cm[9]*(fParamN.GetX()-fRr)*(fParamN.GetX()-fRr);
 475   return (sigmaz>0) ? TMath::Sqrt(sigmaz):100;
 476 }
 477
 478 Double_t AliESDv0::GetSigmaD0(){
 479   //
 480   // Sigma parameterization using covariance matrix
 481   //
 482   // sigma of distance between two tracks in vertex position
 483   // sigma of DCA is proportianal to sigmaD0
 484   // factor 2 difference is explained by the fact that the DCA is calculated at the position
 485   // where the tracks as closest together ( not exact position of the vertex)
 486   //
 487   const Double_t * cp      = fParamP.GetCovariance();
 488   const Double_t * cm      = fParamN.GetCovariance();
 489   Double_t sigmaD0   = cp[0]+cm[0]+cp[2]+cm[2]+fgkParams.fPSigmaOffsetD0*fgkParams.fPSigmaOffsetD0;
 490   sigmaD0           += ((fParamP.GetX()-fRr)*(fParamP.GetX()-fRr))*(cp[5]+cp[9]);
 491   sigmaD0           += ((fParamN.GetX()-fRr)*(fParamN.GetX()-fRr))*(cm[5]+cm[9]);
 492   return (sigmaD0>0)? TMath::Sqrt(sigmaD0):100;
 493 }
 494
 495
 496 Double_t AliESDv0::GetSigmaAP0(){
 497   //
 498   //Sigma parameterization using covariance matrices
 499   //
 500   Double_t prec  = TMath::Sqrt((fNmom[0]+fPmom[0])*(fNmom[0]+fPmom[0])
 501                               +(fNmom[1]+fPmom[1])*(fNmom[1]+fPmom[1])
 502                               +(fNmom[2]+fPmom[2])*(fNmom[2]+fPmom[2]));
 503   Double_t normp = TMath::Sqrt(fPmom[0]*fPmom[0]+fPmom[1]*fPmom[1]+fPmom[2]*fPmom[2])/prec;  // fraction of the momenta
 504   Double_t normm = TMath::Sqrt(fNmom[0]*fNmom[0]+fNmom[1]*fNmom[1]+fNmom[2]*fNmom[2])/prec;
 505   const Double_t * cp      = fParamP.GetCovariance();
 506   const Double_t * cm      = fParamN.GetCovariance();
 507   Double_t sigmaAP0 = fgkParams.fPSigmaOffsetAP0*fgkParams.fPSigmaOffsetAP0;                           // minimal part
 508   sigmaAP0 +=  (cp[5]+cp[9])*(normp*normp)+(cm[5]+cm[9])*(normm*normm);          // angular resolution part
 509   Double_t sigmaAP1 = GetSigmaD0()/(TMath::Abs(fRr)+0.01);                       // vertex position part
 510   sigmaAP0 +=  0.5*sigmaAP1*sigmaAP1;
 511   return (sigmaAP0>0)? TMath::Sqrt(sigmaAP0):100;
 512 }
 513
 514 Double_t AliESDv0::GetEffectiveSigmaD0(){
 515   //
 516   // minimax - effective Sigma parameterization
 517   // p12 effective curvature and v0 radius postion used as parameters
 518   //
 519   Double_t p12 = TMath::Sqrt(fParamP.GetParameter()[4]*fParamP.GetParameter()[4]+
 520                              fParamN.GetParameter()[4]*fParamN.GetParameter()[4]);
 521   Double_t sigmaED0= TMath::Max(TMath::Sqrt(fRr)-fgkParams.fPSigmaRminDE,0.0)*fgkParams.fPSigmaCoefDE*p12*p12;
 522   sigmaED0*= sigmaED0;
 523   sigmaED0*= sigmaED0;
 524   sigmaED0 = TMath::Sqrt(sigmaED0+fgkParams.fPSigmaOffsetDE*fgkParams.fPSigmaOffsetDE);
 525   return (sigmaED0<fgkParams.fPSigmaMaxDE) ? sigmaED0: fgkParams.fPSigmaMaxDE;
 526 }
 527
 528
 529 Double_t AliESDv0::GetEffectiveSigmaAP0(){
 530   //
 531   // effective Sigma parameterization of point angle resolution
 532   //
 533   Double_t p12 = TMath::Sqrt(fParamP.GetParameter()[4]*fParamP.GetParameter()[4]+
 534                              fParamN.GetParameter()[4]*fParamN.GetParameter()[4]);
 535   Double_t sigmaAPE= fgkParams.fPSigmaBase0APE;
 536   sigmaAPE+= fgkParams.fPSigmaR0APE/(fgkParams.fPSigmaR1APE+fRr);
 537   sigmaAPE*= (fgkParams.fPSigmaP0APE+fgkParams.fPSigmaP1APE*p12);
 538   sigmaAPE = TMath::Min(sigmaAPE,fgkParams.fPSigmaMaxAPE);
 539   return sigmaAPE;
 540 }
 541
 542
 543 Double_t  AliESDv0::GetMinimaxSigmaAP0(){
 544   //
 545   // calculate mini-max effective sigma of point angle resolution
 546   //
 547   //compv0->fTree->SetAlias("SigmaAP2","max(min((SigmaAP0+SigmaAPE0)*0.5,1.5*SigmaAPE0),0.5*SigmaAPE0+0.003)");
 548   Double_t    effectiveSigma = GetEffectiveSigmaAP0();
 549   Double_t    sigmaMMAP = 0.5*(GetSigmaAP0()+effectiveSigma);
 550   sigmaMMAP  = TMath::Min(sigmaMMAP, fgkParams.fPMaxFractionAP0*effectiveSigma);
 551   sigmaMMAP  = TMath::Max(sigmaMMAP, fgkParams.fPMinFractionAP0*effectiveSigma+fgkParams.fPMinAP0);
 552   return sigmaMMAP;
 553 }
 554 Double_t  AliESDv0::GetMinimaxSigmaD0(){
 555   //
 556   // calculate mini-max sigma of dca resolution
 557   //
 558   //compv0->fTree->SetAlias("SigmaD2","max(min((SigmaD0+SigmaDE0)*0.5,1.5*SigmaDE0),0.5*SigmaDE0)");
 559   Double_t    effectiveSigma = GetEffectiveSigmaD0();
 560   Double_t    sigmaMMD0 = 0.5*(GetSigmaD0()+effectiveSigma);
 561   sigmaMMD0  = TMath::Min(sigmaMMD0, fgkParams.fPMaxFractionD0*effectiveSigma);
 562   sigmaMMD0  = TMath::Max(sigmaMMD0, fgkParams.fPMinFractionD0*effectiveSigma+fgkParams.fPMinD0);
 563   return sigmaMMD0;
 564 }
 565
 566
 567 Double_t AliESDv0::GetLikelihoodAP(Int_t mode0, Int_t mode1){
 568   //
 569   // get likelihood for point angle
 570   //
 571   Double_t sigmaAP = 0.007;            //default sigma
 572   switch (mode0){
 573   case 0:
 574     sigmaAP = GetSigmaAP0();           // mode 0  - covariance matrix estimates used
 575     break;
 576   case 1:
 577     sigmaAP = GetEffectiveSigmaAP0();  // mode 1 - effective sigma used
 578     break;
 579   case 2:
 580     sigmaAP = GetMinimaxSigmaAP0();    // mode 2 - minimax sigma
 581     break;
 582   }
 583   Double_t apNorm = TMath::Min(TMath::ACos(fPointAngle)/sigmaAP,50.);
 584   //normalized point angle, restricted - because of overflow problems in Exp
 585   Double_t likelihood = 0;
 586   switch(mode1){
 587   case 0:
 588     likelihood = TMath::Exp(-0.5*apNorm*apNorm);
 589     // one component
 590     break;
 591   case 1:
 592     likelihood = (TMath::Exp(-0.5*apNorm*apNorm)+0.5* TMath::Exp(-0.25*apNorm*apNorm))/1.5;
 593     // two components
 594     break;
 595   case 2:
 596     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;
 597     // three components
 598     break;
 599   }
 600   return likelihood;
 601 }
 602
 603 Double_t AliESDv0::GetLikelihoodD(Int_t mode0, Int_t mode1){
 604   //
 605   // get likelihood for DCA
 606   //
 607   Double_t sigmaD = 0.03;            //default sigma
 608   switch (mode0){
 609   case 0:
 610     sigmaD = GetSigmaD0();           // mode 0  - covariance matrix estimates used
 611     break;
 612   case 1:
 613     sigmaD = GetEffectiveSigmaD0();  // mode 1 - effective sigma used
 614     break;
 615   case 2:
 616     sigmaD = GetMinimaxSigmaD0();    // mode 2 - minimax sigma
 617     break;
 618   }
 619
 620   //Bo:  Double_t dNorm = TMath::Min(fDist2/sigmaD,50.);
 621   Double_t dNorm = TMath::Min(fDcaV0Daughters/sigmaD,50.);//Bo:
 622   //normalized point angle, restricted - because of overflow problems in Exp
 623   Double_t likelihood = 0;
 624   switch(mode1){
 625   case 0:
 626     likelihood = TMath::Exp(-2.*dNorm);
 627     // one component
 628     break;
 629   case 1:
 630     likelihood = (TMath::Exp(-2.*dNorm)+0.5* TMath::Exp(-dNorm))/1.5;
 631     // two components
 632     break;
 633   case 2:
 634     likelihood = (TMath::Exp(-2.*dNorm)+0.5* TMath::Exp(-dNorm)+0.25*TMath::Exp(-0.5*dNorm))/1.75;
 635     // three components
 636     break;
 637   }
 638   return likelihood;
 639
 640 }
 641
 642 Double_t AliESDv0::GetLikelihoodC(Int_t mode0, Int_t /*mode1*/) const {
 643   //
 644   // get likelihood for Causality
 645   // !!!  Causality variables defined in AliITStrackerMI !!!
 646   //      when more information was available
 647   //
 648   Double_t likelihood = 0.5;
 649   Double_t minCausal  = TMath::Min(fCausality[0],fCausality[1]);
 650   Double_t maxCausal  = TMath::Max(fCausality[0],fCausality[1]);
 651   //  minCausal           = TMath::Max(minCausal,0.5*maxCausal);
 652   //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");
 653
 654   switch(mode0){
 655   case 0:
 656     //normalization
 657     likelihood = TMath::Power((1.05-2*(0.8-TMath::Exp(-maxCausal))),4.);
 658     break;
 659   case 1:
 660     likelihood = TMath::Power(1.05-(2*(0.8-TMath::Exp(-maxCausal))+(2*(0.8-TMath::Exp(-minCausal))))*0.5,4.);
 661     break;
 662   }
 663   return likelihood;
 664
 665 }
 666
 667 void AliESDv0::SetCausality(Float_t pb0, Float_t pb1, Float_t pa0, Float_t pa1)
 668 {
 669   //
 670   // set probabilities
 671   //
 672   fCausality[0] = pb0;     // probability - track 0 exist before vertex
 673   fCausality[1] = pb1;     // probability - track 1 exist before vertex
 674   fCausality[2] = pa0;     // probability - track 0 exist close after vertex
 675   fCausality[3] = pa1;     // probability - track 1 exist close after vertex
 676 }
 677 void  AliESDv0::SetClusters(const Int_t *clp, const Int_t *clm)
 678 {
 679   //
 680   // Set its clusters indexes
 681   //
 682   for (Int_t i=0;i<6;i++) fClusters[0][i] = clp[i];
 683   for (Int_t i=0;i<6;i++) fClusters[1][i] = clm[i];
 684 }
 685
 686 Double_t AliESDv0::GetEffMass(UInt_t p1, UInt_t p2) const{
 687   //
 688   // calculate effective mass
 689   //
 690   const Double_t kpmass[5] = {TDatabasePDG::Instance()->GetParticle(kElectron)->Mass(),
 691                              TDatabasePDG::Instance()->GetParticle(kMuonMinus)->Mass(),
 692                              TDatabasePDG::Instance()->GetParticle(kPiPlus)->Mass(),
 693                              TDatabasePDG::Instance()->GetParticle(kKPlus)->Mass(),
 694                              TDatabasePDG::Instance()->GetParticle(kProton)->Mass()};
 695   /*
 696   if (p1>4) return -1;
 697   if (p2>4) return -1;
 698   Double_t mass1 = kpmass[p1];
 699   Double_t mass2 = kpmass[p2];
 700   const Double_t *m1 = fPmom;
 701   const Double_t *m2 = fNmom;
 702   //
 703   //if (fRP[p1]+fRM[p2]<fRP[p2]+fRM[p1]){
 704   //  m1 = fPM;
 705   //  m2 = fPP;
 706   //}
 707   //
 708   Double_t e1    = TMath::Sqrt(mass1*mass1+
 709                               m1[0]*m1[0]+
 710                               m1[1]*m1[1]+
 711                               m1[2]*m1[2]);
 712   Double_t e2    = TMath::Sqrt(mass2*mass2+
 713                               m2[0]*m2[0]+
 714                               m2[1]*m2[1]+
 715                               m2[2]*m2[2]);
 716   Double_t mass =
 717     (m2[0]+m1[0])*(m2[0]+m1[0])+
 718     (m2[1]+m1[1])*(m2[1]+m1[1])+
 719     (m2[2]+m1[2])*(m2[2]+m1[2]);
 720
 721   mass = (e1+e2)*(e1+e2)-mass;
 722   if (mass < 0.) mass = 0.;
 723   return (TMath::Sqrt(mass));
 724   */
 725   if(p1>4 || p2>4) return -1;
 726   Double_t e12   = kpmass[p1]*kpmass[p1]+fPmom[0]*fPmom[0]+fPmom[1]*fPmom[1]+fPmom[2]*fPmom[2];
 727   Double_t e22   = kpmass[p2]*kpmass[p2]+fNmom[0]*fNmom[0]+fNmom[1]*fNmom[1]+fNmom[2]*fNmom[2];
 728   Double_t cmass = TMath::Sqrt(TMath::Max(kpmass[p1]*kpmass[p1]+kpmass[p2]*kpmass[p2]
 729                                           +2.*(TMath::Sqrt(e12*e22)-fPmom[0]*fNmom[0]-fPmom[1]*fNmom[1]-fPmom[2]*fNmom[2]),0.));
 730   return cmass;
 731
 732 }