New version of the V0 finder (M.Ivanov)
[u/mrichter/AliRoot.git] / 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 }
56
57
58 void AliESDV0MI::SetCausality(Float_t pb0, Float_t pb1, Float_t pa0, Float_t pa1)
59 {
60   //
61   // set probabilities
62   //
63   fCausality[0] = pb0;     // probability - track 0 exist before vertex
64   fCausality[1] = pb1;     // probability - track 1 exist before vertex
65   fCausality[2] = pa0;     // probability - track 0 exist close after vertex
66   fCausality[3] = pa1;     // probability - track 1 exist close after vertex
67 }
68
69 void AliESDV0MI::SetP(const AliExternalTrackParam & paramp)  {
70   //
71   // set track +
72   //
73   fParamP   = paramp;
74 }
75
76 void AliESDV0MI::SetM(const AliExternalTrackParam & paramm){
77   //
78   //set track -
79   //
80   fParamM = paramm;
81 }
82   
83 void AliESDV0MI::SetRp(const Double_t *rp){
84   //
85   // set pid +
86   //
87   for (Int_t i=0;i<5;i++) fRP[i]=rp[i];
88 }
89
90 void AliESDV0MI::SetRm(const Double_t *rm){
91   //
92   // set pid -
93   //
94   for (Int_t i=0;i<5;i++) fRM[i]=rm[i];
95 }
96
97
98 void  AliESDV0MI::UpdatePID(Double_t pidp[5], Double_t pidm[5])
99 {
100   //
101   // set PID hypothesy
102   //
103   // norm PID to 1
104   Float_t sump =0;
105   Float_t summ =0;
106   for (Int_t i=0;i<5;i++){
107     fRP[i]=pidp[i];
108     sump+=fRP[i];
109     fRM[i]=pidm[i];
110     summ+=fRM[i];
111   }
112   for (Int_t i=0;i<5;i++){
113     fRP[i]/=sump;
114     fRM[i]/=summ;
115   }
116 }
117
118 Float_t AliESDV0MI::GetProb(UInt_t p1, UInt_t p2){
119   //
120   //
121   //
122   //
123   return TMath::Max(fRP[p1]+fRM[p2], fRP[p2]+fRM[p1]);
124 }
125
126 Float_t AliESDV0MI::GetEffMass(UInt_t p1, UInt_t p2){
127   //
128   // calculate effective mass
129   //
130   const Float_t kpmass[5] = {5.10000000000000037e-04,1.05660000000000004e-01,1.39570000000000000e-01,
131                       4.93599999999999983e-01, 9.38270000000000048e-01};
132   if (p1>4) return -1;
133   if (p2>4) return -1;
134   Float_t mass1 = kpmass[p1]; 
135   Float_t mass2 = kpmass[p2];   
136   Double_t *m1 = fPP;
137   Double_t *m2 = fPM;
138   //
139   if (fRP[p1]+fRM[p2]<fRP[p2]+fRM[p1]){
140     m1 = fPM;
141     m2 = fPP;
142   }
143   //
144   Float_t e1    = TMath::Sqrt(mass1*mass1+
145                               m1[0]*m1[0]+
146                               m1[1]*m1[1]+
147                               m1[2]*m1[2]);
148   Float_t e2    = TMath::Sqrt(mass2*mass2+
149                               m2[0]*m2[0]+
150                               m2[1]*m2[1]+
151                               m2[2]*m2[2]);  
152   Float_t mass =  
153     (m2[0]+m1[0])*(m2[0]+m1[0])+
154     (m2[1]+m1[1])*(m2[1]+m1[1])+
155     (m2[2]+m1[2])*(m2[2]+m1[2]);
156   
157   mass = TMath::Sqrt((e1+e2)*(e1+e2)-mass);
158   return mass;
159 }
160
161 void  AliESDV0MI::Update(Float_t vertex[3])
162 {
163   //
164   // updates Kink Info
165   //
166   //  Float_t distance1,distance2;
167   Float_t distance2;
168   //
169   AliHelix phelix(fParamP);
170   AliHelix mhelix(fParamM);    
171   //
172   //find intersection linear
173   //
174   Double_t phase[2][2],radius[2];
175   Int_t  points = phelix.GetRPHIintersections(mhelix, phase, radius,200);
176   Double_t delta1=10000,delta2=10000;  
177   /*
178   if (points<=0) return;
179   if (points>0){
180     phelix.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
181     phelix.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
182     phelix.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
183   }
184   if (points==2){    
185     phelix.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
186     phelix.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
187     phelix.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
188   }
189   distance1 = TMath::Min(delta1,delta2);
190   */
191   //
192   //find intersection parabolic
193   //
194   points = phelix.GetRPHIintersections(mhelix, phase, radius);
195   delta1=10000,delta2=10000;  
196   Double_t d1=1000.,d2=10000.;
197   if (points<=0) return;
198   if (points>0){
199     phelix.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
200     phelix.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
201     Double_t xd[3],xm[3];
202     phelix.Evaluate(phase[0][0],xd);
203     mhelix.Evaluate(phase[0][1],xm);
204     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]);
205   }
206   if (points==2){    
207     phelix.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
208     phelix.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
209     Double_t xd[3],xm[3];
210     phelix.Evaluate(phase[1][0],xd);
211     mhelix.Evaluate(phase[1][1],xm);
212     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]);
213   }
214   //
215   distance2 = TMath::Min(delta1,delta2);
216   if (delta1<delta2){
217     //get V0 info
218     Double_t xd[3],xm[3];
219     phelix.Evaluate(phase[0][0],xd);
220     mhelix.Evaluate(phase[0][1], xm);
221     fXr[0] = 0.5*(xd[0]+xm[0]);
222     fXr[1] = 0.5*(xd[1]+xm[1]);
223     fXr[2] = 0.5*(xd[2]+xm[2]);
224     //
225     phelix.GetMomentum(phase[0][0],fPP);
226     mhelix.GetMomentum(phase[0][1],fPM);
227     phelix.GetAngle(phase[0][0],mhelix,phase[0][1],fAngle);
228     fRr = TMath::Sqrt(fXr[0]*fXr[0]+fXr[1]*fXr[1]);
229   }
230   else{
231     Double_t xd[3],xm[3];
232     phelix.Evaluate(phase[1][0],xd);
233     mhelix.Evaluate(phase[1][1], xm);
234     fXr[0] = 0.5*(xd[0]+xm[0]);
235     fXr[1] = 0.5*(xd[1]+xm[1]);
236     fXr[2] = 0.5*(xd[2]+xm[2]);
237     //
238     phelix.GetMomentum(phase[1][0], fPP);
239     mhelix.GetMomentum(phase[1][1], fPM);
240     phelix.GetAngle(phase[1][0],mhelix,phase[1][1],fAngle);
241     fRr = TMath::Sqrt(fXr[0]*fXr[0]+fXr[1]*fXr[1]);
242   }
243   fDist1 = TMath::Sqrt(TMath::Min(d1,d2));
244   fDist2 = TMath::Sqrt(distance2);      
245   //            
246   //
247   Double_t v[3] = {fXr[0]-vertex[0],fXr[1]-vertex[1],fXr[2]-vertex[2]};
248   Double_t p[3] = {fPP[0]+fPM[0], fPP[1]+fPM[1],fPP[2]+fPM[2]};
249   Double_t vnorm2 = v[0]*v[0]+v[1]*v[1];
250   Double_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
251   vnorm2 = TMath::Sqrt(vnorm2);
252   Double_t pnorm2 = p[0]*p[0]+p[1]*p[1];
253   Double_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
254   pnorm2 = TMath::Sqrt(pnorm2);  
255   fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
256   fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);  
257   fPointAngle   = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
258   //
259 }
260