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 //    Origin: Marian Ivanov marian.ivanov@cern.ch
  20 //-------------------------------------------------------------------------
  21
  22 #include <Riostream.h>
  23 #include <TMath.h>
  24
  25 #include "AliESDV0MI.h"
  26 #include "AliHelix.h"
  27
  28
  29 ClassImp(AliESDV0MI)
  30
  31 AliESDV0MI::AliESDV0MI() :
  32   AliESDv0(),
  33   fParamP(),
  34   fParamM(),
  35   fID(0),
  36   fDist1(-1),
  37   fDist2(-1),
  38   fRr(-1),
  39   fStatus(0),
  40   fRow0(-1),
  41   fDistNorm(0),
  42   fDistSigma(0),
  43   fChi2Before(0),
  44   fNBefore(0),
  45   fChi2After(0),
  46   fNAfter(0),
  47   fPointAngleFi(0),
  48   fPointAngleTh(0),
  49   fPointAngle(0)
  50 {
  51   //
  52   //Dafault constructor
  53   //
  54   for (Int_t i=0;i<4;i++){fCausality[i]=0;}
  55   for (Int_t i=0;i<6;i++){fClusters[0][i]=0; fClusters[1][i]=0;}
  56   for (Int_t i=0;i<2;i++){fNormDCAPrim[0]=0;fNormDCAPrim[1]=0;}
  57 }
  58
  59
  60 void AliESDV0MI::SetCausality(Float_t pb0, Float_t pb1, Float_t pa0, Float_t pa1)
  61 {
  62   //
  63   // set probabilities
  64   //
  65   fCausality[0] = pb0;     // probability - track 0 exist before vertex
  66   fCausality[1] = pb1;     // probability - track 1 exist before vertex
  67   fCausality[2] = pa0;     // probability - track 0 exist close after vertex
  68   fCausality[3] = pa1;     // probability - track 1 exist close after vertex
  69 }
  70 void  AliESDV0MI::SetClusters(Int_t *clp, Int_t *clm)
  71 {
  72   //
  73   // Set its clusters indexes
  74   //
  75   for (Int_t i=0;i<6;i++) fClusters[0][i] = clp[i];
  76   for (Int_t i=0;i<6;i++) fClusters[1][i] = clm[i];
  77 }
  78
  79
  80 void AliESDV0MI::SetP(const AliExternalTrackParam & paramp)  {
  81   //
  82   // set track +
  83   //
  84   fParamP   = paramp;
  85 }
  86
  87 void AliESDV0MI::SetM(const AliExternalTrackParam & paramm){
  88   //
  89   //set track -
  90   //
  91   fParamM = paramm;
  92 }
  93
  94 void AliESDV0MI::SetRp(const Double_t *rp){
  95   //
  96   // set pid +
  97   //
  98   for (Int_t i=0;i<5;i++) fRP[i]=rp[i];
  99 }
 100
 101 void AliESDV0MI::SetRm(const Double_t *rm){
 102   //
 103   // set pid -
 104   //
 105   for (Int_t i=0;i<5;i++) fRM[i]=rm[i];
 106 }
 107
 108
 109 void  AliESDV0MI::UpdatePID(Double_t pidp[5], Double_t pidm[5])
 110 {
 111   //
 112   // set PID hypothesy
 113   //
 114   // norm PID to 1
 115   Float_t sump =0;
 116   Float_t summ =0;
 117   for (Int_t i=0;i<5;i++){
 118     fRP[i]=pidp[i];
 119     sump+=fRP[i];
 120     fRM[i]=pidm[i];
 121     summ+=fRM[i];
 122   }
 123   for (Int_t i=0;i<5;i++){
 124     fRP[i]/=sump;
 125     fRM[i]/=summ;
 126   }
 127 }
 128
 129 Float_t AliESDV0MI::GetProb(UInt_t p1, UInt_t p2){
 130   //
 131   //
 132   //
 133   //
 134   return TMath::Max(fRP[p1]+fRM[p2], fRP[p2]+fRM[p1]);
 135 }
 136
 137 Float_t AliESDV0MI::GetEffMass(UInt_t p1, UInt_t p2){
 138   //
 139   // calculate effective mass
 140   //
 141   const Float_t kpmass[5] = {5.10000000000000037e-04,1.05660000000000004e-01,1.39570000000000000e-01,
 142                       4.93599999999999983e-01, 9.38270000000000048e-01};
 143   if (p1>4) return -1;
 144   if (p2>4) return -1;
 145   Float_t mass1 = kpmass[p1];
 146   Float_t mass2 = kpmass[p2];
 147   Double_t *m1 = fPP;
 148   Double_t *m2 = fPM;
 149   //
 150   //if (fRP[p1]+fRM[p2]<fRP[p2]+fRM[p1]){
 151   //  m1 = fPM;
 152   //  m2 = fPP;
 153   //}
 154   //
 155   Float_t e1    = TMath::Sqrt(mass1*mass1+
 156                               m1[0]*m1[0]+
 157                               m1[1]*m1[1]+
 158                               m1[2]*m1[2]);
 159   Float_t e2    = TMath::Sqrt(mass2*mass2+
 160                               m2[0]*m2[0]+
 161                               m2[1]*m2[1]+
 162                               m2[2]*m2[2]);
 163   Float_t mass =
 164     (m2[0]+m1[0])*(m2[0]+m1[0])+
 165     (m2[1]+m1[1])*(m2[1]+m1[1])+
 166     (m2[2]+m1[2])*(m2[2]+m1[2]);
 167
 168   mass = TMath::Sqrt((e1+e2)*(e1+e2)-mass);
 169   return mass;
 170 }
 171
 172 void  AliESDV0MI::Update(Float_t vertex[3])
 173 {
 174   //
 175   // updates Kink Info
 176   //
 177   //  Float_t distance1,distance2;
 178   Float_t distance2;
 179   //
 180   AliHelix phelix(fParamP);
 181   AliHelix mhelix(fParamM);
 182   //
 183   //find intersection linear
 184   //
 185   Double_t phase[2][2],radius[2];
 186   Int_t  points = phelix.GetRPHIintersections(mhelix, phase, radius,200);
 187   Double_t delta1=10000,delta2=10000;
 188   /*
 189   if (points<=0) return;
 190   if (points>0){
 191     phelix.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 192     phelix.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 193     phelix.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 194   }
 195   if (points==2){
 196     phelix.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 197     phelix.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 198     phelix.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 199   }
 200   distance1 = TMath::Min(delta1,delta2);
 201   */
 202   //
 203   //find intersection parabolic
 204   //
 205   points = phelix.GetRPHIintersections(mhelix, phase, radius);
 206   delta1=10000,delta2=10000;
 207   Double_t d1=1000.,d2=10000.;
 208   Double_t err[3],angles[3];
 209   if (points<=0) return;
 210   if (points>0){
 211     phelix.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 212     phelix.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 213     if (TMath::Abs(fParamP.X()-TMath::Sqrt(radius[0])<3) && TMath::Abs(fParamM.X()-TMath::Sqrt(radius[0])<3)){
 214       // if we are close to vertex use error parama
 215       //
 216       err[1] = fParamP.GetCovariance()[0]+fParamM.GetCovariance()[0]+0.05*0.05
 217         +0.3*(fParamP.GetCovariance()[2]+fParamM.GetCovariance()[2]);
 218       err[2] = fParamP.GetCovariance()[2]+fParamM.GetCovariance()[2]+0.05*0.05
 219         +0.3*(fParamP.GetCovariance()[0]+fParamM.GetCovariance()[0]);
 220
 221       phelix.GetAngle(phase[0][0],mhelix,phase[0][1],angles);
 222       Double_t tfi  = TMath::Abs(TMath::Tan(angles[0]));
 223       Double_t tlam = TMath::Abs(TMath::Tan(angles[1]));
 224       err[0] = err[1]/((0.2+tfi)*(0.2+tfi))+err[2]/((0.2+tlam)*(0.2+tlam));
 225       err[0] = ((err[1]*err[2]/((0.2+tfi)*(0.2+tfi)*(0.2+tlam)*(0.2+tlam))))/err[0];
 226       phelix.ParabolicDCA2(mhelix,phase[0][0],phase[0][1],radius[0],delta1,err);
 227     }
 228     Double_t xd[3],xm[3];
 229     phelix.Evaluate(phase[0][0],xd);
 230     mhelix.Evaluate(phase[0][1],xm);
 231     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]);
 232   }
 233   if (points==2){
 234     phelix.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 235     phelix.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 236     if (TMath::Abs(fParamP.X()-TMath::Sqrt(radius[1])<3) && TMath::Abs(fParamM.X()-TMath::Sqrt(radius[1])<3)){
 237       // if we are close to vertex use error paramatrization
 238       //
 239       err[1] = fParamP.GetCovariance()[0]+fParamM.GetCovariance()[0]+0.05*0.05
 240         +0.3*(fParamP.GetCovariance()[2]+fParamM.GetCovariance()[2]);
 241       err[2] = fParamP.GetCovariance()[2]+fParamM.GetCovariance()[2]+0.05*0.05
 242         +0.3*(fParamP.GetCovariance()[0]+fParamM.GetCovariance()[0]);
 243
 244       phelix.GetAngle(phase[1][0],mhelix,phase[1][1],angles);
 245       Double_t tfi  = TMath::Abs(TMath::Tan(angles[0]));
 246       Double_t tlam = TMath::Abs(TMath::Tan(angles[1]));
 247       err[0] = err[1]/((0.2+tfi)*(0.2+tfi))+err[2]/((0.2+tlam)*(0.2+tlam));
 248       err[0] = ((err[1]*err[2]/((0.2+tfi)*(0.2+tfi)*(0.2+tlam)*(0.2+tlam))))/err[0];
 249       phelix.ParabolicDCA2(mhelix,phase[1][0],phase[1][1],radius[1],delta2,err);
 250     }
 251     Double_t xd[3],xm[3];
 252     phelix.Evaluate(phase[1][0],xd);
 253     mhelix.Evaluate(phase[1][1],xm);
 254     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]);
 255   }
 256   //
 257   distance2 = TMath::Min(delta1,delta2);
 258   if (delta1<delta2){
 259     //get V0 info
 260     Double_t xd[3],xm[3];
 261     phelix.Evaluate(phase[0][0],xd);
 262     mhelix.Evaluate(phase[0][1], xm);
 263     fXr[0] = 0.5*(xd[0]+xm[0]);
 264     fXr[1] = 0.5*(xd[1]+xm[1]);
 265     fXr[2] = 0.5*(xd[2]+xm[2]);
 266
 267     Float_t wy = fParamP.GetCovariance()[0]/(fParamP.GetCovariance()[0]+fParamM.GetCovariance()[0]);
 268     Float_t wz = fParamP.GetCovariance()[2]/(fParamP.GetCovariance()[2]+fParamM.GetCovariance()[2]);
 269     fXr[0] = 0.5*( (1.-wy)*xd[0]+ wy*xm[0] + (1.-wz)*xd[0]+ wz*xm[0] );
 270     fXr[1] = (1.-wy)*xd[1]+ wy*xm[1];
 271     fXr[2] = (1.-wz)*xd[2]+ wz*xm[2];
 272     //
 273     phelix.GetMomentum(phase[0][0],fPP);
 274     mhelix.GetMomentum(phase[0][1],fPM);
 275     phelix.GetAngle(phase[0][0],mhelix,phase[0][1],fAngle);
 276     fRr = TMath::Sqrt(fXr[0]*fXr[0]+fXr[1]*fXr[1]);
 277   }
 278   else{
 279     Double_t xd[3],xm[3];
 280     phelix.Evaluate(phase[1][0],xd);
 281     mhelix.Evaluate(phase[1][1], xm);
 282     fXr[0] = 0.5*(xd[0]+xm[0]);
 283     fXr[1] = 0.5*(xd[1]+xm[1]);
 284     fXr[2] = 0.5*(xd[2]+xm[2]);
 285     Float_t wy = fParamP.GetCovariance()[0]/(fParamP.GetCovariance()[0]+fParamM.GetCovariance()[0]);
 286     Float_t wz = fParamP.GetCovariance()[2]/(fParamP.GetCovariance()[2]+fParamM.GetCovariance()[2]);
 287     fXr[0] = 0.5*( (1.-wy)*xd[0]+ wy*xm[0] + (1.-wz)*xd[0]+ wz*xm[0] );
 288     fXr[1] = (1.-wy)*xd[1]+ wy*xm[1];
 289     fXr[2] = (1.-wz)*xd[2]+ wz*xm[2];
 290     //
 291     phelix.GetMomentum(phase[1][0], fPP);
 292     mhelix.GetMomentum(phase[1][1], fPM);
 293     phelix.GetAngle(phase[1][0],mhelix,phase[1][1],fAngle);
 294     fRr = TMath::Sqrt(fXr[0]*fXr[0]+fXr[1]*fXr[1]);
 295   }
 296   fDist1 = TMath::Sqrt(TMath::Min(d1,d2));
 297   fDist2 = TMath::Sqrt(distance2);
 298   //
 299   //
 300   Double_t v[3] = {fXr[0]-vertex[0],fXr[1]-vertex[1],fXr[2]-vertex[2]};
 301   Double_t p[3] = {fPP[0]+fPM[0], fPP[1]+fPM[1],fPP[2]+fPM[2]};
 302   Double_t vnorm2 = v[0]*v[0]+v[1]*v[1];
 303   Double_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
 304   vnorm2 = TMath::Sqrt(vnorm2);
 305   Double_t pnorm2 = p[0]*p[0]+p[1]*p[1];
 306   Double_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
 307   pnorm2 = TMath::Sqrt(pnorm2);
 308   fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
 309   fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);
 310   fPointAngle   = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
 311   //
 312 }
 313