CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
AnalyticalPropagator.cc
Go to the documentation of this file.
2 
6 
8 
20 
22 
23 #include <cmath>
24 
25 using namespace SurfaceSideDefinition;
26 
27 std::pair<TrajectoryStateOnSurface,double>
29  const Plane& plane) const
30 {
31  // check curvature
32  double rho = fts.transverseCurvature();
33 
34  // propagate parameters
35  GlobalPoint x;
37  double s;
38 
39  // check if already on plane
40  const float maxDistToPlane(0.1e-4);
41  const float numericalPrecision(5.e-7);
42  float maxDz = numericalPrecision*plane.position().mag();
43  if ( fabs(plane.localZ(fts.position()))>(maxDistToPlane>maxDz?maxDistToPlane:maxDz) ) {
44  // propagate
45  bool parametersOK = this->propagateParametersOnPlane(fts, plane, x, p, s);
46  // check status and deltaPhi limit
47  float dphi2 = s*rho;
48  dphi2 = dphi2*dphi2*fts.momentum().perp2()/fts.momentum().mag2();
49  if ( !parametersOK || dphi2>theMaxDPhi2 ) return TsosWP(TrajectoryStateOnSurface(),0.);
50  }
51  else {
52  LogDebug("AnalyticalPropagator")<<"not going anywhere. Already on surface.\n"
53  <<"plane.localZ(fts.position()): "<<plane.localZ(fts.position())<<"\n"
54  <<"maxDistToPlane: "<<maxDistToPlane<<"\n"
55  <<"maxDz: "<<maxDz<<"\n"
56  <<"plane.position().mag(): "<<plane.position().mag();
57  x = fts.position();
58  p = fts.momentum();
59  s = 0.;
60  }
61  //
62  // Compute propagated state and check change in curvature
63  //
64  GlobalTrajectoryParameters gtp(x,p,fts.charge(),theField);
65  if ( fabs(rho)>1.e-10 && fabs((gtp.transverseCurvature()-rho)/rho)>theMaxDBzRatio )
66  return TsosWP(TrajectoryStateOnSurface(),0.);
67  //
68  // construct TrajectoryStateOnSurface
69  //
70  return propagatedStateWithPath(fts,plane,gtp,s);
71 }
72 
73 
74 std::pair<TrajectoryStateOnSurface,double>
76  const Cylinder& cylinder) const
77 {
78  // check curvature
79  double rho = fts.transverseCurvature();
80 
81  // propagate parameters
82  GlobalPoint x;
84  double s = 0;
85 
86  bool parametersOK = this->propagateParametersOnCylinder(fts, cylinder, x, p, s);
87  // check status and deltaPhi limit
88  float dphi2 = s*rho;
89  dphi2 = dphi2*dphi2*fts.momentum().perp2()/fts.momentum().mag2();
90  if ( !parametersOK || dphi2>theMaxDPhi2 ) return TsosWP(TrajectoryStateOnSurface(),0.);
91  //
92  // Compute propagated state and check change in curvature
93  //
94  GlobalTrajectoryParameters gtp(x,p,fts.charge(),theField);
95  if ( fabs(rho)>1.e-10 && fabs((gtp.transverseCurvature()-rho)/rho)>theMaxDBzRatio )
96  return TsosWP(TrajectoryStateOnSurface(),0.);
97  //
98  // create result TSOS on TangentPlane (local parameters & errors are better defined)
99  //
100  ReferenceCountingPointer<TangentPlane> plane(cylinder.tangentPlane(x)); // need to be here until tsos is created!
101  return propagatedStateWithPath(fts,*plane,gtp,s);
102 }
103 
104 std::pair<TrajectoryStateOnSurface,double>
106  const Surface& surface,
107  const GlobalTrajectoryParameters& gtp,
108  const double& s) const
109 {
110  //
111  // for forward propagation: state is before surface,
112  // for backward propagation: state is after surface
113  //
114  SurfaceSide side = PropagationDirectionFromPath()(s,propagationDirection())==alongMomentum
116  //
117  //
118  // error propagation (if needed) and conversion to a TrajectoryStateOnSurface
119  //
120  if (fts.hasError()) {
121  //
122  // compute jacobian
123  //
124  AnalyticalCurvilinearJacobian analyticalJacobian(fts.parameters(), gtp.position(), gtp.momentum(), s);
125  const AlgebraicMatrix55 &jacobian = analyticalJacobian.jacobian();
126  CurvilinearTrajectoryError cte(ROOT::Math::Similarity(jacobian, fts.curvilinearError().matrix()));
127  return TsosWP(TrajectoryStateOnSurface(gtp,cte,surface,side),s);
128  }
129  else {
130  //
131  // return state without errors
132  //
133  return TsosWP(TrajectoryStateOnSurface(gtp,surface,side),s);
134  }
135 }
136 
138  const FreeTrajectoryState& fts, const Cylinder& cylinder,
139  GlobalPoint& x, GlobalVector& p, double& s) const
140 {
141 
142  GlobalPoint const & sp = cylinder.position();
143  if (sp.x()!=0. || sp.y()!=0.) {
144  throw PropagationException("Cannot propagate to an arbitrary cylinder");
145  }
146  // preset output
147  x = fts.position();
148  p = fts.momentum();
149  s = 0;
150  // (transverse) curvature
151  double rho = fts.transverseCurvature();
152  //
153  // Straight line approximation? |rho|<1.e-10 equivalent to ~ 1um
154  // difference in transversal position at 10m.
155  //
156  if( fabs(rho)<1.e-10 )
157  return propagateWithLineCrossing(fts.position(),p,cylinder,x,s);
158  //
159  // Helix case
160  //
161  // check for possible intersection
162  const double tolerance = 1.e-4; // 1 micron distance
163  double rdiff = x.perp() - cylinder.radius();
164  if ( fabs(rdiff) < tolerance ) return true;
165  //
166  // Instantiate HelixBarrelCylinderCrossing and get solutions
167  //
168  HelixBarrelCylinderCrossing cylinderCrossing(fts.position(),fts.momentum(),rho,
169  propagationDirection(),cylinder);
170  if ( !cylinderCrossing.hasSolution() ) return false;
171  // path length
172  s = cylinderCrossing.pathLength();
173  // point
174  x = cylinderCrossing.position();
175  // direction (renormalised)
176  p = cylinderCrossing.direction().unit()*fts.momentum().mag();
177  return true;
178 }
179 
180 bool
182  const Plane& plane,
183  GlobalPoint& x,
184  GlobalVector& p,
185  double& s) const
186 {
187  // initialisation of position, momentum and path length
188  x = fts.position();
189  p = fts.momentum();
190  s = 0;
191  // (transverse) curvature
192  double rho = fts.transverseCurvature();
193  //
194  // Straight line approximation? |rho|<1.e-10 equivalent to ~ 1um
195  // difference in transversal position at 10m.
196  //
197  if( fabs(rho)<1.e-10 )
198  return propagateWithLineCrossing(fts.position(),p,plane,x,s);
199  //
200  // Helix case
201  //
202  GlobalVector u = plane.normalVector();
203  const double small = 1.e-6; // for orientation of planes
204  //
205  // Frame-independant point and vector are created explicitely to
206  // avoid confusing gcc (refuses to compile with temporary objects
207  // in the constructor).
208  //
211  if (fabs(u.z()) < small) {
212  // barrel plane:
213  // instantiate HelixBarrelPlaneCrossing, get vector of solutions and check for existance
214  HelixBarrelPlaneCrossingByCircle planeCrossing(helixPos,helixDir,rho,propagationDirection());
215  return propagateWithHelixCrossing(planeCrossing,plane,fts.momentum().mag(),x,p,s);
216  }
217  if (fabs(u.x()) < small && fabs(u.y()) < small) {
218  // forward plane:
219  // instantiate HelixForwardPlaneCrossing, get vector of solutions and check for existance
220  HelixForwardPlaneCrossing planeCrossing(helixPos,helixDir,rho,propagationDirection());
221  return propagateWithHelixCrossing(planeCrossing,plane,fts.momentum().mag(),x,p,s);
222  }
223  else {
224  // arbitrary plane:
225  // instantiate HelixArbitraryPlaneCrossing, get vector of solutions and check for existance
226  HelixArbitraryPlaneCrossing planeCrossing(helixPos,helixDir,rho,propagationDirection());
227  return propagateWithHelixCrossing(planeCrossing,plane,fts.momentum().mag(),x,p,s);
228  }
229 }
230 
231 bool
233  const GlobalVector& p0,
234  const Plane& plane,
235  GlobalPoint& x, double& s) const {
236  //
237  // Instantiate auxiliary object for finding intersection.
238  // Frame-independant point and vector are created explicitely to
239  // avoid confusing gcc (refuses to compile with temporary objects
240  // in the constructor).
241  //
244  StraightLinePlaneCrossing planeCrossing(pos,dir,propagationDirection());
245  //
246  // get solution
247  //
248  std::pair<bool,double> propResult = planeCrossing.pathLength(plane);
249  if ( !propResult.first ) return false;
250  s = propResult.second;
251  // point (reconverted to GlobalPoint)
252  x = GlobalPoint(planeCrossing.position(s));
253  //
254  return true;
255 }
256 
257 bool
259  const GlobalVector& p0,
260  const Cylinder& cylinder,
261  GlobalPoint& x, double& s) const {
262  //
263  // Instantiate auxiliary object for finding intersection.
264  // Frame-independant point and vector are created explicitely to
265  // avoid confusing gcc (refuses to compile with temporary objects
266  // in the constructor).
267  //
268  StraightLineBarrelCylinderCrossing cylCrossing(x0,p0,propagationDirection());
269  //
270  // get solution
271  //
272  std::pair<bool,double> propResult = cylCrossing.pathLength(cylinder);
273  if ( !propResult.first ) return false;
274  s = propResult.second;
275  // point (reconverted to GlobalPoint)
276  x = cylCrossing.position(s);
277  //
278  return true;
279 }
280 
281 bool
283  const Plane& plane,
284  const float pmag,
285  GlobalPoint& x,
286  GlobalVector& p,
287  double& s) const {
288  // get solution
289  std::pair<bool,double> propResult = planeCrossing.pathLength(plane);
290  if ( !propResult.first ) return false;
291  s = propResult.second;
292  // point (reconverted to GlobalPoint)
293  HelixPlaneCrossing::PositionType xGen = planeCrossing.position(s);
294  x = GlobalPoint(xGen);
295  // direction (reconverted to GlobalVector, renormalised)
296  HelixPlaneCrossing::DirectionType pGen = planeCrossing.direction(s);
297  pGen *= pmag/pGen.mag();
298  p = GlobalVector(pGen);
299  //
300  return true;
301 }
#define LogDebug(id)
T mag2() const
Definition: PV3DBase.h:60
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
virtual ReferenceCountingPointer< TangentPlane > tangentPlane(const GlobalPoint &) const
tangent plane to surface from global point
Definition: Cylinder.cc:23
const GlobalTrajectoryParameters & parameters() const
virtual DirectionType direction(double s) const =0
Definition: DDAxes.h:10
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
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
bool propagateWithLineCrossing(const GlobalPoint &, const GlobalVector &, const Plane &, GlobalPoint &, double &) const
straight line parameter propagation to a plane
T perp2() const
Definition: PV3DBase.h:65
TrackCharge charge() const
const CurvilinearTrajectoryError & curvilinearError() const
Definition: Plane.h:17
bool propagateParametersOnPlane(const FreeTrajectoryState &fts, const Plane &plane, GlobalPoint &x, GlobalVector &p, double &s) const
parameter propagation to plane (returns position, momentum and path length)
std::pair< TrajectoryStateOnSurface, double > propagatedStateWithPath(const FreeTrajectoryState &fts, const Surface &surface, const GlobalTrajectoryParameters &gtp, const double &s) const
propagation of errors (if needed) and generation of a new TSOS
T mag() const
Definition: PV3DBase.h:61
std::pair< TrajectoryStateOnSurface, double > propagateWithPath(const FreeTrajectoryState &fts, const Plane &plane) const
propagation to plane with path length
bool propagateWithHelixCrossing(HelixPlaneCrossing &, const Plane &, const float, GlobalPoint &, GlobalVector &, double &s) const
helix parameter propagation to a plane using HelixPlaneCrossing
virtual const AlgebraicMatrix55 & jacobian() const
Common base class.
GlobalVector momentum() const
std::pair< bool, double > pathLength(const Plane &plane) const
PositionType position(double s) const
Vector3DBase unit() const
Definition: Vector3DBase.h:57
bool propagateParametersOnCylinder(const FreeTrajectoryState &fts, const Cylinder &cylinder, GlobalPoint &x, GlobalVector &p, double &s) const
parameter propagation to cylinder (returns position, momentum and path length)
GlobalPoint position() const
virtual std::pair< bool, double > pathLength(const Plane &)=0
double transverseCurvature() const
virtual PositionType position(double s) const =0
const AlgebraicSymMatrix55 & matrix() const
std::pair< TrajectoryStateOnSurface, double > TsosWP
string s
Definition: asciidump.py:422
dbl *** dir
Definition: mlp_gen.cc:35
Definition: DDAxes.h:10
std::pair< bool, double > pathLength(const Cylinder &cyl) const
T x() const
Definition: PV3DBase.h:56
const PositionType & position() const
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepStd< double, 5, 5 > > AlgebraicMatrix55
Global3DVector GlobalVector
Definition: GlobalVector.h:10