CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HelixArbitraryPlaneCrossing2Order.cc
Go to the documentation of this file.
2 
5 
6 #include <cmath>
7 #include <cfloat>
8 
10  const DirectionType& direction,
11  const float curvature,
12  const PropagationDirection propDir) :
13  theX0(point.x()),
14  theY0(point.y()),
15  theZ0(point.z()),
16  theRho(curvature),
17  thePropDir(propDir)
18 {
19  //
20  // Components of direction vector (with correct normalisation)
21  //
22  double px = direction.x();
23  double py = direction.y();
24  double pz = direction.z();
25  double pt2 = px*px+py*py;
26  double p2 = pt2+pz*pz;
27  double pI = 1./sqrt(p2);
28  double ptI = 1./sqrt(pt2);
29  theCosPhi0 = px*ptI;
30  theSinPhi0 = py*ptI;
31  theCosTheta = pz*pI;
32  theSinThetaI = p2*pI*ptI; // (1/(pt/p)) = p/pt = p*ptI and p = p2/p = p2*pI
33 }
34 
35 //
36 // Propagation status and path length to intersection
37 //
38 std::pair<bool,double>
40  //
41  // get local z-vector in global co-ordinates and
42  // distance to starting point
43  //
44  GlobalVector normalToPlane = plane.normalVector();
45  double nPx = normalToPlane.x();
46  double nPy = normalToPlane.y();
47  double nPz = normalToPlane.z();
48  double cP = plane.localZ(GlobalPoint(theX0,theY0,theZ0));
49  //
50  // coefficients of 2nd order equation to obtain intersection point
51  // in approximation (without curvature-related factors).
52  //
53  double ceq1 = theRho*(nPx*theSinPhi0-nPy*theCosPhi0);
54  double ceq2 = nPx*theCosPhi0 + nPy*theSinPhi0 + nPz*theCosTheta*theSinThetaI;
55  double ceq3 = cP;
56  //
57  // Check for degeneration to linear equation (zero
58  // curvature, forward plane or direction perp. to plane)
59  //
60  double dS1,dS2;
61  if ( fabs(ceq1)>FLT_MIN ) {
62  double deq1 = ceq2*ceq2;
63  double deq2 = ceq1*ceq3;
64  if ( fabs(deq1)<FLT_MIN || fabs(deq2/deq1)>1.e-6 ) {
65  //
66  // Standard solution for quadratic equations
67  //
68  double deq = deq1+2*deq2;
69  if ( deq<0. ) return std::pair<bool,double>(false,0);
70  double ceq = -0.5*(ceq2+(ceq2>0?1:-1)*sqrt(deq));
71  dS1 = -2*(ceq/ceq1)*theSinThetaI;
72  dS2 = (ceq3/ceq)*theSinThetaI;
73  }
74  else {
75  //
76  // Solution by expansion of sqrt(1+deq)
77  //
78  double ceq = (ceq2/ceq1)*theSinThetaI;
79  double deq = deq2/deq1;
80  deq *= (1-0.5*deq);
81  dS1 = -ceq*deq;
82  dS2 = ceq*(2+deq);
83  }
84  }
85  else {
86  //
87  // Special case: linear equation
88  //
89  dS1 = dS2 = -(ceq3/ceq2)*theSinThetaI;
90  }
91  //
92  // Choose and return solution
93  //
94  return solutionByDirection(dS1,dS2);
95 }
96 
97 //
98 // Position after a step of path length s (2nd order)
99 //
102  // use double precision result
104  return PositionType(pos.x(),pos.y(),pos.z());
105 }
106 //
107 // Position after a step of path length s (2nd order) (in double precision)
108 //
111  // based on path length in the transverse plane
112  double st = s/theSinThetaI;
113  return PositionTypeDouble(theX0+(theCosPhi0-(st*0.5*theRho)*theSinPhi0)*st,
114  theY0+(theSinPhi0+(st*0.5*theRho)*theCosPhi0)*st,
116 }
117 //
118 // Direction after a step of path length 2 (2nd order) (in double precision)
119 //
122  // use double precision result
124  return DirectionType(dir.x(),dir.y(),dir.z());
125 }
126 //
127 // Direction after a step of path length 2 (2nd order)
128 //
131  // based on delta phi
132  double dph = s*theRho/theSinThetaI;
134  theSinPhi0+(theCosPhi0-0.5*dph*theSinPhi0)*dph,
136 }
137 //
138 // Choice of solution according to propagation direction
139 //
140 std::pair<bool,double>
142  const double dS2) const {
143  bool valid = false;
144  double path = 0;
145  if ( thePropDir == anyDirection ) {
146  valid = true;
147  path = smallestPathLength(dS1,dS2);
148  }
149  else {
150  // use same logic for both directions (invert if necessary)
151  double propSign = thePropDir==alongMomentum ? 1 : -1;
152  double s1(propSign*dS1);
153  double s2(propSign*dS2);
154  // sort
155  if ( s1 > s2 ) {
156  double aux = s1;
157  s1 = s2;
158  s2 = aux;
159  }
160  // choose solution (if any with positive sign)
161  if ( s1<0 && s2>=0 ) {
162  // First solution in backward direction: choose second one.
163  valid = true;
164  path = propSign*s2;
165  }
166  else if ( s1>=0 ) {
167  // First solution in forward direction: choose it (s2 is further away!).
168  valid = true;
169  path = propSign*s1;
170  }
171  }
172  return std::pair<bool,double>(valid,path);
173 }
HelixArbitraryPlaneCrossing2Order(const PositionType &point, const DirectionType &direction, const float curvature, const PropagationDirection propDir=alongMomentum)
T y() const
Cartesian y coordinate.
T x() const
Cartesian x coordinate.
DirectionTypeDouble directionInDouble(double s) const
GlobalVector normalVector() const
Definition: Plane.h:47
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
virtual std::pair< bool, double > pathLength(const Plane &)
T y() const
Definition: PV3DBase.h:57
float localZ(const GlobalPoint &gp) const
Fast access to distance from plane for a point.
Definition: Plane.h:52
PropagationDirection
double smallestPathLength(const double firstPathLength, const double secondPathLength) const
Definition: Plane.h:17
double double double z
tuple s2
Definition: indexGen.py:106
int path() const
Definition: HLTadd.h:3
T curvature(T InversePt, const edm::EventSetup &iSetup)
T z() const
Cartesian z coordinate.
T sqrt(T t)
Definition: SSEVec.h:28
T z() const
Definition: PV3DBase.h:58
PositionTypeDouble positionInDouble(double s) const
std::pair< bool, double > solutionByDirection(const double dS1, const double dS2) const
Basic3DVector< float > PositionType
the helix is passed to the constructor and does not appear in the interface
Basic3DVector< float > DirectionType
double p2[4]
Definition: TauolaWrapper.h:90
virtual PositionType position(double s) const
virtual DirectionType direction(double s) const
string s
Definition: asciidump.py:422
dbl *** dir
Definition: mlp_gen.cc:35
Definition: DDAxes.h:10
T x() const
Definition: PV3DBase.h:56
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
Definition: invegas.h:5