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