Merge branch 'master' into TPCdev
[u/mrichter/AliRoot.git] / STEER / STEER / AliV0.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 //-------------------------------------------------------------------------
17 //
18 //    Implementation of the V0 vertex class
19 //    Numerical part - AliHelix functionality used             
20 //
21 //    Origin: Marian Ivanov marian.ivanov@cern.ch
22 //-------------------------------------------------------------------------
23 #include <TMath.h>
24
25 #include "AliV0.h"
26 #include "AliHelix.h"
27
28
29 ClassImp(AliV0)
30
31 void  AliV0::Update(Float_t vertex[3])
32 {
33   //
34   // updates Kink Info
35   //
36   //  Float_t distance1,distance2;
37   Float_t distance2;
38   //
39   AliHelix phelix(fParamP);
40   AliHelix mhelix(fParamN);    
41   //
42   //find intersection linear
43   //
44   Double_t phase[2][2]={{0}},radius[2]={0};
45   Int_t  points = phelix.GetRPHIintersections(mhelix, phase, radius,200);
46   Double_t delta1=10000,delta2=10000;  
47   /*
48   if (points<=0) return;
49   if (points>0){
50     phelix.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
51     phelix.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
52     phelix.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
53   }
54   if (points==2){    
55     phelix.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
56     phelix.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
57     phelix.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
58   }
59   distance1 = TMath::Min(delta1,delta2);
60   */
61   //
62   //find intersection parabolic
63   //
64   points = phelix.GetRPHIintersections(mhelix, phase, radius);
65   delta1=10000,delta2=10000;  
66   Double_t d1=1000.,d2=10000.;
67   Double_t err[3],angles[3];
68   if (points<=0) return;
69   if (points>0){
70     phelix.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
71     phelix.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
72     if (TMath::Abs(fParamP.GetX()-TMath::Sqrt(radius[0])<3) && TMath::Abs(fParamN.GetX()-TMath::Sqrt(radius[0])<3)){
73       // if we are close to vertex use error parama
74       //
75       err[1] = fParamP.GetCovariance()[0]+fParamN.GetCovariance()[0]+0.05*0.05
76         +0.3*(fParamP.GetCovariance()[2]+fParamN.GetCovariance()[2]);
77       err[2] = fParamP.GetCovariance()[2]+fParamN.GetCovariance()[2]+0.05*0.05
78         +0.3*(fParamP.GetCovariance()[0]+fParamN.GetCovariance()[0]);
79       
80       phelix.GetAngle(phase[0][0],mhelix,phase[0][1],angles);
81       Double_t tfi  = TMath::Abs(TMath::Tan(angles[0]));
82       Double_t tlam = TMath::Abs(TMath::Tan(angles[1]));
83       err[0] = err[1]/((0.2+tfi)*(0.2+tfi))+err[2]/((0.2+tlam)*(0.2+tlam));
84       err[0] = ((err[1]*err[2]/((0.2+tfi)*(0.2+tfi)*(0.2+tlam)*(0.2+tlam))))/err[0];
85       phelix.ParabolicDCA2(mhelix,phase[0][0],phase[0][1],radius[0],delta1,err);
86     }
87     Double_t xd[3],xm[3];
88     phelix.Evaluate(phase[0][0],xd);
89     mhelix.Evaluate(phase[0][1],xm);
90     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]);
91   }
92   if (points==2){    
93     phelix.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
94     phelix.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
95     if (TMath::Abs(fParamP.GetX()-TMath::Sqrt(radius[1])<3) && TMath::Abs(fParamN.GetX()-TMath::Sqrt(radius[1])<3)){
96       // if we are close to vertex use error paramatrization
97       //
98       err[1] = fParamP.GetCovariance()[0]+fParamN.GetCovariance()[0]+0.05*0.05
99         +0.3*(fParamP.GetCovariance()[2]+fParamN.GetCovariance()[2]);
100       err[2] = fParamP.GetCovariance()[2]+fParamN.GetCovariance()[2]+0.05*0.05
101         +0.3*(fParamP.GetCovariance()[0]+fParamN.GetCovariance()[0]);
102       
103       phelix.GetAngle(phase[1][0],mhelix,phase[1][1],angles);
104       Double_t tfi  = TMath::Abs(TMath::Tan(angles[0]));
105       Double_t tlam = TMath::Abs(TMath::Tan(angles[1]));
106       err[0] = err[1]/((0.2+tfi)*(0.2+tfi))+err[2]/((0.2+tlam)*(0.2+tlam));     
107       err[0] = ((err[1]*err[2]/((0.2+tfi)*(0.2+tfi)*(0.2+tlam)*(0.2+tlam))))/err[0];
108       phelix.ParabolicDCA2(mhelix,phase[1][0],phase[1][1],radius[1],delta2,err);
109     }
110     Double_t xd[3],xm[3];
111     phelix.Evaluate(phase[1][0],xd);
112     mhelix.Evaluate(phase[1][1],xm);
113     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]);
114   }
115   //
116   distance2 = TMath::Min(delta1,delta2);
117   if (delta1<delta2){
118     //get V0 info
119     Double_t xd[3],xm[3];
120     phelix.Evaluate(phase[0][0],xd);
121     mhelix.Evaluate(phase[0][1], xm);
122     fPos[0] = 0.5*(xd[0]+xm[0]);
123     fPos[1] = 0.5*(xd[1]+xm[1]);
124     fPos[2] = 0.5*(xd[2]+xm[2]);
125
126     Float_t wy = fParamP.GetCovariance()[0]/(fParamP.GetCovariance()[0]+fParamN.GetCovariance()[0]);
127     Float_t wz = fParamP.GetCovariance()[2]/(fParamP.GetCovariance()[2]+fParamN.GetCovariance()[2]);
128     fPos[0] = 0.5*( (1.-wy)*xd[0]+ wy*xm[0] + (1.-wz)*xd[0]+ wz*xm[0] );
129     fPos[1] = (1.-wy)*xd[1]+ wy*xm[1];
130     fPos[2] = (1.-wz)*xd[2]+ wz*xm[2];
131     //
132     phelix.GetMomentum(phase[0][0],fPmom);
133     mhelix.GetMomentum(phase[0][1],fNmom);
134     phelix.GetAngle(phase[0][0],mhelix,phase[0][1],fAngle);
135     fRr = TMath::Sqrt(fPos[0]*fPos[0]+fPos[1]*fPos[1]);
136   }
137   else{
138     Double_t xd[3],xm[3];
139     phelix.Evaluate(phase[1][0],xd);
140     mhelix.Evaluate(phase[1][1], xm);
141     fPos[0] = 0.5*(xd[0]+xm[0]);
142     fPos[1] = 0.5*(xd[1]+xm[1]);
143     fPos[2] = 0.5*(xd[2]+xm[2]);
144     Float_t wy = fParamP.GetCovariance()[0]/(fParamP.GetCovariance()[0]+fParamN.GetCovariance()[0]);
145     Float_t wz = fParamP.GetCovariance()[2]/(fParamP.GetCovariance()[2]+fParamN.GetCovariance()[2]);
146     fPos[0] = 0.5*( (1.-wy)*xd[0]+ wy*xm[0] + (1.-wz)*xd[0]+ wz*xm[0] );
147     fPos[1] = (1.-wy)*xd[1]+ wy*xm[1];
148     fPos[2] = (1.-wz)*xd[2]+ wz*xm[2];
149     //
150     phelix.GetMomentum(phase[1][0], fPmom);
151     mhelix.GetMomentum(phase[1][1], fNmom);
152     phelix.GetAngle(phase[1][0],mhelix,phase[1][1],fAngle);
153     fRr = TMath::Sqrt(fPos[0]*fPos[0]+fPos[1]*fPos[1]);
154   }
155   //Bo:  fDist1 = TMath::Sqrt(TMath::Min(d1,d2));
156   //Bo:  fDist2 = TMath::Sqrt(distance2);
157   fDcaV0Daughters = TMath::Sqrt(distance2);//Bo:
158   //            
159   //
160   Double_t v[3] = {fPos[0]-vertex[0],fPos[1]-vertex[1],fPos[2]-vertex[2]};
161   Double_t p[3] = {fPmom[0]+fNmom[0], fPmom[1]+fNmom[1],fPmom[2]+fNmom[2]};
162   Double_t vnorm2 = v[0]*v[0]+v[1]*v[1];
163   if (TMath::Abs(v[2])>100000) return;
164   Double_t vnorm3 = TMath::Sqrt(TMath::Abs(v[2]*v[2]+vnorm2));
165   vnorm2 = TMath::Sqrt(vnorm2);
166   Double_t pnorm2 = p[0]*p[0]+p[1]*p[1];
167   Double_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
168   pnorm2 = TMath::Sqrt(pnorm2);  
169   if (vnorm2*pnorm2>0) {
170     fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
171   }
172   if (vnorm3*pnorm3>0) {
173     fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);  
174     fPointAngle   = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
175   }
176   //
177 }
178