MathTools/src/Helix.cc

// $Id: Helix.cc,v 1.1 2008/09/17 04:01:50 loizides Exp $

#include "MitCommon/MathTools/interface/Helix.h"
#include <TSystem.h>

using namespace mithep;

//--------------------------------------------------------------------------------------------------
Helix::Helix(const TVector3 &MomentumGev, const TVector3 &PositionCm,
             double q, double BFieldTesla)
{
  double CotTheta,W,Z0,D0;
  Angle  Phi0;

  if (BFieldTesla != 0.0 && q != 0.0) {
    double CurvatureConstant
                    = 0.0029979;
    double Helicity = -1.0*fabs(BFieldTesla)*fabs(q)/(BFieldTesla*q);
    double Radius   = fabs(MomentumGev.Perp()/
                           (CurvatureConstant*BFieldTesla*q));
    
    if (Radius==0.0)
      W = HUGE_VAL;
    else
      W = Helicity/Radius;
    Angle  phi1     = MomentumGev.Phi();
    double x        = PositionCm.x(),
           y        = PositionCm.y(),
           z        = PositionCm.z(); 
    double sinPhi1  = sin(phi1),
           cosPhi1  = cos(phi1);
    double gamma    = atan((x*cosPhi1 + y*sinPhi1)/(x*sinPhi1-y*cosPhi1 -1/W));
    Phi0            = phi1+gamma;
    D0              = ((1/W + y*cosPhi1 - x*sinPhi1) /cos(gamma) - 1/W);
    CotTheta        = MomentumGev.z()/MomentumGev.Perp();
    Z0              = z + gamma*CotTheta/W;
  }
  // For the special case that we have no Helix, never happens for us .. right?
  else {
    std::cout << "MAJOR ERROR in HELIX!!!! Should not happen. STOP!!!!\n";
    gSystem->Exit(123);
  }
}

//--------------------------------------------------------------------------------------------------
bool Helix::operator == (const Helix &right) const
{
  return
    fCotTheta  == right.fCotTheta  &&
    fCurvature == right.fCurvature &&
    fZ0        == right.fZ0        &&
    fD0        == right.fD0        &&
    fPhi0      == right.fPhi0;
}

//--------------------------------------------------------------------------------------------------
bool Helix::operator != (const Helix & right) const
{
  return !((*this)==right);
}

//--------------------------------------------------------------------------------------------------
Helix::~Helix()
{
  delete fVParameters;
}

//--------------------------------------------------------------------------------------------------
TVector3 Helix::Position(double s) const
{
  fCacheSinesAndCosines(s);
  if (s == 0.0 || fCurvature == 0.0) {
    return TVector3(-fD0*fSinPhi0+s*fCosPhi0*fSinTheta,
                     fD0*fCosPhi0+s*fSinPhi0*fSinTheta,
                     fZ0+s*fCosTheta);
  } else {
    return TVector3((fCosPhi0*fSs-fSinPhi0*(2.0*fCurvature*fD0+1.0-fCc))
                    /(2.0*fCurvature),
                    (fSinPhi0*fSs+fCosPhi0*(2.0*fCurvature*fD0+1.0-fCc))
                    /(2.0*fCurvature),   
                    fS*fCosTheta + fZ0);
  }
}

//--------------------------------------------------------------------------------------------------
TVector3 Helix::Direction(double s) const
{
  fCacheSinesAndCosines(s);
  if (s == 0.0) {
    return TVector3(fCosPhi0*fSinTheta,fSinPhi0*fSinTheta,fCosTheta);
  }
  double xtan = fSinTheta*(fCosPhi0*fCc -fSinPhi0*fSs); 
  double ytan = fSinTheta*(fCosPhi0*fSs +fSinPhi0*fCc); 
  double ztan = fCosTheta;
  return TVector3(xtan,ytan,ztan);
}

//--------------------------------------------------------------------------------------------------
double Helix::PathLengthAtRhoEquals(double rho) const
{
  return (SinTheta()?(L2DAtR(rho)/SinTheta()):0.0);
}

//--------------------------------------------------------------------------------------------------
double Helix::InverseRadius() const
{
  return fCurvature*2.0;
}

//--------------------------------------------------------------------------------------------------
double Helix::Radius() const
{
  return fabs(1.0/InverseRadius());
}

//--------------------------------------------------------------------------------------------------
SignedAngle Helix::TurningAngle(double s) const
{
  return s/Radius();
}

//--------------------------------------------------------------------------------------------------
double Helix::Curvature() const
{
  return fCurvature;
}

//--------------------------------------------------------------------------------------------------
double Helix::Helicity() const
{
  return fCurvature>0 ? 1.0 : -1.0 ;
}

//--------------------------------------------------------------------------------------------------
double Helix::CotTheta() const
{
  return fCotTheta;
}

