CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HelixArbitraryPlaneCrossing.cc
Go to the documentation of this file.
4 
5 #include <cmath>
6 #include <iostream>
8 
10  const DirectionType& direction,
11  const float curvature,
12  const PropagationDirection propDir) :
13  theQuadraticCrossingFromStart(point,direction,curvature,propDir),
14  theX0(point.x()),
15  theY0(point.y()),
16  theZ0(point.z()),
17  theRho(curvature),
18  thePropDir(propDir),
19  theCachedS(0),
20  theCachedDPhi(0.),
21  theCachedSDPhi(0.),
22  theCachedCDPhi(1.)
23 {
24  //
25  // Components of direction vector (with correct normalisation)
26  //
27  double px = direction.x();
28  double py = direction.y();
29  double pz = direction.z();
30  double pt2 = px*px+py*py;
31  double p2 = pt2+pz*pz;
32  double pI = 1./sqrt(p2);
33  double ptI = 1./sqrt(pt2);
34  theCosPhi0 = px*ptI;
35  theSinPhi0 = py*ptI;
36  theCosTheta = pz*pI;
37  theSinTheta = pt2*ptI*pI;
38 }
39 //
40 // Propagation status and path length to intersection
41 //
42 std::pair<bool,double>
44  //
45  // Constants used for control of convergence
46  //
47  const int maxIterations(100);
48  //
49  // maximum distance to plane (taking into account numerical precision)
50  //
51  float maxNumDz = theNumericalPrecision*plane.position().mag();
52  float safeMaxDist = (theMaxDistToPlane>maxNumDz?theMaxDistToPlane:maxNumDz);
53  //
54  // Prepare internal value of the propagation direction and position / direction vectors for iteration
55  //
59  //
60  // Prepare iterations: count and total pathlength
61  //
62  unsigned int iteration(maxIterations+1);
63  double dSTotal(0.);
64  //
65  bool first = true;
66  while ( notAtSurface(plane,xnew,safeMaxDist) ) {
67  //
68  // return empty solution vector if no convergence after maxIterations iterations
69  //
70  if unlikely( --iteration == 0 ) {
71  edm::LogInfo("HelixArbitraryPlaneCrossing") << "pathLength : no convergence";
72  return std::pair<bool,double>(false,0);
73  }
74 
75  //
76  // Use existing 2nd order object at first pass, create temporary object
77  // for subsequent passes.
78  //
79  std::pair<bool,double> deltaS2;
80  if unlikely( first ) {
81  first = false;
83  }
84  else {
85  HelixArbitraryPlaneCrossing2Order quadraticCrossing(xnew.x(),xnew.y(),xnew.z(),
86  pnew.x(),pnew.y(),
88  theRho,
89  anyDirection);
90 
91  deltaS2 = quadraticCrossing.pathLength(plane);
92  }
93 
94  if unlikely( !deltaS2.first ) return deltaS2;
95  //
96  // Calculate and sort total pathlength (max. 2 solutions)
97  //
98  dSTotal += deltaS2.second;
100  if ( propDir == anyDirection ) {
101  propDir = newDir;
102  }
103  else {
104  if unlikely( newDir!=propDir ) return std::pair<bool,double>(false,0);
105  }
106  //
107  // Step forward by dSTotal.
108  //
109  xnew = positionInDouble(dSTotal);
110  pnew = directionInDouble(dSTotal);
111  }
112  //
113  // Return result
114  //
115  return std::pair<bool,double>(true,dSTotal);
116 }
117 //
118 // Position on helix after a step of path length s.
119 //
122  // use result in double precision
124  return PositionType(pos.x(),pos.y(),pos.z());
125 }
126 //
127 // Position on helix after a step of path length s in double precision.
128 //
131  //
132  // Calculate delta phi (if not already available)
133  //
134  if unlikely( s!=theCachedS ) {
135  theCachedS = s;
139  }
140  //
141  // Calculate with appropriate formulation of full helix formula or with
142  // 2nd order approximation.
143  //
144 // if ( fabs(theCachedDPhi)>1.e-1 ) {
145  if ( std::abs(theCachedDPhi)>1.e-4 ) {
146  // "standard" helix formula
147  double o = 1./theRho;
149  theY0+( theCosPhi0*(1.-theCachedCDPhi)+theSinPhi0*theCachedSDPhi)*o,
151  }
152 // else if ( fabs(theCachedDPhi)>theNumericalPrecision ) {
153 // // full helix formula, but avoiding (1-cos(deltaPhi)) for small angles
154 // return PositionTypeDouble(theX0+(-theSinPhi0*theCachedSDPhi*theCachedSDPhi/(1.+theCachedCDPhi)+
155 // theCosPhi0*theCachedSDPhi)/theRho,
156 // theY0+(theCosPhi0*theCachedSDPhi*theCachedSDPhi/(1.+theCachedCDPhi)+
157 // theSinPhi0*theCachedSDPhi)/theRho,
158 // theZ0+theCachedS*theCosTheta);
159 // }
160  else {
161  // Use 2nd order.
163  }
164 }
165 //
166 // Direction vector on helix after a step of path length s.
167 //
170  // use result in double precision
172  return DirectionType(dir.x(),dir.y(),dir.z());
173 }
174 //
175 // Direction vector on helix after a step of path length s in double precision.
176 //
179  //
180  // Calculate delta phi (if not already available)
181  //
182  if unlikely( s!=theCachedS ) {
183  theCachedS = s;
187  }
188 
189  if ( std::abs(theCachedDPhi)>1.e-4 ) {
190  // full helix formula
192  theSinPhi0*theCachedCDPhi+theCosPhi0*theCachedSDPhi,
194  }
195  else {
196  // 2nd order
198  }
199 }
200 // Iteration control: continue if distance to plane > theMaxDistToPlane. Includes
201 // protection for numerical precision (Surfaces work with single precision).
203  const PositionTypeDouble& point,
204  const float maxDist) const {
205  float dz = plane.localZ(Surface::GlobalPoint(point.x(),point.y(),point.z()));
206  return std::abs(dz)>maxDist;
207 }
208 
DirectionTypeDouble directionInDouble(double s) const
HelixArbitraryPlaneCrossing2Order theQuadraticCrossingFromStart
T y() const
Cartesian y coordinate.
T x() const
Cartesian x coordinate.
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
DirectionTypeDouble directionInDouble(double s) const
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
virtual std::pair< bool, double > pathLength(const Plane &)
Basic3DVector< double > PositionTypeDouble
float localZ(const GlobalPoint &gp) const
Fast access to distance from plane for a point.
Definition: Plane.h:52
const PropagationDirection thePropDir
#define abs(x)
Definition: mlp_lapack.h:159
virtual DirectionType direction(double s) const
PropagationDirection
virtual PositionType position(double s) const
Definition: Plane.h:17
double double double z
#define unlikely(x)
Definition: Likely.h:21
T curvature(T InversePt, const edm::EventSetup &iSetup)
tuple iteration
Definition: align_cfg.py:5
T z() const
Cartesian z coordinate.
bool notAtSurface(const Plane &, const PositionTypeDouble &, const float) const dso_internal
T sqrt(T t)
Definition: SSEVec.h:46
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
virtual std::pair< bool, double > pathLength(const Plane &plane)
PositionTypeDouble positionInDouble(double s) const
double f[11][100]
PositionTypeDouble positionInDouble(double s) const
Basic3DVector< float > PositionType
the helix is passed to the constructor and does not appear in the interface
Basic3DVector< float > DirectionType
bool first
Definition: L1TdeRCT.cc:94
double p2[4]
Definition: TauolaWrapper.h:90
Basic3DVector< double > DirectionTypeDouble
dbl *** dir
Definition: mlp_gen.cc:35
Definition: DDAxes.h:10
HelixArbitraryPlaneCrossing(const PositionType &point, const DirectionType &direction, const float curvature, const PropagationDirection propDir=alongMomentum)
const PositionType & position() const
*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