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
#	Content
1	// $Id: Helix.cc,v 1.1 2008/09/17 04:01:50 loizides Exp $
2
3	#include "MitCommon/MathTools/interface/Helix.h"
4	#include <TSystem.h>
5
6	using namespace mithep;
7
8	//--------------------------------------------------------------------------------------------------
9	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	std::cout << "MAJOR ERROR in HELIX!!!! Should not happen. STOP!!!!\n";
41	gSystem->Exit(123);
42	}
43	}
44
45	//--------------------------------------------------------------------------------------------------
46	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	//--------------------------------------------------------------------------------------------------
57	bool Helix::operator != (const Helix & right) const
58	{
59	return !((*this)==right);
60	}
61
62	//--------------------------------------------------------------------------------------------------
63	Helix::~Helix()
64	{
65	delete fVParameters;
66	}
67
68	//--------------------------------------------------------------------------------------------------
69	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	//--------------------------------------------------------------------------------------------------
86	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	//--------------------------------------------------------------------------------------------------
99	double Helix::PathLengthAtRhoEquals(double rho) const
100	{
101	return (SinTheta()?(L2DAtR(rho)/SinTheta()):0.0);
102	}
103
104	//--------------------------------------------------------------------------------------------------
105	double Helix::InverseRadius() const
106	{
107	return fCurvature*2.0;
108	}
109
110	//--------------------------------------------------------------------------------------------------
111	double Helix::Radius() const
112	{
113	return fabs(1.0/InverseRadius());
114	}
115
116	//--------------------------------------------------------------------------------------------------
117	SignedAngle Helix::TurningAngle(double s) const
118	{
119	return s/Radius();
120	}
121
122	//--------------------------------------------------------------------------------------------------
123	double Helix::Curvature() const
124	{
125	return fCurvature;
126	}
127
128	//--------------------------------------------------------------------------------------------------
129	double Helix::Helicity() const
130	{
131	return fCurvature>0 ? 1.0 : -1.0 ;
132	}
133
134	//--------------------------------------------------------------------------------------------------
135	double Helix::CotTheta() const
136	{
137	return fCotTheta;
138	}
139
140	//--------------------------------------------------------------------------------------------------
141	Angle Helix::Phi0() const
142	{
143	return fPhi0;
144	}
145
146	//--------------------------------------------------------------------------------------------------
147	double Helix::D0() const
148	{
149	return fD0;
150	}
151
152	//--------------------------------------------------------------------------------------------------
153	double Helix::SignLz() const
154	{
155	return (fD0>0) ? -1.0 : 1.0;
156	}
157
158	//--------------------------------------------------------------------------------------------------
159	double Helix::Z0() const
160	{
161	return fZ0;
162	}
163
164	//--------------------------------------------------------------------------------------------------
165	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	//--------------------------------------------------------------------------------------------------
175	double Helix::SinPhi0() const
176	{
177	fRefreshCache();
178	return fSinPhi0;
179	}
180
181	//--------------------------------------------------------------------------------------------------
182	double Helix::CosPhi0() const
183	{
184	fRefreshCache();
185	return fCosPhi0;
186	}
187
188	//--------------------------------------------------------------------------------------------------
189	double Helix::SinTheta() const
190	{
191	fRefreshCache();
192	return fSinTheta;
193	}
194
195	//--------------------------------------------------------------------------------------------------
196	double Helix::CosTheta() const
197	{
198	fRefreshCache();
199	return fCosTheta;
200	}
201
202	//--------------------------------------------------------------------------------------------------
203	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	//--------------------------------------------------------------------------------------------------
218	double Helix::ZAtR(double rho) const
219	{
220	return fZ0 + CotTheta()*L2DAtR(rho);
221	}
222
223	//--------------------------------------------------------------------------------------------------
224	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	//--------------------------------------------------------------------------------------------------
252	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	}