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