MathTools/src/Helix.cc

// $Id: Helix.cc,v 1.2 2008/09/04 13:55:30 loizides Exp $

#include "MitCommon/MathTools/interface/Helix.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 !!!! STOP !!!!\n";
  //  Line             auxLine(PositionCm,MomentumGev.unit());  
  //  Z0             = auxLine.Z0();
  //  Phi0           = auxLine.Phi0();
  //  CotTheta       = auxLine.CotTheta();
  //  D0             = auxLine.D0();
  //  W              = 0.0;
  //  
  //  fCotTheta      = CotTheta;
  //  fCurvature     = W/2;
  //  fZ0            = Z0;
  //  fD0            = D0;
  //  fPhi0          = Phi0;
  //  fIsStale       = 1;
  //  fS             = -999.999;
  //  fAa            = -999.999; 
  //  fSs            = -999.999; 
  //  fCc            = -999.999; 
  //  fSinPhi0       = 1.0;
  //  fCosPhi0       = 1.0;
  //  fSinTheta      = 1.0;
  //  fCosTheta      = 1.0;
  //  fVParameters   = 0;
  //  fCenterIsValid = false;
  //  fMx            = 0.0;
  //  fMy            = 0.0;
  }
}

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;
}

//void Helix::Location(Trajectory::Location &loc, double s) const
//{
//  fCacheSinesAndCosines(s);
//  double cP0sT = fCosPhi0*fSinTheta, sP0sT = fSinPhi0*fSinTheta;
//  if (s && fCurvature) {
//    loc.SetLocation(s,
//                    (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),
//                    s*fCosTheta + fZ0,
//                    cP0sT*fCc-sP0sT*fSs,
//                    cP0sT*fSs+sP0sT*fCc,
//                    fCosTheta
//                    );
//  }
//  else {
//    loc.SetLocation(s,
//                    -fD0*fSinPhi0+s*cP0sT,
//                    fD0*fCosPhi0+s*sP0sT,
//                    fZ0+s*fCosTheta,
//                    cP0sT,sP0sT,fCosTheta
//                    );
//  }
//}

