STEER/AliESDV0MI.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 //
  20 //    Implementation of the ESD V0MI vertex class
  21 //            This class is part of the Event Data Summary
  22 //            set of classes and contains information about
  23 //            V0 kind vertexes generated by a neutral particle
  24 //    Numerical part - AliHelix functionality used
  25 //
  26 //    Likelihoods for Point angle, DCA and Causality defined => can be used as cut parameters
  27 //    HIGHLY recomended
  28 //
  29 //    Quality information can be used as additional cut variables
  30 //
  31 //    Origin: Marian Ivanov marian.ivanov@cern.ch
  32 //-------------------------------------------------------------------------
  33
  34 #include <Riostream.h>
  35 #include <TMath.h>
  36
  37 #include "AliESDV0MI.h"
  38 #include "AliHelix.h"
  39
  40
  41 ClassImp(AliESDV0MI)
  42
  43 AliESDV0MIParams  AliESDV0MI::fgkParams;
  44
  45
  46 AliESDV0MI::AliESDV0MI() :
  47   AliESDv0(),
  48   fParamP(),
  49   fParamM(),
  50   fID(0),
  51   fDist1(-1),
  52   fDist2(-1),
  53   fRr(-1),
  54   fStatus(0),
  55   fRow0(-1),
  56   fDistNorm(0),
  57   fDistSigma(0),
  58   fChi2Before(0),
  59   fNBefore(0),
  60   fChi2After(0),
  61   fNAfter(0),
  62   fPointAngleFi(0),
  63   fPointAngleTh(0),
  64   fPointAngle(0)
  65 {
  66   //
  67   //Dafault constructor
  68   //
  69   for (Int_t i=0;i<5;i++){
  70     fRP[i]=fRM[i]=0;
  71   }
  72   fLab[0]=fLab[1]=-1;
  73   fIndex[0]=fIndex[1]=-1;
  74   for (Int_t i=0;i<6;i++){fClusters[0][i]=0; fClusters[1][i]=0;}
  75   fNormDCAPrim[0]=fNormDCAPrim[1]=0;
  76   for (Int_t i=0;i<3;i++){fPP[i]=fPM[i]=fXr[i]=fAngle[i]=0;}
  77   for (Int_t i=0;i<3;i++){fOrder[i]=0;}
  78   for (Int_t i=0;i<4;i++){fCausality[i]=0;}
  79 }
  80
  81 Double_t AliESDV0MI::GetSigmaY(){
  82   //
  83   // return sigmay in y  at vertex position  using covariance matrix
  84   //
  85   const Double_t * cp  = fParamP.GetCovariance();
  86   const Double_t * cm  = fParamM.GetCovariance();
  87   Double_t sigmay = cp[0]+cm[0]+ cp[5]*(fParamP.X()-fRr)*(fParamP.X()-fRr)+ cm[5]*(fParamM.X()-fRr)*(fParamM.X()-fRr);
  88   return (sigmay>0) ? TMath::Sqrt(sigmay):100;
  89 }
  90
  91 Double_t AliESDV0MI::GetSigmaZ(){
  92   //
  93   // return sigmay in y  at vertex position  using covariance matrix
  94   //
  95   const Double_t * cp  = fParamP.GetCovariance();
  96   const Double_t * cm  = fParamM.GetCovariance();
  97   Double_t sigmaz = cp[2]+cm[2]+ cp[9]*(fParamP.X()-fRr)*(fParamP.X()-fRr)+ cm[9]*(fParamM.X()-fRr)*(fParamM.X()-fRr);
  98   return (sigmaz>0) ? TMath::Sqrt(sigmaz):100;
  99 }
 100
 101 Double_t AliESDV0MI::GetSigmaD0(){
 102   //
 103   // Sigma parameterization using covariance matrix
 104   //
 105   // sigma of distance between two tracks in vertex position
 106   // sigma of DCA is proportianal to sigmaD0
 107   // factor 2 difference is explained by the fact that the DCA is calculated at the position
 108   // where the tracks as closest together ( not exact position of the vertex)
 109   //
 110   const Double_t * cp      = fParamP.GetCovariance();
 111   const Double_t * cm      = fParamM.GetCovariance();
 112   Double_t sigmaD0   = cp[0]+cm[0]+cp[2]+cm[2]+fgkParams.fPSigmaOffsetD0*fgkParams.fPSigmaOffsetD0;
 113   sigmaD0           += ((fParamP.X()-fRr)*(fParamP.X()-fRr))*(cp[5]+cp[9]);
 114   sigmaD0           += ((fParamM.X()-fRr)*(fParamM.X()-fRr))*(cm[5]+cm[9]);
 115   return (sigmaD0>0)? TMath::Sqrt(sigmaD0):100;
 116 }
 117
 118
 119 Double_t AliESDV0MI::GetSigmaAP0(){
 120   //
 121   //Sigma parameterization using covariance matrices
 122   //
 123   Double_t prec  = TMath::Sqrt((fPM[0]+fPP[0])*(fPM[0]+fPP[0])
 124                                +(fPM[1]+fPP[1])*(fPM[1]+fPP[1])
 125                                +(fPM[2]+fPP[2])*(fPM[2]+fPP[2]));
 126   Double_t normp = TMath::Sqrt(fPP[0]*fPP[0]+fPP[1]*fPP[1]+fPP[2]*fPP[2])/prec;  // fraction of the momenta
 127   Double_t normm = TMath::Sqrt(fPM[0]*fPM[0]+fPM[1]*fPM[1]+fPM[2]*fPM[2])/prec;
 128   const Double_t * cp      = fParamP.GetCovariance();
 129   const Double_t * cm      = fParamM.GetCovariance();
 130   Double_t sigmaAP0 = fgkParams.fPSigmaOffsetAP0*fgkParams.fPSigmaOffsetAP0;                           // minimal part
 131   sigmaAP0 +=  (cp[5]+cp[9])*(normp*normp)+(cm[5]+cm[9])*(normm*normm);          // angular resolution part
 132   Double_t sigmaAP1 = GetSigmaD0()/(TMath::Abs(fRr)+0.01);                       // vertex position part
 133   sigmaAP0 +=  0.5*sigmaAP1*sigmaAP1;
 134   return (sigmaAP0>0)? TMath::Sqrt(sigmaAP0):100;
 135 }
 136
 137 Double_t AliESDV0MI::GetEffectiveSigmaD0(){
 138   //
 139   // minimax - effective Sigma parameterization
 140   // p12 effective curvature and v0 radius postion used as parameters
 141   //
 142   Double_t p12 = TMath::Sqrt(fParamP.GetParameter()[4]*fParamP.GetParameter()[4]+
 143                              fParamM.GetParameter()[4]*fParamM.GetParameter()[4]);
 144   Double_t sigmaED0= TMath::Max(TMath::Sqrt(fRr)-fgkParams.fPSigmaRminDE,0.0)*fgkParams.fPSigmaCoefDE*p12*p12;
 145   sigmaED0*= sigmaED0;
 146   sigmaED0*= sigmaED0;
 147   sigmaED0 = TMath::Sqrt(sigmaED0+fgkParams.fPSigmaOffsetDE*fgkParams.fPSigmaOffsetDE);
 148   return (sigmaED0<fgkParams.fPSigmaMaxDE) ? sigmaED0: fgkParams.fPSigmaMaxDE;
 149 }
 150
 151
 152 Double_t AliESDV0MI::GetEffectiveSigmaAP0(){
 153   //
 154   // effective Sigma parameterization of point angle resolution
 155   //
 156   Double_t p12 = TMath::Sqrt(fParamP.GetParameter()[4]*fParamP.GetParameter()[4]+
 157                              fParamM.GetParameter()[4]*fParamM.GetParameter()[4]);
 158   Double_t sigmaAPE= fgkParams.fPSigmaBase0APE;
 159   sigmaAPE+= fgkParams.fPSigmaR0APE/(fgkParams.fPSigmaR1APE+fRr);
 160   sigmaAPE*= (fgkParams.fPSigmaP0APE+fgkParams.fPSigmaP1APE*p12);
 161   sigmaAPE = TMath::Min(sigmaAPE,fgkParams.fPSigmaMaxAPE);
 162   return sigmaAPE;
 163 }
 164
 165
 166 Double_t  AliESDV0MI::GetMinimaxSigmaAP0(){
 167   //
 168   // calculate mini-max effective sigma of point angle resolution
 169   //
 170   //compv0->fTree->SetAlias("SigmaAP2","max(min((SigmaAP0+SigmaAPE0)*0.5,1.5*SigmaAPE0),0.5*SigmaAPE0+0.003)");
 171   Double_t    effectiveSigma = GetEffectiveSigmaAP0();
 172   Double_t    sigmaMMAP = 0.5*(GetSigmaAP0()+effectiveSigma);
 173   sigmaMMAP  = TMath::Min(sigmaMMAP, fgkParams.fPMaxFractionAP0*effectiveSigma);
 174   sigmaMMAP  = TMath::Max(sigmaMMAP, fgkParams.fPMinFractionAP0*effectiveSigma+fgkParams.fPMinAP0);
 175   return sigmaMMAP;
 176 }
 177 Double_t  AliESDV0MI::GetMinimaxSigmaD0(){
 178   //
 179   // calculate mini-max sigma of dca resolution
 180   //
 181   //compv0->fTree->SetAlias("SigmaD2","max(min((SigmaD0+SigmaDE0)*0.5,1.5*SigmaDE0),0.5*SigmaDE0)");
 182   Double_t    effectiveSigma = GetEffectiveSigmaD0();
 183   Double_t    sigmaMMD0 = 0.5*(GetSigmaD0()+effectiveSigma);
 184   sigmaMMD0  = TMath::Min(sigmaMMD0, fgkParams.fPMaxFractionD0*effectiveSigma);
 185   sigmaMMD0  = TMath::Max(sigmaMMD0, fgkParams.fPMinFractionD0*effectiveSigma+fgkParams.fPMinD0);
 186   return sigmaMMD0;
 187 }
 188
 189
 190 Double_t AliESDV0MI::GetLikelihoodAP(Int_t mode0, Int_t mode1){
 191   //
 192   // get likelihood for point angle
 193   //
 194   Double_t sigmaAP = 0.007;            //default sigma
 195   switch (mode0){
 196   case 0:
 197     sigmaAP = GetSigmaAP0();           // mode 0  - covariance matrix estimates used
 198     break;
 199   case 1:
 200     sigmaAP = GetEffectiveSigmaAP0();  // mode 1 - effective sigma used
 201     break;
 202   case 2:
 203     sigmaAP = GetMinimaxSigmaAP0();    // mode 2 - minimax sigma
 204     break;
 205   }
 206   Double_t apNorm = TMath::Min(TMath::ACos(fPointAngle)/sigmaAP,50.);
 207   //normalized point angle, restricted - because of overflow problems in Exp
 208   Double_t likelihood = 0;
 209   switch(mode1){
 210   case 0:
 211     likelihood = TMath::Exp(-0.5*apNorm*apNorm);
 212     // one component
 213     break;
 214   case 1:
 215     likelihood = (TMath::Exp(-0.5*apNorm*apNorm)+0.5* TMath::Exp(-0.25*apNorm*apNorm))/1.5;
 216     // two components
 217     break;
 218   case 2:
 219     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;
 220     // three components
 221     break;
 222   }
 223   return likelihood;
 224 }
 225
 226 Double_t AliESDV0MI::GetLikelihoodD(Int_t mode0, Int_t mode1){
 227   //
 228   // get likelihood for DCA
 229   //
 230   Double_t sigmaD = 0.03;            //default sigma
 231   switch (mode0){
 232   case 0:
 233     sigmaD = GetSigmaD0();           // mode 0  - covariance matrix estimates used
 234     break;
 235   case 1:
 236     sigmaD = GetEffectiveSigmaD0();  // mode 1 - effective sigma used
 237     break;
 238   case 2:
 239     sigmaD = GetMinimaxSigmaD0();    // mode 2 - minimax sigma
 240     break;
 241   }
 242   Double_t dNorm = TMath::Min(fDist2/sigmaD,50.);
 243   //normalized point angle, restricted - because of overflow problems in Exp
 244   Double_t likelihood = 0;
 245   switch(mode1){
 246   case 0:
 247     likelihood = TMath::Exp(-2.*dNorm);
 248     // one component
 249     break;
 250   case 1:
 251     likelihood = (TMath::Exp(-2.*dNorm)+0.5* TMath::Exp(-dNorm))/1.5;
 252     // two components
 253     break;
 254   case 2:
 255     likelihood = (TMath::Exp(-2.*dNorm)+0.5* TMath::Exp(-dNorm)+0.25*TMath::Exp(-0.5*dNorm))/1.75;
 256     // three components
 257     break;
 258   }
 259   return likelihood;
 260
 261 }
 262
 263 Double_t AliESDV0MI::GetLikelihoodC(Int_t mode0, Int_t /*mode1*/){
 264   //
 265   // get likelihood for Causality
 266   // !!!  Causality variables defined in AliITStrackerMI !!!
 267   //      when more information was available
 268   //
 269   Double_t likelihood = 0.5;
 270   Double_t minCausal  = TMath::Min(fCausality[0],fCausality[1]);
 271   Double_t maxCausal  = TMath::Max(fCausality[0],fCausality[1]);
 272   //  minCausal           = TMath::Max(minCausal,0.5*maxCausal);
 273   //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");
 274
 275   switch(mode0){
 276   case 0:
 277     //normalization
 278     likelihood = TMath::Power((1.05-2*(0.8-TMath::Exp(-maxCausal))),4.);
 279     break;
 280   case 1:
 281     likelihood = TMath::Power(1.05-(2*(0.8-TMath::Exp(-maxCausal))+(2*(0.8-TMath::Exp(-minCausal))))*0.5,4.);
 282     break;
 283   }
 284   return likelihood;
 285
 286 }
 287
 288 void AliESDV0MI::SetCausality(Float_t pb0, Float_t pb1, Float_t pa0, Float_t pa1)
 289 {
 290   //
 291   // set probabilities
 292   //
 293   fCausality[0] = pb0;     // probability - track 0 exist before vertex
 294   fCausality[1] = pb1;     // probability - track 1 exist before vertex
 295   fCausality[2] = pa0;     // probability - track 0 exist close after vertex
 296   fCausality[3] = pa1;     // probability - track 1 exist close after vertex
 297 }
 298 void  AliESDV0MI::SetClusters(Int_t *clp, Int_t *clm)
 299 {
 300   //
 301   // Set its clusters indexes
 302   //
 303   for (Int_t i=0;i<6;i++) fClusters[0][i] = clp[i];
 304   for (Int_t i=0;i<6;i++) fClusters[1][i] = clm[i];
 305 }
 306
 307
 308 void AliESDV0MI::SetP(const AliExternalTrackParam & paramp)  {
 309   //
 310   // set track +
 311   //
 312   fParamP   = paramp;
 313 }
 314
 315 void AliESDV0MI::SetM(const AliExternalTrackParam & paramm){
 316   //
 317   //set track -
 318   //
 319   fParamM = paramm;
 320 }
 321
 322 void AliESDV0MI::SetRp(const Double_t *rp){
 323   //
 324   // set pid +
 325   //
 326   for (Int_t i=0;i<5;i++) fRP[i]=rp[i];
 327 }
 328
 329 void AliESDV0MI::SetRm(const Double_t *rm){
 330   //
 331   // set pid -
 332   //
 333   for (Int_t i=0;i<5;i++) fRM[i]=rm[i];
 334 }
 335
 336
 337 void  AliESDV0MI::UpdatePID(Double_t pidp[5], Double_t pidm[5])
 338 {
 339   //
 340   // set PID hypothesy
 341   //
 342   // norm PID to 1
 343   Float_t sump =0;
 344   Float_t summ =0;
 345   for (Int_t i=0;i<5;i++){
 346     fRP[i]=pidp[i];
 347     sump+=fRP[i];
 348     fRM[i]=pidm[i];
 349     summ+=fRM[i];
 350   }
 351   for (Int_t i=0;i<5;i++){
 352     fRP[i]/=sump;
 353     fRM[i]/=summ;
 354   }
 355 }
 356
 357 Float_t AliESDV0MI::GetProb(UInt_t p1, UInt_t p2){
 358   //
 359   //
 360   //
 361   //
 362   return TMath::Max(fRP[p1]+fRM[p2], fRP[p2]+fRM[p1]);
 363 }
 364
 365 Float_t AliESDV0MI::GetEffMass(UInt_t p1, UInt_t p2){
 366   //
 367   // calculate effective mass
 368   //
 369   const Float_t kpmass[5] = {5.10000000000000037e-04,1.05660000000000004e-01,1.39570000000000000e-01,
 370                       4.93599999999999983e-01, 9.38270000000000048e-01};
 371   if (p1>4) return -1;
 372   if (p2>4) return -1;
 373   Float_t mass1 = kpmass[p1];
 374   Float_t mass2 = kpmass[p2];
 375   Double_t *m1 = fPP;
 376   Double_t *m2 = fPM;
 377   //
 378   //if (fRP[p1]+fRM[p2]<fRP[p2]+fRM[p1]){
 379   //  m1 = fPM;
 380   //  m2 = fPP;
 381   //}
 382   //
 383   Float_t e1    = TMath::Sqrt(mass1*mass1+
 384                               m1[0]*m1[0]+
 385                               m1[1]*m1[1]+
 386                               m1[2]*m1[2]);
 387   Float_t e2    = TMath::Sqrt(mass2*mass2+
 388                               m2[0]*m2[0]+
 389                               m2[1]*m2[1]+
 390                               m2[2]*m2[2]);
 391   Float_t mass =
 392     (m2[0]+m1[0])*(m2[0]+m1[0])+
 393     (m2[1]+m1[1])*(m2[1]+m1[1])+
 394     (m2[2]+m1[2])*(m2[2]+m1[2]);
 395
 396   mass = TMath::Sqrt((e1+e2)*(e1+e2)-mass);
 397   return mass;
 398 }
 399
 400 void  AliESDV0MI::Update(Float_t vertex[3])
 401 {
 402   //
 403   // updates Kink Info
 404   //
 405   //  Float_t distance1,distance2;
 406   Float_t distance2;
 407   //
 408   AliHelix phelix(fParamP);
 409   AliHelix mhelix(fParamM);
 410   //
 411   //find intersection linear
 412   //
 413   Double_t phase[2][2],radius[2];
 414   Int_t  points = phelix.GetRPHIintersections(mhelix, phase, radius,200);
 415   Double_t delta1=10000,delta2=10000;
 416   /*
 417   if (points<=0) return;
 418   if (points>0){
 419     phelix.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 420     phelix.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 421     phelix.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 422   }
 423   if (points==2){
 424     phelix.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 425     phelix.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 426     phelix.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 427   }
 428   distance1 = TMath::Min(delta1,delta2);
 429   */
 430   //
 431   //find intersection parabolic
 432   //
 433   points = phelix.GetRPHIintersections(mhelix, phase, radius);
 434   delta1=10000,delta2=10000;
 435   Double_t d1=1000.,d2=10000.;
 436   Double_t err[3],angles[3];
 437   if (points<=0) return;
 438   if (points>0){
 439     phelix.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 440     phelix.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 441     if (TMath::Abs(fParamP.X()-TMath::Sqrt(radius[0])<3) && TMath::Abs(fParamM.X()-TMath::Sqrt(radius[0])<3)){
 442       // if we are close to vertex use error parama
 443       //
 444       err[1] = fParamP.GetCovariance()[0]+fParamM.GetCovariance()[0]+0.05*0.05
 445         +0.3*(fParamP.GetCovariance()[2]+fParamM.GetCovariance()[2]);
 446       err[2] = fParamP.GetCovariance()[2]+fParamM.GetCovariance()[2]+0.05*0.05
 447         +0.3*(fParamP.GetCovariance()[0]+fParamM.GetCovariance()[0]);
 448
 449       phelix.GetAngle(phase[0][0],mhelix,phase[0][1],angles);
 450       Double_t tfi  = TMath::Abs(TMath::Tan(angles[0]));
 451       Double_t tlam = TMath::Abs(TMath::Tan(angles[1]));
 452       err[0] = err[1]/((0.2+tfi)*(0.2+tfi))+err[2]/((0.2+tlam)*(0.2+tlam));
 453       err[0] = ((err[1]*err[2]/((0.2+tfi)*(0.2+tfi)*(0.2+tlam)*(0.2+tlam))))/err[0];
 454       phelix.ParabolicDCA2(mhelix,phase[0][0],phase[0][1],radius[0],delta1,err);
 455     }
 456     Double_t xd[3],xm[3];
 457     phelix.Evaluate(phase[0][0],xd);
 458     mhelix.Evaluate(phase[0][1],xm);
 459     d1 = (xd[0]-xm[0])*(xd[0]-xm[0])+(xd[1]-xm[1])*(xd[1]-xm[1])+(xd[2]-xm[2])*(xd[2]-xm[2]);
 460   }
 461   if (points==2){
 462     phelix.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 463     phelix.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 464     if (TMath::Abs(fParamP.X()-TMath::Sqrt(radius[1])<3) && TMath::Abs(fParamM.X()-TMath::Sqrt(radius[1])<3)){
 465       // if we are close to vertex use error paramatrization
 466       //
 467       err[1] = fParamP.GetCovariance()[0]+fParamM.GetCovariance()[0]+0.05*0.05
 468         +0.3*(fParamP.GetCovariance()[2]+fParamM.GetCovariance()[2]);
 469       err[2] = fParamP.GetCovariance()[2]+fParamM.GetCovariance()[2]+0.05*0.05
 470         +0.3*(fParamP.GetCovariance()[0]+fParamM.GetCovariance()[0]);
 471
 472       phelix.GetAngle(phase[1][0],mhelix,phase[1][1],angles);
 473       Double_t tfi  = TMath::Abs(TMath::Tan(angles[0]));
 474       Double_t tlam = TMath::Abs(TMath::Tan(angles[1]));
 475       err[0] = err[1]/((0.2+tfi)*(0.2+tfi))+err[2]/((0.2+tlam)*(0.2+tlam));
 476       err[0] = ((err[1]*err[2]/((0.2+tfi)*(0.2+tfi)*(0.2+tlam)*(0.2+tlam))))/err[0];
 477       phelix.ParabolicDCA2(mhelix,phase[1][0],phase[1][1],radius[1],delta2,err);
 478     }
 479     Double_t xd[3],xm[3];
 480     phelix.Evaluate(phase[1][0],xd);
 481     mhelix.Evaluate(phase[1][1],xm);
 482     d2 = (xd[0]-xm[0])*(xd[0]-xm[0])+(xd[1]-xm[1])*(xd[1]-xm[1])+(xd[2]-xm[2])*(xd[2]-xm[2]);
 483   }
 484   //
 485   distance2 = TMath::Min(delta1,delta2);
 486   if (delta1<delta2){
 487     //get V0 info
 488     Double_t xd[3],xm[3];
 489     phelix.Evaluate(phase[0][0],xd);
 490     mhelix.Evaluate(phase[0][1], xm);
 491     fXr[0] = 0.5*(xd[0]+xm[0]);
 492     fXr[1] = 0.5*(xd[1]+xm[1]);
 493     fXr[2] = 0.5*(xd[2]+xm[2]);
 494
 495     Float_t wy = fParamP.GetCovariance()[0]/(fParamP.GetCovariance()[0]+fParamM.GetCovariance()[0]);
 496     Float_t wz = fParamP.GetCovariance()[2]/(fParamP.GetCovariance()[2]+fParamM.GetCovariance()[2]);
 497     fXr[0] = 0.5*( (1.-wy)*xd[0]+ wy*xm[0] + (1.-wz)*xd[0]+ wz*xm[0] );
 498     fXr[1] = (1.-wy)*xd[1]+ wy*xm[1];
 499     fXr[2] = (1.-wz)*xd[2]+ wz*xm[2];
 500     //
 501     phelix.GetMomentum(phase[0][0],fPP);
 502     mhelix.GetMomentum(phase[0][1],fPM);
 503     phelix.GetAngle(phase[0][0],mhelix,phase[0][1],fAngle);
 504     fRr = TMath::Sqrt(fXr[0]*fXr[0]+fXr[1]*fXr[1]);
 505   }
 506   else{
 507     Double_t xd[3],xm[3];
 508     phelix.Evaluate(phase[1][0],xd);
 509     mhelix.Evaluate(phase[1][1], xm);
 510     fXr[0] = 0.5*(xd[0]+xm[0]);
 511     fXr[1] = 0.5*(xd[1]+xm[1]);
 512     fXr[2] = 0.5*(xd[2]+xm[2]);
 513     Float_t wy = fParamP.GetCovariance()[0]/(fParamP.GetCovariance()[0]+fParamM.GetCovariance()[0]);
 514     Float_t wz = fParamP.GetCovariance()[2]/(fParamP.GetCovariance()[2]+fParamM.GetCovariance()[2]);
 515     fXr[0] = 0.5*( (1.-wy)*xd[0]+ wy*xm[0] + (1.-wz)*xd[0]+ wz*xm[0] );
 516     fXr[1] = (1.-wy)*xd[1]+ wy*xm[1];
 517     fXr[2] = (1.-wz)*xd[2]+ wz*xm[2];
 518     //
 519     phelix.GetMomentum(phase[1][0], fPP);
 520     mhelix.GetMomentum(phase[1][1], fPM);
 521     phelix.GetAngle(phase[1][0],mhelix,phase[1][1],fAngle);
 522     fRr = TMath::Sqrt(fXr[0]*fXr[0]+fXr[1]*fXr[1]);
 523   }
 524   fDist1 = TMath::Sqrt(TMath::Min(d1,d2));
 525   fDist2 = TMath::Sqrt(distance2);
 526   //
 527   //
 528   Double_t v[3] = {fXr[0]-vertex[0],fXr[1]-vertex[1],fXr[2]-vertex[2]};
 529   Double_t p[3] = {fPP[0]+fPM[0], fPP[1]+fPM[1],fPP[2]+fPM[2]};
 530   Double_t vnorm2 = v[0]*v[0]+v[1]*v[1];
 531   if (TMath::Abs(v[2])>100000) return;
 532   Double_t vnorm3 = TMath::Sqrt(TMath::Abs(v[2]*v[2]+vnorm2));
 533   vnorm2 = TMath::Sqrt(vnorm2);
 534   Double_t pnorm2 = p[0]*p[0]+p[1]*p[1];
 535   Double_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
 536   pnorm2 = TMath::Sqrt(pnorm2);
 537   fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
 538   fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);
 539   fPointAngle   = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
 540   //
 541 }
 542