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
#	User	Rev	Content
1	loizides	1.3	// $Id: Helix.cc,v 1.2 2009/03/20 13:33:19 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	loizides	1.3	ClassImp(mithep::Helix)
7
8	loizides	1.1	using namespace mithep;
9
10	loizides	1.2	//--------------------------------------------------------------------------------------------------
11	loizides	1.1	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	loizides	1.2	std::cout << "MAJOR ERROR in HELIX!!!! Should not happen. STOP!!!!\n";
43			gSystem->Exit(123);
44	loizides	1.1	}
45			}
46
47	loizides	1.2	//--------------------------------------------------------------------------------------------------
48	loizides	1.1	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	loizides	1.2	//--------------------------------------------------------------------------------------------------
59	loizides	1.1	bool Helix::operator != (const Helix & right) const
60			{
61			return !((*this)==right);
62			}
63
64	loizides	1.2	//--------------------------------------------------------------------------------------------------
65	loizides	1.1	Helix::~Helix()
66			{
67			delete fVParameters;
68			}
69
70	loizides	1.2	//--------------------------------------------------------------------------------------------------
71	loizides	1.1	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	loizides	1.2	//--------------------------------------------------------------------------------------------------
88	loizides	1.1	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	loizides	1.2	//--------------------------------------------------------------------------------------------------
101	loizides	1.1	double Helix::PathLengthAtRhoEquals(double rho) const
102			{
103			return (SinTheta()?(L2DAtR(rho)/SinTheta()):0.0);
104			}
105
106	loizides	1.2	//--------------------------------------------------------------------------------------------------
107	loizides	1.1	double Helix::InverseRadius() const
108			{
109			return fCurvature*2.0;
110			}
111
112	loizides	1.2	//--------------------------------------------------------------------------------------------------
113	loizides	1.1	double Helix::Radius() const
114			{
115			return fabs(1.0/InverseRadius());
116			}
117
118	loizides	1.2	//--------------------------------------------------------------------------------------------------
119	loizides	1.1	SignedAngle Helix::TurningAngle(double s) const
120			{
121			return s/Radius();
122			}
123
124	loizides	1.2	//--------------------------------------------------------------------------------------------------
125	loizides	1.1	double Helix::Curvature() const
126			{
127			return fCurvature;
128			}
129
130	loizides	1.2	//--------------------------------------------------------------------------------------------------
131	loizides	1.1	double Helix::Helicity() const
132			{
133			return fCurvature>0 ? 1.0 : -1.0 ;
134			}
135
136	loizides	1.2	//--------------------------------------------------------------------------------------------------
137	loizides	1.1	double Helix::CotTheta() const
138			{
139			return fCotTheta;
140			}
141
142	loizides	1.2	//--------------------------------------------------------------------------------------------------
143	loizides	1.1	Angle Helix::Phi0() const
144			{
145			return fPhi0;
146			}
147
148	loizides	1.2	//--------------------------------------------------------------------------------------------------
149	loizides	1.1	double Helix::D0() const
150			{
151			return fD0;
152			}
153
154	loizides	1.2	//--------------------------------------------------------------------------------------------------
155	loizides	1.1	double Helix::SignLz() const
156			{
157			return (fD0>0) ? -1.0 : 1.0;
158			}
159
160	loizides	1.2	//--------------------------------------------------------------------------------------------------
161	loizides	1.1	double Helix::Z0() const
162			{
163			return fZ0;
164			}
165
166	loizides	1.2	//--------------------------------------------------------------------------------------------------
167	loizides	1.1	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	loizides	1.2	//--------------------------------------------------------------------------------------------------
177	loizides	1.1	double Helix::SinPhi0() const
178			{
179			fRefreshCache();
180			return fSinPhi0;
181			}
182	loizides	1.2
183			//--------------------------------------------------------------------------------------------------
184	loizides	1.1	double Helix::CosPhi0() const
185			{
186			fRefreshCache();
187			return fCosPhi0;
188			}
189	loizides	1.2
190			//--------------------------------------------------------------------------------------------------
191	loizides	1.1	double Helix::SinTheta() const
192			{
193			fRefreshCache();
194			return fSinTheta;
195			}
196	loizides	1.2
197			//--------------------------------------------------------------------------------------------------
198	loizides	1.1	double Helix::CosTheta() const
199			{
200			fRefreshCache();
201			return fCosTheta;
202			}
203
204	loizides	1.2	//--------------------------------------------------------------------------------------------------
205	loizides	1.1	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	loizides	1.2	//--------------------------------------------------------------------------------------------------
220	loizides	1.1	double Helix::ZAtR(double rho) const
221			{
222			return fZ0 + CotTheta()*L2DAtR(rho);
223			}
224
225	loizides	1.2	//--------------------------------------------------------------------------------------------------
226	loizides	1.1	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	loizides	1.2	//--------------------------------------------------------------------------------------------------
254	loizides	1.1	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			}