Additional protection
[u/mrichter/AliRoot.git] / STEER / AliESDV0MI.cxx
CommitLineData
51ad6848 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
29ClassImp(AliESDV0MI)
30
31AliESDV0MI::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
43void AliESDV0MI::SetP(const AliExternalTrackParam & paramp) {
44 //
45 // set mother
46 //
47 fParamP = paramp;
48}
49
50void AliESDV0MI::SetM(const AliExternalTrackParam & paramm){
51 //
52 //set daughter
53 //
54 fParamM = paramm;
55
56}
57
58void 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
78Float_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
86Float_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
121void 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;
b07a036b 136
137 if (points<=0) return;
51ad6848 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.;
b07a036b 155 if (points<=0) return;
51ad6848 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