MathTools/src/Helix.cc

// $Id: Helix.cc,v 1.2 2009/03/20 13:33:19 loizides Exp $

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

ClassImp(mithep::Helix)

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