df825985608c1a1d214c56422afb6ee4d189b25c
[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   //
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 kpmass[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 = kpmass[p1]; 
95   Float_t mass2 = kpmass[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