]>
Commit | Line | Data |
---|---|---|
d0e92d9a | 1 | /////////////////////////////////////////////////////////////////////////// |
2 | // // | |
3 | // AliFmHelix: a helper helix class // | |
4 | // Includes all the operations and specifications of the helix. Can be // | |
5 | // used to determine path lengths, distance of closest approach etc. // | |
6 | // // | |
7 | /////////////////////////////////////////////////////////////////////////// | |
8 | #include <math.h> | |
9 | #include "AliFmHelix.h" | |
10 | #include "AliFmPhysicalHelix.h" | |
11 | #include "PhysicalConstants.h" | |
12 | #include "SystemOfUnits.h" | |
13 | #ifdef __ROOT__ | |
14 | ClassImpT(AliFmPhysicalHelix,double); | |
15 | #endif | |
16 | AliFmPhysicalHelix::AliFmPhysicalHelix(){} | |
17 | ||
18 | AliFmPhysicalHelix::~AliFmPhysicalHelix() { /* nop */ } | |
19 | ||
20 | AliFmPhysicalHelix::AliFmPhysicalHelix(const AliFmThreeVector<double>& p, | |
21 | const AliFmThreeVector<double>& o, | |
22 | double B, double q) | |
23 | { | |
24 | // Constructor from given parameters | |
25 | fH = (q*B <= 0) ? 1 : -1; | |
26 | if(p.y() == 0 && p.x() == 0) | |
27 | SetPhase((M_PI/4)*(1-2.*fH)); | |
28 | else | |
29 | SetPhase(atan2(p.y(),p.x())-fH*M_PI/2); | |
30 | SetDipAngle(atan2(p.z(),p.perp())); | |
31 | fOrigin = o; | |
32 | ||
33 | #ifndef ST_NO_NAMESPACES | |
34 | { | |
35 | using namespace units; | |
36 | #endif | |
37 | SetCurvature(fabs((c_light*nanosecond/meter*q*B/tesla)/ | |
38 | (abs(p.mag())/GeV*fCosDipAngle)/meter)); | |
39 | #ifndef ST_NO_NAMESPACES | |
40 | } | |
41 | #endif | |
42 | } | |
43 | ||
44 | AliFmPhysicalHelix::AliFmPhysicalHelix(double c, double d, double phase, | |
45 | const AliFmThreeVector<double>& o, int h) | |
46 | : AliFmHelix(c, d, phase, o, h) { /* nop */} | |
47 | ||
48 | ||
49 | AliFmThreeVector<double> AliFmPhysicalHelix::Momentum(double B) const | |
50 | { | |
51 | // momentum for given magnetic field | |
52 | if (fSingularity) | |
53 | return(AliFmThreeVector<double>(0,0,0)); | |
54 | else { | |
55 | #ifndef ST_NO_NAMESPACES | |
56 | { | |
57 | using namespace units; | |
58 | #endif | |
59 | double pt = GeV*fabs(c_light*nanosecond/meter*B/tesla)/(fabs(fCurvature)*meter); | |
60 | ||
61 | return (AliFmThreeVector<double>(pt*cos(fPhase+fH*M_PI/2), // pos part pos field | |
62 | pt*sin(fPhase+fH*M_PI/2), | |
63 | pt*tan(fDipAngle))); | |
64 | #ifndef ST_NO_NAMESPACES | |
65 | } | |
66 | #endif | |
67 | } | |
68 | } | |
69 | ||
70 | AliFmThreeVector<double> AliFmPhysicalHelix::MomentumAt(double S, double B) const | |
71 | { | |
72 | // Obtain phase-shifted momentum from phase-shift of origin | |
73 | double xc = this->XCenter(); | |
74 | double yc = this->YCenter(); | |
75 | double rx = (Y(S)-yc)/(fOrigin.y()-yc); | |
76 | double ry = (X(S)-xc)/(fOrigin.x()-xc); | |
77 | return (this->Momentum(B)).pseudoProduct(rx,ry,1.0); | |
78 | } | |
79 | ||
80 | int AliFmPhysicalHelix::Charge(double B) const | |
81 | { | |
82 | // charge | |
83 | return (B > 0 ? -fH : fH); | |
84 | } | |
85 | ||
86 | double AliFmPhysicalHelix::GeometricSignedDistance(double x, double y) | |
87 | { | |
88 | // Geometric signed distance | |
89 | double thePath = this->PathLength(x,y); | |
90 | AliFmThreeVector<double> tDCA2dPosition = this->At(thePath); | |
91 | tDCA2dPosition.setZ(0); | |
92 | AliFmThreeVector<double> position(x,y,0); | |
93 | AliFmThreeVector<double> tDCAVec = (tDCA2dPosition-position); | |
94 | AliFmThreeVector<double> momVec; | |
95 | // Deal with straight tracks | |
96 | if (this->fSingularity) { | |
97 | momVec = this->At(1)- this->At(0); | |
98 | momVec.setZ(0); | |
99 | } | |
100 | else { | |
101 | momVec = this->MomentumAt(thePath,1./tesla); // Don't care about Bmag. Helicity is what matters. | |
102 | momVec.setZ(0); | |
103 | } | |
104 | ||
105 | double cross = tDCAVec.x()*momVec.y() - tDCAVec.y()*momVec.x(); | |
106 | double theSign = (cross>=0) ? 1. : -1.; | |
107 | return theSign*tDCAVec.perp(); | |
108 | } | |
109 | ||
110 | double AliFmPhysicalHelix::CurvatureSignedDistance(double x, double y) | |
111 | { | |
112 | // Protect against fH = 0 or zero field | |
113 | if (this->fSingularity || abs(this->fH)<=0) { | |
114 | return (this->GeometricSignedDistance(x,y)); | |
115 | } | |
116 | else { | |
117 | return (this->GeometricSignedDistance(x,y))/(this->fH); | |
118 | } | |
119 | ||
120 | } | |
121 | ||
122 | double AliFmPhysicalHelix::GeometricSignedDistance(const AliFmThreeVector<double>& pos) | |
123 | { | |
124 | // Geometric distance | |
125 | double sdca2d = this->GeometricSignedDistance(pos.x(),pos.y()); | |
126 | double theSign = (sdca2d>=0) ? 1. : -1.; | |
127 | return (this->Distance(pos))*theSign; | |
128 | } | |
129 | ||
130 | double AliFmPhysicalHelix::CurvatureSignedDistance(const AliFmThreeVector<double>& pos) | |
131 | { | |
132 | // Distance with the sign dependent on curvature sigm | |
133 | double sdca2d = this->CurvatureSignedDistance(pos.x(),pos.y()); | |
134 | double theSign = (sdca2d>=0) ? 1. : -1.; | |
135 | return (this->Distance(pos))*theSign; | |
136 | } |