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>
7 
9  const DirectionType& direction,
10  const float curvature,
11  const PropagationDirection propDir) :
12  theQuadraticCrossingFromStart(point,direction,curvature,propDir),
13  theX0(point.x()),
14  theY0(point.y()),
15  theZ0(point.z()),
16  theRho(curvature),
17  thePropDir(propDir),
18  theCachedS(0),
19  theCachedDPhi(0.),
20  theCachedSDPhi(0.),
21  theCachedCDPhi(1.)
22 {
23  //
24  // Components of direction vector (with correct normalisation)
25  //
26  double px = direction.x();
27  double py = direction.y();
28  double pz = direction.z();
29  double pt2 = px*px+py*py;
30  double p2 = pt2+pz*pz;
31  double pI = 1./sqrt(p2);
32  double ptI = 1./sqrt(pt2);
33  theCosPhi0 = px*ptI;
34  theSinPhi0 = py*ptI;
35  theCosTheta = pz*pI;
36  theSinTheta = pt2*ptI*pI;
37 }
38 //
39 // Propagation status and path length to intersection
40 //
41 std::pair<bool,double>
43  //
44  // Constants used for control of convergence
45  //
46  const int maxIterations(100);
47  //
48  // maximum distance to plane (taking into account numerical precision)
49  //
50  float maxNumDz = theNumericalPrecision*plane.position().mag();
51  float safeMaxDist = (theMaxDistToPlane>maxNumDz?theMaxDistToPlane:maxNumDz);
52  //
53  // Prepare internal value of the propagation direction and position / direction vectors for iteration
54  //
58  //
59  // Prepare iterations: count and total pathlength
60  //
61  int iteration(maxIterations);
62  double dSTotal(0.);
63  //
64  bool first = true;
65  while ( notAtSurface(plane,xnew,safeMaxDist) ) {
66  //
67  // return empty solution vector if no convergence after maxIterations iterations
68  //
69  if ( --iteration<0 ) {
70  edm::LogInfo("HelixArbitraryPlaneCrossing") << "pathLength : no convergence";
71  return std::pair<bool,double>(false,0);
72  }
73 
74  //
75  // Use existing 2nd order object at first pass, create temporary object
76  // for subsequent passes.
77  //
78  std::pair<bool,double> deltaS2;
79  if ( first ) {
80  first = false;
82  }
83  else {
84  HelixArbitraryPlaneCrossing2Order quadraticCrossing(xnew.x(),xnew.y(),xnew.z(),
85  pnew.x(),pnew.y(),
87  theRho,
88  anyDirection);
89 
90  deltaS2 = quadraticCrossing.pathLength(plane);
91  }
92 
93  if ( !deltaS2.first ) return deltaS2;
94  //
95  // Calculate and sort total pathlength (max. 2 solutions)
96  //
97  dSTotal += deltaS2.second;
99  if ( propDir == anyDirection ) {
100  propDir = newDir;
101  }
102  else {
103  if ( newDir!=propDir ) return std::pair<bool,double>(false,0);
104  }
105  //
106  // Step forward by dSTotal.
107  //
108  xnew = positionInDouble(dSTotal);
109  pnew = directionInDouble(dSTotal);
110  }
111  //
112  // Return result
113  //
114  return std::pair<bool,double>(true,dSTotal);
115 }
116 //
117 // Position on helix after a step of path length s.
118 //
121  // use result in double precision
123  return PositionType(pos.x(),pos.y(),pos.z());
124 }
125 //
126 // Position on helix after a step of path length s in double precision.
127 //
130  //
131  // Calculate delta phi (if not already available)
132  //
133  if ( s!=theCachedS ) {
134  theCachedS = s;
138  }
139  //
140  // Calculate with appropriate formulation of full helix formula or with
141  // 2nd order approximation.
142  //
143 // if ( fabs(theCachedDPhi)>1.e-1 ) {
144  if ( std::abs(theCachedDPhi)>1.e-4 ) {
145  // "standard" helix formula
146  double o = 1./theRho;
148  theY0+( theCosPhi0*(1.-theCachedCDPhi)+theSinPhi0*theCachedSDPhi)*o,
150  }
151 // else if ( fabs(theCachedDPhi)>theNumericalPrecision ) {
152 // // full helix formula, but avoiding (1-cos(deltaPhi)) for small angles
153 // return PositionTypeDouble(theX0+(-theSinPhi0*theCachedSDPhi*theCachedSDPhi/(1.+theCachedCDPhi)+
154 // theCosPhi0*theCachedSDPhi)/theRho,
155 // theY0+(theCosPhi0*theCachedSDPhi*theCachedSDPhi/(1.+theCachedCDPhi)+
156 // theSinPhi0*theCachedSDPhi)/theRho,
157 // theZ0+theCachedS*theCosTheta);
158 // }
159  else {
160  // Use 2nd order.
162  }
163 }
164 //
165 // Direction vector on helix after a step of path length s.
166 //
169  // use result in double precision
171  return DirectionType(dir.x(),dir.y(),dir.z());
172 }
173 //
174 // Direction vector on helix after a step of path length s in double precision.
175 //
178  //
179  // Calculate delta phi (if not already available)
180  //
181  if ( s!=theCachedS ) {
182  theCachedS = s;
186  }
187 
188  if ( std::abs(theCachedDPhi)>1.e-4 ) {
189  // full helix formula
191  theSinPhi0*theCachedCDPhi+theCosPhi0*theCachedSDPhi,
193  }
194  else {
195  // 2nd order
197  }
198 }
199 // Iteration control: continue if distance to plane > theMaxDistToPlane. Includes
200 // protection for numerical precision (Surfaces work with single precision).
202  const PositionTypeDouble& point,
203  const float maxDist) const {
204  float dz = plane.localZ(Surface::GlobalPoint(point.x(),point.y(),point.z()));
205  return std::abs(dz)>maxDist;
206 }
207 
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
T curvature(T InversePt, const edm::EventSetup &iSetup)
tuple iteration
Definition: align_cfg.py:5
Definition: DDAxes.h:10
T z() const
Cartesian z coordinate.
T sqrt(T t)
Definition: SSEVec.h:28
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:79
double p2[4]
Definition: TauolaWrapper.h:90
Basic3DVector< double > DirectionTypeDouble
string s
Definition: asciidump.py:422
dbl *** dir
Definition: mlp_gen.cc:35
bool notAtSurface(const Plane &, const PositionTypeDouble &, const float) const
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