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 #include <TPDGCode.h>
  25 #include "AliESDV0MI.h"
  26 #include "AliHelix.h"
  27
  28
  29 ClassImp(AliESDV0MI)
  30
  31 AliESDV0MI::AliESDV0MI(){
  32   //
  33   //Dafault constructor
  34   //
  35   fID =0;
  36   fDist1  =-1;
  37   fDist2  =-1;
  38   fRr     =-1;
  39   fStatus = 0;
  40   fRow0   =-1;
  41 }
  42
  43 void AliESDV0MI::SetP(const AliExternalTrackParam & paramp)  {
  44   //
  45   // set mother
  46   //
  47   fParamP   = paramp;
  48 }
  49
  50 void AliESDV0MI::SetM(const AliExternalTrackParam & paramm){
  51   //
  52   //set daughter
  53   //
  54   fParamM = paramm;
  55
  56 }
  57
  58 void  AliESDV0MI::UpdatePID(Double_t pidp[5], Double_t pidm[5])
  59 {
  60   //
  61   // set PID hypothesy
  62   //
  63   // norm PID to 1
  64   Float_t sump =0;
  65   Float_t summ =0;
  66   for (Int_t i=0;i<5;i++){
  67     fRP[i]=pidp[i];
  68     sump+=fRP[i];
  69     fRM[i]=pidm[i];
  70     summ+=fRM[i];
  71   }
  72   for (Int_t i=0;i<5;i++){
  73     fRP[i]/=sump;
  74     fRM[i]/=summ;
  75   }
  76 }
  77
  78 Float_t AliESDV0MI::GetProb(UInt_t p1, UInt_t p2){
  79   //
  80   //
  81   //
  82   //
  83   return TMath::Max(fRP[p1]+fRM[p2], fRP[p2]+fRM[p1]);
  84 }
  85
  86 Float_t AliESDV0MI::GetEffMass(UInt_t p1, UInt_t p2){
  87   //
  88   // calculate effective mass
  89   //
  90   const Float_t pmass[5] = {5.10000000000000037e-04,1.05660000000000004e-01,1.39570000000000000e-01,
  91                       4.93599999999999983e-01, 9.38270000000000048e-01};
  92   if (p1>4) return -1;
  93   if (p2>4) return -1;
  94   Float_t mass1 = pmass[p1];
  95   Float_t mass2 = pmass[p2];
  96   Double_t *m1 = fPP;
  97   Double_t *m2 = fPM;
  98   //
  99   if (fRP[p1]+fRM[p2]<fRP[p2]+fRM[p1]){
 100     m1 = fPM;
 101     m2 = fPP;
 102   }
 103   //
 104   Float_t e1    = TMath::Sqrt(mass1*mass1+
 105                               m1[0]*m1[0]+
 106                               m1[1]*m1[1]+
 107                               m1[2]*m1[2]);
 108   Float_t e2    = TMath::Sqrt(mass2*mass2+
 109                               m2[0]*m2[0]+
 110                               m2[1]*m2[1]+
 111                               m2[2]*m2[2]);
 112   Float_t mass =
 113     (m2[0]+m1[0])*(m2[0]+m1[0])+
 114     (m2[1]+m1[1])*(m2[1]+m1[1])+
 115     (m2[2]+m1[2])*(m2[2]+m1[2]);
 116
 117   mass = TMath::Sqrt((e1+e2)*(e1+e2)-mass);
 118   return mass;
 119 }
 120
 121 void  AliESDV0MI::Update(Float_t vertex[3])
 122 {
 123   //
 124   // updates Kink Info
 125   //
 126   Float_t distance1,distance2;
 127   //
 128   AliHelix phelix(fParamP);
 129   AliHelix mhelix(fParamM);
 130   //
 131   //find intersection linear
 132   //
 133   Double_t phase[2][2],radius[2];
 134   Int_t  points = phelix.GetRPHIintersections(mhelix, phase, radius,200);
 135   Double_t delta1=10000,delta2=10000;
 136
 137   if (points<=0) return;
 138   if (points>0){
 139     phelix.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 140     phelix.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 141     phelix.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 142   }
 143   if (points==2){
 144     phelix.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 145     phelix.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 146     phelix.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 147   }
 148   distance1 = TMath::Min(delta1,delta2);
 149   //
 150   //find intersection parabolic
 151   //
 152   points = phelix.GetRPHIintersections(mhelix, phase, radius);
 153   delta1=10000,delta2=10000;
 154   Double_t d1=1000.,d2=10000.;
 155   if (points<=0) return;
 156   if (points>0){
 157     phelix.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 158     phelix.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 159     Double_t xd[3],xm[3];
 160     phelix.Evaluate(phase[0][0],xd);
 161     mhelix.Evaluate(phase[0][1],xm);
 162     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]);
 163   }
 164   if (points==2){
 165     phelix.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 166     phelix.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 167     Double_t xd[3],xm[3];
 168     phelix.Evaluate(phase[1][0],xd);
 169     mhelix.Evaluate(phase[1][1],xm);
 170     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]);
 171   }
 172   //
 173   distance2 = TMath::Min(delta1,delta2);
 174   if (delta1<delta2){
 175     //get V0 info
 176     Double_t xd[3],xm[3];
 177     phelix.Evaluate(phase[0][0],xd);
 178     mhelix.Evaluate(phase[0][1], xm);
 179     fXr[0] = 0.5*(xd[0]+xm[0]);
 180     fXr[1] = 0.5*(xd[1]+xm[1]);
 181     fXr[2] = 0.5*(xd[2]+xm[2]);
 182     //
 183     phelix.GetMomentum(phase[0][0],fPP);
 184     mhelix.GetMomentum(phase[0][1],fPM);
 185     phelix.GetAngle(phase[0][0],mhelix,phase[0][1],fAngle);
 186     fRr = TMath::Sqrt(fXr[0]*fXr[0]+fXr[1]*fXr[1]);
 187   }
 188   else{
 189     Double_t xd[3],xm[3];
 190     phelix.Evaluate(phase[1][0],xd);
 191     mhelix.Evaluate(phase[1][1], xm);
 192     fXr[0] = 0.5*(xd[0]+xm[0]);
 193     fXr[1] = 0.5*(xd[1]+xm[1]);
 194     fXr[2] = 0.5*(xd[2]+xm[2]);
 195     //
 196     phelix.GetMomentum(phase[1][0], fPP);
 197     mhelix.GetMomentum(phase[1][1], fPM);
 198     phelix.GetAngle(phase[1][0],mhelix,phase[1][1],fAngle);
 199     fRr = TMath::Sqrt(fXr[0]*fXr[0]+fXr[1]*fXr[1]);
 200   }
 201   fDist1 = TMath::Sqrt(TMath::Min(d1,d2));
 202   fDist2 = TMath::Sqrt(distance2);
 203   //
 204   //
 205   Float_t v[3] = {fXr[0]-vertex[0],fXr[1]-vertex[1],fXr[2]-vertex[2]};
 206   Float_t p[3] = {fPP[0]+fPM[0], fPP[1]+fPM[1],fPP[2]+fPM[2]};
 207   Float_t vnorm2 = v[0]*v[0]+v[1]*v[1];
 208   Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
 209   vnorm2 = TMath::Sqrt(vnorm2);
 210   Float_t pnorm2 = p[0]*p[0]+p[1]*p[1];
 211   Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
 212   pnorm2 = TMath::Sqrt(pnorm2);
 213   fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
 214   fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);
 215   fPointAngle   = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
 216   //
 217 }
 218