//--------------------------------------------------------------------------------------------------
Angle Helix::Phi0() const
{
  return fPhi0;
}

//--------------------------------------------------------------------------------------------------
double Helix::D0() const
{
  return fD0;
}

//--------------------------------------------------------------------------------------------------
double Helix::SignLz() const
{
  return (fD0>0) ? -1.0 : 1.0;
}

//--------------------------------------------------------------------------------------------------
double Helix::Z0() const
{
  return fZ0;
}

//--------------------------------------------------------------------------------------------------
TVector3 Helix::SecondDerivative(double s) const
{
  double phi1    = fPhi0+s*2.0*fCurvature*fSinTheta;
  double sp1     = sin(phi1);
  double xsecond = -fSinTheta*sp1*2.0*fCurvature*fSinTheta;
  double ysecond = fSinTheta*sqrt(1.0-sp1*sp1)*2.0*fCurvature*fSinTheta;
  return TVector3(xsecond,ysecond,0.0);
}

//--------------------------------------------------------------------------------------------------
double Helix::SinPhi0() const
{
  fRefreshCache();
  return fSinPhi0;
}

//--------------------------------------------------------------------------------------------------
double Helix::CosPhi0() const
{
  fRefreshCache();
  return fCosPhi0;
}

//--------------------------------------------------------------------------------------------------
double Helix::SinTheta() const
{
  fRefreshCache();
  return fSinTheta;
}

//--------------------------------------------------------------------------------------------------
double Helix::CosTheta() const
{
  fRefreshCache();
  return fCosTheta;
}

//--------------------------------------------------------------------------------------------------
Angle Helix::PhiAtR(double rho) const
{
  double c = Curvature();
  double d = D0();
  double a = c/(1.0+2.0*c*d);
  double b = d - a*d*d;
  double arcsin = a*rho+b/rho;
  if (arcsin>1.0 || arcsin<-1.0) {
    return (arcsin > 0.) ? M_PI : -M_PI;
  }
  Angle phi = Phi0() + asin(arcsin);
  return phi;
}

//--------------------------------------------------------------------------------------------------
double Helix::ZAtR(double rho) const
{
  return fZ0 + CotTheta()*L2DAtR(rho);
}

//--------------------------------------------------------------------------------------------------
double Helix::L2DAtR(double rho) const {
  double L2D;

  double c = Curvature();
  double d = D0();

  if (c!=0.0) {
    double rad = (rho*rho-d*d)/(1.0+2.0*c*d);
    double rprime;

    if (rad<0.0)
      rprime = 0.0;
    else
      rprime = sqrt( rad );

    if (c*rprime > 1.0 || c*rprime < -1.0)
      L2D = c*rprime > 0. ? M_PI/c : -M_PI/c;
    else
      L2D = asin(c*rprime)/c;

  }
  else {
    L2D = rho;
  }
  return L2D;
}

