#include "AliHelix.h"
#include "AliKalmanTrack.h"
-#include "AliExternalTrackParam.h"
+#include "AliTracker.h"
#include "TMath.h"
ClassImp(AliHelix)
alpha=t.GetAlpha();
//
//circle parameters
- fHelix[4]=fHelix[4]/t.GetConvConst(); // C
+ //PH Sometimes fP4 and fHelix[4] are very big and the calculation
+ //PH of the Sqrt cannot be done. To be investigated...
+ fHelix[4]=fHelix[4]/(-1000/0.299792458/AliTracker::GetBz()); // C
cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
Double_t xc, yc, rc;
rc = 1/fHelix[4];
xc = x-fHelix[2]*rc;
- yc = fHelix[0]+TMath::Sqrt(1-(x-xc)*(x-xc)*fHelix[4]*fHelix[4])/fHelix[4];
+ Double_t dummy = 1-(x-xc)*(x-xc)*fHelix[4]*fHelix[4];
+ if (dummy<0) {
+ AliError(Form("The argument of the Sqrt is %f => set to 0\n",dummy));
+ dummy = 0;
+ }
+ yc = fHelix[0]+TMath::Sqrt(dummy)/fHelix[4];
fHelix[6] = xc*cs - yc*sn;
fHelix[7] = xc*sn + yc*cs;
Double_t alpha,x,cs,sn;
const Double_t *param =t.GetParameter();
for (Int_t i=0;i<5;i++) fHelix[i]=param[i];
- x = t.X();
- alpha=t.Alpha();
+ x = t.GetX();
+ alpha=t.GetAlpha();
//
//circle parameters
- fHelix[4]=fHelix[4]/AliKalmanTrack::GetConvConst(); // C
+ //PH Sometimes fP4 and fHelix[4] are very big and the calculation
+ //PH of the Sqrt cannot be done. To be investigated...
+ fHelix[4]=fHelix[4]/(-1000/0.299792458/AliTracker::GetBz()); // C
cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
Double_t xc, yc, rc;
rc = 1/fHelix[4];
xc = x-fHelix[2]*rc;
- yc = fHelix[0]+TMath::Sqrt(1-(x-xc)*(x-xc)*fHelix[4]*fHelix[4])/fHelix[4];
+ Double_t dummy = 1-(x-xc)*(x-xc)*fHelix[4]*fHelix[4];
+ if (dummy<0) {
+ AliError(Form("The argument of the Sqrt is %f => set to 0\n",dummy));
+ dummy = 0;
+ }
+ yc = fHelix[0]+TMath::Sqrt(dummy)/fHelix[4];
fHelix[6] = xc*cs - yc*sn;
fHelix[7] = xc*sn + yc*cs;
//
Double_t pt = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
if (TMath::Abs(conversion)<0.00000001)
- conversion = AliKalmanTrack::GetConvConst();
+ conversion = -1000/0.299792458/AliTracker::GetBz();
//
//
fHelix[4] = charge/(conversion*pt); // C
}
-void AliHelix::GetMomentum(Double_t phase, Double_t p[4],Double_t conversion)
+void AliHelix::GetMomentum(Double_t phase, Double_t p[4],Double_t conversion, Double_t *xr)
{
// return momentum at given phase
Double_t x[3],g[3],gg[3];
Evaluate(phase,x,g,gg);
- if (TMath::Abs(conversion)<0.0001) conversion = AliKalmanTrack::GetConvConst();
+ if (TMath::Abs(conversion)<0.0001) conversion = -1000/0.299792458/AliTracker::GetBz();
Double_t mt = TMath::Sqrt(g[0]*g[0]+g[1]*g[1]);
p[0] = fHelix[8]*g[0]/(mt*conversion);
p[1] = fHelix[8]*g[1]/(mt*conversion);
p[2] = fHelix[8]*g[2]/(mt*conversion);
+ if (xr){
+ xr[0] = x[0]; xr[1] = x[1]; xr[2] = x[2];
+ }
}
void AliHelix::GetAngle(Double_t t1, AliHelix &h, Double_t t2, Double_t angle[3])
return 1;
}
+Double_t AliHelix::GetPointAngle(AliHelix &h, Double_t phase[2], const Float_t * vertex)
+{
+ //
+ // get point angle bettwen two helixes
+ //
+ Double_t r0[3],p0[4];
+ Double_t r1[3],p1[4];
+ GetMomentum(phase[0],p0,1,r0);
+ h.GetMomentum(phase[1],p1,1,r1);
+ //
+ Double_t r[3] = {(r0[0]+r1[0])*0.5-vertex[0],(r0[1]+r1[1])*0.5-vertex[1],(r0[2]+r1[2])*0.5-vertex[2]};
+ //intersection point - relative to the prim vertex
+ Double_t p[3] = { p0[0]+p1[0], p0[1]+p1[1],p0[2]+p1[2]};
+ // derivation vector
+ Double_t normr = TMath::Sqrt(r[0]*r[0]+r[1]*r[1]+r[2]*r[2]);
+ Double_t normp = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
+ Double_t pointAngle = (r[0]*p[0]+r[1]*p[1]+r[2]*p[2])/(normr*normp);
+ return pointAngle;
+}
+
Double_t AliHelix::GetPhase(Double_t x, Double_t y )
{
// Double_t * c1 = &fHelix[6];
//Double_t * c2 = &(h.fHelix[6]);
// Double_t c1[3] = {fHelix[5],fHelix[0],fHelix[8]};
+
+ // PH initiaziation in case of return
+ phase[0][0]=phase[0][1]=phase[1][0]=phase[1][1]=0;
+ ri[0]=ri[1]=1000000;
+
Double_t c1[3] = {0,0,fHelix[8]};
Double_t c2[3] = {h.fHelix[5]-fHelix[5],h.fHelix[0]-fHelix[0],h.fHelix[8]};
x0[0] = (d+c1[2]-c2[2])*c2[0]/(2*d)+ fHelix[5];
y0[0] = (d+c1[2]-c2[2])*c2[1]/(2*d)+ fHelix[0];
// return 0;
- phase[0][0] = GetPhase(x0[0],y0[0]);
- phase[0][1] = h.GetPhase(x0[0],y0[0]);
- ri[0] = x0[0]*x0[0]+y0[0]*y0[0];
+ phase[1][0] = phase[0][0] = GetPhase(x0[0],y0[0]);
+ phase[1][1] = phase[0][1] = h.GetPhase(x0[0],y0[0]);
+ ri[1] = ri[0] = x0[0]*x0[0]+y0[0]*y0[0];
return 1;
}
if ( (d+c2[2])<c1[2]){
//
Double_t xx = c2[0]+ c2[0]*c2[2]/d+ fHelix[5];
Double_t yy = c2[1]+ c2[1]*c2[2]/d+ fHelix[0];
- phase[0][1] = h.GetPhase(xx,yy);
+ phase[1][1] = phase[0][1] = h.GetPhase(xx,yy);
//
Double_t xx2 = c2[0]*c1[2]/d+ fHelix[5];
Double_t yy2 = c2[1]*c1[2]/d+ fHelix[0];
- phase[0][0] = GetPhase(xx2,yy2);
- ri[0] = xx*xx+yy*yy;
+ phase[1][0] = phase[0][0] = GetPhase(xx2,yy2);
+ ri[1] = ri[0] = xx*xx+yy*yy;
return 1;
}
//
Double_t xx = -c2[0]*c1[2]/d+ fHelix[5];
Double_t yy = -c2[1]*c1[2]/d+ fHelix[0];
- phase[0][1] = GetPhase(xx,yy);
+ phase[1][1] = phase[0][1] = GetPhase(xx,yy);
//
Double_t xx2 = c2[0]- c2[0]*c2[2]/d+ fHelix[5];
Double_t yy2 = c2[1]- c2[1]*c2[2]/d+ fHelix[0];
- phase[0][0] = h.GetPhase(xx2,yy2);
- ri[0] = xx*xx+yy*yy;
+ phase[1][0] = phase[0][0] = h.GetPhase(xx2,yy2);
+ ri[1] = ri[0] = xx*xx+yy*yy;
return 1;
}