//Trajectory::Location* Helix::newIntersectionWith (const HepPlane3D& plane) const
//{
//  if (!SinTheta())                  // fastest way out of a screwy situation.
//    return 0;
//
//  double deltaS=1.0,s=0.0;
//  s = fabs(plane.d())/SinTheta();
//  HepVector3D normal = plane.normal();
//  Trajectory::Location *ploc = new Trajectory::Location ;
//  for (int iteration=0;iteration<100;iteration++) {
//    if (fabs(deltaS)>0.0001){
//      Location(*ploc,s);
//      deltaS = ((-(plane.distance(ploc->position())))
//              /(normal.dot(ploc->direction()))); 
//      s+=deltaS;
//    }
//    else {
//      if (s>0)
//      return ploc;
//    }
//  }
//  delete ploc;
//  return 0;
//}
Revision:	1.1
Committed:	Wed Sep 17 04:01:50 2008 UTC (16 years, 7 months ago) by loizides
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_008pre2, Mit_008pre1, Mit_006b, Mit_006a, Mit_006, Mit_005, Mit_004
Log Message:	Moved MitVertex contents to MitCommon. MitVertex therefore is obsolute and should not be touched anymore.
#	User	Rev	Content
1	loizides	1.1	// $Id: Helix.cc,v 1.2 2008/09/04 13:55:30 loizides Exp $
2
3			#include "MitCommon/MathTools/interface/Helix.h"
4
5			using namespace mithep;
6
7			Helix::Helix(const TVector3 &MomentumGev, const TVector3 &PositionCm,
8			double q, double BFieldTesla)
9			{
10			double CotTheta,W,Z0,D0;
11			Angle Phi0;
12
13			if (BFieldTesla != 0.0 && q != 0.0) {
14			double CurvatureConstant
15			= 0.0029979;
16			double Helicity = -1.0fabs(BFieldTesla)fabs(q)/(BFieldTesla*q);
17			double Radius = fabs(MomentumGev.Perp()/
18			(CurvatureConstantBFieldTeslaq));
19
20			if (Radius==0.0)
21			W = HUGE_VAL;
22			else
23			W = Helicity/Radius;
24			Angle phi1 = MomentumGev.Phi();
25			double x = PositionCm.x(),
26			y = PositionCm.y(),
27			z = PositionCm.z();
28			double sinPhi1 = sin(phi1),
29			cosPhi1 = cos(phi1);
30			double gamma = atan((xcosPhi1 + ysinPhi1)/(xsinPhi1-ycosPhi1 -1/W));
31			Phi0 = phi1+gamma;
32			D0 = ((1/W + ycosPhi1 - xsinPhi1) /cos(gamma) - 1/W);
33			CotTheta = MomentumGev.z()/MomentumGev.Perp();
34			Z0 = z + gamma*CotTheta/W;
35			}
36			// For the special case that we have no Helix, never happens for us .. right?
37			else {
38			std::cout << "MAJOR ERROR in HELIX !!!! STOP !!!!\n";
39			// Line auxLine(PositionCm,MomentumGev.unit());
40			// Z0 = auxLine.Z0();
41			// Phi0 = auxLine.Phi0();
42			// CotTheta = auxLine.CotTheta();
43			// D0 = auxLine.D0();
44			// W = 0.0;
45			//
46			// fCotTheta = CotTheta;
47			// fCurvature = W/2;
48			// fZ0 = Z0;
49			// fD0 = D0;
50			// fPhi0 = Phi0;
51			// fIsStale = 1;
52			// fS = -999.999;
53			// fAa = -999.999;
54			// fSs = -999.999;
55			// fCc = -999.999;
56			// fSinPhi0 = 1.0;
57			// fCosPhi0 = 1.0;
58			// fSinTheta = 1.0;
59			// fCosTheta = 1.0;
60			// fVParameters = 0;
61			// fCenterIsValid = false;
62			// fMx = 0.0;
63			// fMy = 0.0;
64			}
65			}
66
67			bool Helix::operator == (const Helix &right) const
68			{
69			return
70			fCotTheta == right.fCotTheta &&
71			fCurvature == right.fCurvature &&
72			fZ0 == right.fZ0 &&
73			fD0 == right.fD0 &&
74			fPhi0 == right.fPhi0;
75			}
76
77			bool Helix::operator != (const Helix & right) const
78			{
79			return !((*this)==right);
80			}
81
82			Helix::~Helix()
83			{
84			delete fVParameters;
85			}
86
87
88
89			TVector3 Helix::Position(double s) const
90			{
91			fCacheSinesAndCosines(s);
92			if (s == 0.0 \|\| fCurvature == 0.0) {
93			return TVector3(-fD0fSinPhi0+sfCosPhi0*fSinTheta,
94			fD0fCosPhi0+sfSinPhi0*fSinTheta,
95			fZ0+s*fCosTheta);
96			} else {
97			return TVector3((fCosPhi0fSs-fSinPhi0(2.0fCurvaturefD0+1.0-fCc))
98			/(2.0*fCurvature),
99			(fSinPhi0fSs+fCosPhi0(2.0fCurvaturefD0+1.0-fCc))
100			/(2.0*fCurvature),
101			fS*fCosTheta + fZ0);
102			}
103			}
104
105			TVector3 Helix::Direction(double s) const
106			{
107			fCacheSinesAndCosines(s);
108			if (s == 0.0) {
109			return TVector3(fCosPhi0fSinTheta,fSinPhi0fSinTheta,fCosTheta);
110			}
111			double xtan = fSinTheta(fCosPhi0fCc -fSinPhi0*fSs);
112			double ytan = fSinTheta(fCosPhi0fSs +fSinPhi0*fCc);
113			double ztan = fCosTheta;
114			return TVector3(xtan,ytan,ztan);
115			}
116
117			double Helix::PathLengthAtRhoEquals(double rho) const
118			{
119			return (SinTheta()?(L2DAtR(rho)/SinTheta()):0.0);
120			}
121
122			double Helix::InverseRadius() const
123			{
124			return fCurvature*2.0;
125			}
126
127			double Helix::Radius() const
128			{
129			return fabs(1.0/InverseRadius());
130			}
131
132			SignedAngle Helix::TurningAngle(double s) const
133			{
134			return s/Radius();
135			}
136
137			double Helix::Curvature() const
138			{
139			return fCurvature;
140			}
141
142			double Helix::Helicity() const
143			{
144			return fCurvature>0 ? 1.0 : -1.0 ;
145			}
146
147			double Helix::CotTheta() const
148			{
149			return fCotTheta;
150			}
151
152			Angle Helix::Phi0() const
153			{
154			return fPhi0;
155			}
156
157			double Helix::D0() const
158			{
159			return fD0;
160			}
161
162			double Helix::SignLz() const
163			{
164			return (fD0>0) ? -1.0 : 1.0;
165			}
166
167			double Helix::Z0() const
168			{
169			return fZ0;
170			}
171
172
173			TVector3 Helix::SecondDerivative(double s) const
174			{
175			double phi1 = fPhi0+s2.0fCurvature*fSinTheta;
176			double sp1 = sin(phi1);
177			double xsecond = -fSinThetasp12.0fCurvaturefSinTheta;
178			double ysecond = fSinThetasqrt(1.0-sp1sp1)2.0fCurvature*fSinTheta;
179			return TVector3(xsecond,ysecond,0.0);
180			}
181
182			double Helix::SinPhi0() const
183			{
184			fRefreshCache();
185			return fSinPhi0;
186			}
187			double Helix::CosPhi0() const
188			{
189			fRefreshCache();
190			return fCosPhi0;
191			}
192			double Helix::SinTheta() const
193			{
194			fRefreshCache();
195			return fSinTheta;
196			}
197			double Helix::CosTheta() const
198			{
199			fRefreshCache();
200			return fCosTheta;
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			double Helix::ZAtR(double rho) const
218			{
219			return fZ0 + CotTheta()*L2DAtR(rho);
220			}
221
222			double Helix::L2DAtR(double rho) const {
223			double L2D;
224
225			double c = Curvature();
226			double d = D0();
227
228			if (c!=0.0) {
229			double rad = (rhorho-dd)/(1.0+2.0cd);
230			double rprime;
231
232			if (rad<0.0)
233			rprime = 0.0;
234			else
235			rprime = sqrt( rad );
236
237			if (crprime > 1.0 \|\| crprime < -1.0)
238			L2D = c*rprime > 0. ? M_PI/c : -M_PI/c;
239			else
240			L2D = asin(c*rprime)/c;
241
242			}
243			else {
244			L2D = rho;
245			}
246			return L2D;
247			}
248
249			double Helix::CosAlphaAtR(double rho) const {
250			double c = Curvature();
251			double d = D0();
252			double a = c/(1.0+2.0cd);
253			double phi = PhiAtR(rho);
254
255			double x = rho*cos(phi);
256			double y = rho*sin(phi);
257
258			double dir_x0 = (1.0+2.0cd) * (1.0-2.0ay);
259			double dir_y0 = (1.0+2.0cd) * 2.0ax;
260
261			double phi0 = Phi0();
262			double dir_x = dir_x0cos(phi0) - dir_y0sin(phi0);
263			double dir_y = dir_x0sin(phi0) + dir_y0cos(phi0);
264
265			double cosalpha = (xdir_x + ydir_y)/rho;
266			return cosalpha;
267			}
268
269			//void Helix::Location(Trajectory::Location &loc, double s) const
270			//{
271			// fCacheSinesAndCosines(s);
272			// double cP0sT = fCosPhi0fSinTheta, sP0sT = fSinPhi0fSinTheta;
273			// if (s && fCurvature) {
274			// loc.SetLocation(s,
275			// (fCosPhi0*fSs-
276			// fSinPhi0(2.0fCurvature*fD0+1.0-fCc))
277			// /(2.0*fCurvature),
278			// (fSinPhi0*fSs+
279			// fCosPhi0(2.0fCurvature*fD0+1.0-fCc))
280			// /(2.0*fCurvature),
281			// s*fCosTheta + fZ0,
282			// cP0sTfCc-sP0sTfSs,
283			// cP0sTfSs+sP0sTfCc,
284			// fCosTheta
285			// );
286			// }
287			// else {
288			// loc.SetLocation(s,
289			// -fD0fSinPhi0+scP0sT,
290			// fD0fCosPhi0+ssP0sT,
291			// fZ0+s*fCosTheta,
292			// cP0sT,sP0sT,fCosTheta
293			// );
294			// }
295			//}
296
297			//Trajectory::Location* Helix::newIntersectionWith (const HepPlane3D& plane) const
298			//{
299			// if (!SinTheta()) // fastest way out of a screwy situation.
300			// return 0;
301			//
302			// double deltaS=1.0,s=0.0;
303			// s = fabs(plane.d())/SinTheta();
304			// HepVector3D normal = plane.normal();
305			// Trajectory::Location *ploc = new Trajectory::Location ;
306			// for (int iteration=0;iteration<100;iteration++) {
307			// if (fabs(deltaS)>0.0001){
308			// Location(*ploc,s);
309			// deltaS = ((-(plane.distance(ploc->position())))
310			// /(normal.dot(ploc->direction())));
311			// s+=deltaS;
312			// }
313			// else {
314			// if (s>0)
315			// return ploc;
316			// }
317			// }
318			// delete ploc;
319			// return 0;
320			//}