//--------------------------------------------------------------------------------------------------
double Helix::CosAlphaAtR(double rho) const {
  double c   = Curvature();
  double d   = D0();
  double a   = c/(1.0+2.0*c*d);
  double phi = PhiAtR(rho);

  double x   = rho*cos(phi);
  double y   = rho*sin(phi);

  double dir_x0 = (1.0+2.0*c*d) * (1.0-2.0*a*y);
  double dir_y0 = (1.0+2.0*c*d) * 2.0*a*x;

  double phi0 = Phi0();
  double dir_x = dir_x0*cos(phi0) - dir_y0*sin(phi0);
  double dir_y = dir_x0*sin(phi0) + dir_y0*cos(phi0);

  double cosalpha = (x*dir_x + y*dir_y)/rho;
  return cosalpha;
}
Revision:	1.2
Committed:	Fri Mar 20 13:33:19 2009 UTC (16 years, 1 month ago) by loizides
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_009c, Mit_009b, Mit_009a, Mit_009
Changes since 1.1:	+33 -83 lines
Log Message:	Cleanup
#	User	Rev	Content
1	loizides	1.2	// $Id: Helix.cc,v 1.1 2008/09/17 04:01:50 loizides Exp $
2	loizides	1.1
3			#include "MitCommon/MathTools/interface/Helix.h"
4	loizides	1.2	#include <TSystem.h>
5	loizides	1.1
6			using namespace mithep;
7
8	loizides	1.2	//--------------------------------------------------------------------------------------------------
9	loizides	1.1	Helix::Helix(const TVector3 &MomentumGev, const TVector3 &PositionCm,
10			double q, double BFieldTesla)
11			{
12			double CotTheta,W,Z0,D0;
13			Angle Phi0;
14
15			if (BFieldTesla != 0.0 && q != 0.0) {
16			double CurvatureConstant
17			= 0.0029979;
18			double Helicity = -1.0fabs(BFieldTesla)fabs(q)/(BFieldTesla*q);
19			double Radius = fabs(MomentumGev.Perp()/
20			(CurvatureConstantBFieldTeslaq));
21
22			if (Radius==0.0)
23			W = HUGE_VAL;
24			else
25			W = Helicity/Radius;
26			Angle phi1 = MomentumGev.Phi();
27			double x = PositionCm.x(),
28			y = PositionCm.y(),
29			z = PositionCm.z();
30			double sinPhi1 = sin(phi1),
31			cosPhi1 = cos(phi1);
32			double gamma = atan((xcosPhi1 + ysinPhi1)/(xsinPhi1-ycosPhi1 -1/W));
33			Phi0 = phi1+gamma;
34			D0 = ((1/W + ycosPhi1 - xsinPhi1) /cos(gamma) - 1/W);
35			CotTheta = MomentumGev.z()/MomentumGev.Perp();
36			Z0 = z + gamma*CotTheta/W;
37			}
38			// For the special case that we have no Helix, never happens for us .. right?
39			else {
40	loizides	1.2	std::cout << "MAJOR ERROR in HELIX!!!! Should not happen. STOP!!!!\n";
41			gSystem->Exit(123);
42	loizides	1.1	}
43			}
44
45	loizides	1.2	//--------------------------------------------------------------------------------------------------
46	loizides	1.1	bool Helix::operator == (const Helix &right) const
47			{
48			return
49			fCotTheta == right.fCotTheta &&
50			fCurvature == right.fCurvature &&
51			fZ0 == right.fZ0 &&
52			fD0 == right.fD0 &&
53			fPhi0 == right.fPhi0;
54			}
55
56	loizides	1.2	//--------------------------------------------------------------------------------------------------
57	loizides	1.1	bool Helix::operator != (const Helix & right) const
58			{
59			return !((*this)==right);
60			}
61
62	loizides	1.2	//--------------------------------------------------------------------------------------------------
63	loizides	1.1	Helix::~Helix()
64			{
65			delete fVParameters;
66			}
67
68	loizides	1.2	//--------------------------------------------------------------------------------------------------
69	loizides	1.1	TVector3 Helix::Position(double s) const
70			{
71			fCacheSinesAndCosines(s);
72			if (s == 0.0 \|\| fCurvature == 0.0) {
73			return TVector3(-fD0fSinPhi0+sfCosPhi0*fSinTheta,
74			fD0fCosPhi0+sfSinPhi0*fSinTheta,
75			fZ0+s*fCosTheta);
76			} else {
77			return TVector3((fCosPhi0fSs-fSinPhi0(2.0fCurvaturefD0+1.0-fCc))
78			/(2.0*fCurvature),
79			(fSinPhi0fSs+fCosPhi0(2.0fCurvaturefD0+1.0-fCc))
80			/(2.0*fCurvature),
81			fS*fCosTheta + fZ0);
82			}
83			}
84
85	loizides	1.2	//--------------------------------------------------------------------------------------------------
86	loizides	1.1	TVector3 Helix::Direction(double s) const
87			{
88			fCacheSinesAndCosines(s);
89			if (s == 0.0) {
90			return TVector3(fCosPhi0fSinTheta,fSinPhi0fSinTheta,fCosTheta);
91			}
92			double xtan = fSinTheta(fCosPhi0fCc -fSinPhi0*fSs);
93			double ytan = fSinTheta(fCosPhi0fSs +fSinPhi0*fCc);
94			double ztan = fCosTheta;
95			return TVector3(xtan,ytan,ztan);
96			}
97
98	loizides	1.2	//--------------------------------------------------------------------------------------------------
99	loizides	1.1	double Helix::PathLengthAtRhoEquals(double rho) const
100			{
101			return (SinTheta()?(L2DAtR(rho)/SinTheta()):0.0);
102			}
103
104	loizides	1.2	//--------------------------------------------------------------------------------------------------
105	loizides	1.1	double Helix::InverseRadius() const
106			{
107			return fCurvature*2.0;
108			}
109
110	loizides	1.2	//--------------------------------------------------------------------------------------------------
111	loizides	1.1	double Helix::Radius() const
112			{
113			return fabs(1.0/InverseRadius());
114			}
115
116	loizides	1.2	//--------------------------------------------------------------------------------------------------
117	loizides	1.1	SignedAngle Helix::TurningAngle(double s) const
118			{
119			return s/Radius();
120			}
121
122	loizides	1.2	//--------------------------------------------------------------------------------------------------
123	loizides	1.1	double Helix::Curvature() const
124			{
125			return fCurvature;
126			}
127
128	loizides	1.2	//--------------------------------------------------------------------------------------------------
129	loizides	1.1	double Helix::Helicity() const
130			{
131			return fCurvature>0 ? 1.0 : -1.0 ;
132			}
133
134	loizides	1.2	//--------------------------------------------------------------------------------------------------
135	loizides	1.1	double Helix::CotTheta() const
136			{
137			return fCotTheta;
138			}
139
140	loizides	1.2	//--------------------------------------------------------------------------------------------------
141	loizides	1.1	Angle Helix::Phi0() const
142			{
143			return fPhi0;
144			}
145
146	loizides	1.2	//--------------------------------------------------------------------------------------------------
147	loizides	1.1	double Helix::D0() const
148			{
149			return fD0;
150			}
151
152	loizides	1.2	//--------------------------------------------------------------------------------------------------
153	loizides	1.1	double Helix::SignLz() const
154			{
155			return (fD0>0) ? -1.0 : 1.0;
156			}
157
158	loizides	1.2	//--------------------------------------------------------------------------------------------------
159	loizides	1.1	double Helix::Z0() const
160			{
161			return fZ0;
162			}
163
164	loizides	1.2	//--------------------------------------------------------------------------------------------------
165	loizides	1.1	TVector3 Helix::SecondDerivative(double s) const
166			{
167			double phi1 = fPhi0+s2.0fCurvature*fSinTheta;
168			double sp1 = sin(phi1);
169			double xsecond = -fSinThetasp12.0fCurvaturefSinTheta;
170			double ysecond = fSinThetasqrt(1.0-sp1sp1)2.0fCurvature*fSinTheta;
171			return TVector3(xsecond,ysecond,0.0);
172			}
173
174	loizides	1.2	//--------------------------------------------------------------------------------------------------
175	loizides	1.1	double Helix::SinPhi0() const
176			{
177			fRefreshCache();
178			return fSinPhi0;
179			}
180	loizides	1.2
181			//--------------------------------------------------------------------------------------------------
182	loizides	1.1	double Helix::CosPhi0() const
183			{
184			fRefreshCache();
185			return fCosPhi0;
186			}
187	loizides	1.2
188			//--------------------------------------------------------------------------------------------------
189	loizides	1.1	double Helix::SinTheta() const
190			{
191			fRefreshCache();
192			return fSinTheta;
193			}
194	loizides	1.2
195			//--------------------------------------------------------------------------------------------------
196	loizides	1.1	double Helix::CosTheta() const
197			{
198			fRefreshCache();
199			return fCosTheta;
200			}
201
202	loizides	1.2	//--------------------------------------------------------------------------------------------------
203	loizides	1.1	Angle Helix::PhiAtR(double rho) const
204			{
205			double c = Curvature();
206			double d = D0();
207			double a = c/(1.0+2.0cd);
208			double b = d - add;
209			double arcsin = a*rho+b/rho;
210			if (arcsin>1.0 \|\| arcsin<-1.0) {
211			return (arcsin > 0.) ? M_PI : -M_PI;
212			}
213			Angle phi = Phi0() + asin(arcsin);
214			return phi;
215			}
216
217	loizides	1.2	//--------------------------------------------------------------------------------------------------
218	loizides	1.1	double Helix::ZAtR(double rho) const
219			{
220			return fZ0 + CotTheta()*L2DAtR(rho);
221			}
222
223	loizides	1.2	//--------------------------------------------------------------------------------------------------
224	loizides	1.1	double Helix::L2DAtR(double rho) const {
225			double L2D;
226
227			double c = Curvature();
228			double d = D0();
229
230			if (c!=0.0) {
231			double rad = (rhorho-dd)/(1.0+2.0cd);
232			double rprime;
233
234			if (rad<0.0)
235			rprime = 0.0;
236			else
237			rprime = sqrt( rad );
238
239			if (crprime > 1.0 \|\| crprime < -1.0)
240			L2D = c*rprime > 0. ? M_PI/c : -M_PI/c;
241			else
242			L2D = asin(c*rprime)/c;
243
244			}
245			else {
246			L2D = rho;
247			}
248			return L2D;
249			}
250
251	loizides	1.2	//--------------------------------------------------------------------------------------------------
252	loizides	1.1	double Helix::CosAlphaAtR(double rho) const {
253			double c = Curvature();
254			double d = D0();
255			double a = c/(1.0+2.0cd);
256			double phi = PhiAtR(rho);
257
258			double x = rho*cos(phi);
259			double y = rho*sin(phi);
260
261			double dir_x0 = (1.0+2.0cd) * (1.0-2.0ay);
262			double dir_y0 = (1.0+2.0cd) * 2.0ax;
263
264			double phi0 = Phi0();
265			double dir_x = dir_x0cos(phi0) - dir_y0sin(phi0);
266			double dir_y = dir_x0sin(phi0) + dir_y0cos(phi0);
267
268			double cosalpha = (xdir_x + ydir_y)/rho;
269			return cosalpha;
270			}