Lines getting the matched track moved to a method in AliCalorimeterUtils. Lines copie...
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemto / AliFmPhysicalHelix.cxx
CommitLineData
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__
14ClassImpT(AliFmPhysicalHelix,double);
15#endif
16AliFmPhysicalHelix::AliFmPhysicalHelix(){}
17
18AliFmPhysicalHelix::~AliFmPhysicalHelix() { /* nop */ }
19
20AliFmPhysicalHelix::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);
69c1c8ff 30 SetDipAngle(atan2(p.z(),p.Perp()));
d0e92d9a 31 fOrigin = o;
32
33#ifndef ST_NO_NAMESPACES
34 {
35 using namespace units;
36#endif
81a888a2 37 SetCurvature(fabs((kCLight*nanosecond/meter*q*B/tesla)/
69c1c8ff 38 (abs(p.Mag())/GeV*fCosDipAngle)/meter));
d0e92d9a 39#ifndef ST_NO_NAMESPACES
40 }
41#endif
42}
43
44AliFmPhysicalHelix::AliFmPhysicalHelix(double c, double d, double phase,
45 const AliFmThreeVector<double>& o, int h)
46 : AliFmHelix(c, d, phase, o, h) { /* nop */}
47
48
49AliFmThreeVector<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
81a888a2 59 double pt = GeV*fabs(kCLight*nanosecond/meter*B/tesla)/(fabs(fCurvature)*meter);
d0e92d9a 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
70AliFmThreeVector<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);
69c1c8ff 77 return (this->Momentum(B)).PseudoProduct(rx,ry,1.0);
d0e92d9a 78}
79
80int AliFmPhysicalHelix::Charge(double B) const
81{
82 // charge
83 return (B > 0 ? -fH : fH);
84}
85
86double 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);
a19edcc9 91 tDCA2dPosition.SetZ(0);
d0e92d9a 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);
a19edcc9 98 momVec.SetZ(0);
d0e92d9a 99 }
100 else {
101 momVec = this->MomentumAt(thePath,1./tesla); // Don't care about Bmag. Helicity is what matters.
a19edcc9 102 momVec.SetZ(0);
d0e92d9a 103 }
104
105 double cross = tDCAVec.x()*momVec.y() - tDCAVec.y()*momVec.x();
106 double theSign = (cross>=0) ? 1. : -1.;
69c1c8ff 107 return theSign*tDCAVec.Perp();
d0e92d9a 108}
109
110double 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
122double 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
130double 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}