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 
101  //try {
102  ReferenceCountingPointer<TangentPlane> plane(cylinder.tangentPlane(x)); // need to be here until tsos is created!
103  return propagatedStateWithPath(fts,*plane,gtp,s);
104  /*
105  } catch(...) {
106  std::cout << "wrong tangent to cylinder " << x
107  << " pos, rad " << cylinder.position() << " " << cylinder.radius()
108  << std::endl;
109  return TsosWP(TrajectoryStateOnSurface(),0.);
110  }
111  */
112 }
113 
114 std::pair<TrajectoryStateOnSurface,double>
116  const Surface& surface,
117  const GlobalTrajectoryParameters& gtp,
118  const double& s) const
119 {
120  //
121  // for forward propagation: state is before surface,
122  // for backward propagation: state is after surface
123  //
124  SurfaceSide side = PropagationDirectionFromPath()(s,propagationDirection())==alongMomentum
126  //
127  //
128  // error propagation (if needed) and conversion to a TrajectoryStateOnSurface
129  //
130  if (fts.hasError()) {
131  //
132  // compute jacobian
133  //
134  AnalyticalCurvilinearJacobian analyticalJacobian(fts.parameters(), gtp.position(), gtp.momentum(), s);
135  const AlgebraicMatrix55 &jacobian = analyticalJacobian.jacobian();
136  CurvilinearTrajectoryError cte(ROOT::Math::Similarity(jacobian, fts.curvilinearError().matrix()));
137  return TsosWP(TrajectoryStateOnSurface(gtp,cte,surface,side),s);
138  }
139  else {
140  //
141  // return state without errors
142  //
143  return TsosWP(TrajectoryStateOnSurface(gtp,surface,side),s);
144  }
145 }
146 
148  const FreeTrajectoryState& fts, const Cylinder& cylinder,
149  GlobalPoint& x, GlobalVector& p, double& s) const
150 {
151 
152  GlobalPoint const & sp = cylinder.position();
153  if (sp.x()!=0. || sp.y()!=0.) {
154  throw PropagationException("Cannot propagate to an arbitrary cylinder");
155  }
156  // preset output
157  x = fts.position();
158  p = fts.momentum();
159  s = 0;
160  // (transverse) curvature
161  double rho = fts.transverseCurvature();
162  //
163  // Straight line approximation? |rho|<1.e-10 equivalent to ~ 1um
164  // difference in transversal position at 10m.
165  //
166  if( fabs(rho)<1.e-10 )
167  return propagateWithLineCrossing(fts.position(),p,cylinder,x,s);
168  //
169  // Helix case
170  //
171  // check for possible intersection
172  const double tolerance = 1.e-4; // 1 micron distance
173  double rdiff = x.perp() - cylinder.radius();
174  if ( fabs(rdiff) < tolerance ) return true;
175  //
176  // Instantiate HelixBarrelCylinderCrossing and get solutions
177  //
178  HelixBarrelCylinderCrossing cylinderCrossing(fts.position(),fts.momentum(),rho,
179  propagationDirection(),cylinder);
180  if ( !cylinderCrossing.hasSolution() ) return false;
181  // path length
182  s = cylinderCrossing.pathLength();
183  // point
184  x = cylinderCrossing.position();
185  // direction (renormalised)
186  p = cylinderCrossing.direction().unit()*fts.momentum().mag();
187  return true;
188 }
189 
190 bool
192  const Plane& plane,
193  GlobalPoint& x,
194  GlobalVector& p,
195  double& s) const
196 {
197  // initialisation of position, momentum and path length
198  x = fts.position();
199  p = fts.momentum();
200  s = 0;
201  // (transverse) curvature
202  double rho = fts.transverseCurvature();
203  //
204  // Straight line approximation? |rho|<1.e-10 equivalent to ~ 1um
205  // difference in transversal position at 10m.
206  //
207  if( fabs(rho)<1.e-10 )
208  return propagateWithLineCrossing(fts.position(),p,plane,x,s);
209  //
210  // Helix case
211  //
212  GlobalVector u = plane.normalVector();
213  const double small = 1.e-6; // for orientation of planes
214  //
215  // Frame-independant point and vector are created explicitely to
216  // avoid confusing gcc (refuses to compile with temporary objects
217  // in the constructor).
218  //
221  if (isOldPropagationType && fabs(u.z()) < small) {
222  // barrel plane:
223  // instantiate HelixBarrelPlaneCrossing, get vector of solutions and check for existance
224  HelixBarrelPlaneCrossingByCircle planeCrossing(helixPos,helixDir,rho,propagationDirection());
225  return propagateWithHelixCrossing(planeCrossing,plane,fts.momentum().mag(),x,p,s);
226  }
227  if (isOldPropagationType && fabs(u.x()) < small && fabs(u.y()) < small) {
228  // forward plane:
229  // instantiate HelixForwardPlaneCrossing, get vector of solutions and check for existance
230  HelixForwardPlaneCrossing planeCrossing(helixPos,helixDir,rho,propagationDirection());
231  return propagateWithHelixCrossing(planeCrossing,plane,fts.momentum().mag(),x,p,s);
232  }
233  else {
234  // arbitrary plane:
235  // instantiate HelixArbitraryPlaneCrossing, get vector of solutions and check for existance
236  if(isOldPropagationType){
237  HelixArbitraryPlaneCrossing planeCrossing(helixPos,helixDir,rho,propagationDirection());
238  return propagateWithHelixCrossing(planeCrossing,plane,fts.momentum().mag(),x,p,s);
239  }else{
240  //--- Alternative implementation to be used for the propagation of the parameters of looping
241  // particles that cross twice the (infinite) surface of the plane. It is not trivial to determine
242  // which of the two intersections has to be returned.
243 
244  //---- FIXME: WHAT FOLLOWS HAS TO BE REWRITTEN IN A CLEANER (AND CPU-OPTIMIZED) WAY ---------
245  LogDebug("AnalyticalPropagator") << "In AnaliticalProp, calling HAPC " << "\n"
246  << "plane is centered in xyz: "
247  << plane.position().x() << " , "
248  << plane.position().y() << " , "
249  << plane.position().z() << "\n";
250 
251 
252  GlobalPoint gp1 = fts.position();
253  GlobalVector gm1 = fts.momentum();
254  double s1 = 0;
255  double rho1 = fts.transverseCurvature();
256  HelixPlaneCrossing::PositionType helixPos1(gp1);
257  HelixPlaneCrossing::DirectionType helixDir1(gm1);
258  LogDebug("AnalyticalPropagator") << "gp1 before calling planeCrossing1: " << gp1 << "\n";
259  HelixArbitraryPlaneCrossing planeCrossing1(helixPos1,helixDir1,rho1,propagationDirection());
260 
263 
264  double tolerance(0.0050);
265  if(propagationDirection()==oppositeToMomentum)
266  tolerance *=-1;
267 
268  bool check1 = propagateWithHelixCrossing(planeCrossing1,plane,fts.momentum().mag(),gp1,gm1,s1);
269  double dphi1 = fabs(fts.momentum().phi()-gm1.phi());
270  LogDebug("AnalyticalPropagator") << "check1, s1, dphi, gp1: "
271  << check1 << " , "
272  << s1 << " , "
273  << dphi1 << " , "
274  << gp1 << "\n";
275 
276  //move forward a bit to avoid that the propagator doesn't propagate because the state is already on surface.
277  //we want to go to the other point of intersection between the helix and the plane
278  xGen = planeCrossing1.position(s1+tolerance);
279  pGen = planeCrossing1.direction(s1+tolerance);
280 
281  /*
282  if(!check1 || s1>170 ){
283  //PropagationDirection newDir = (propagationDirection() == alongMomentum) ? oppositeToMomentum : alongMomentum;
284  PropagationDirection newDir = anyDirection;
285  HelixArbitraryPlaneCrossing planeCrossing1B(helixPos1,helixDir1,rho1,newDir);
286  check1 = propagateWithHelixCrossing(planeCrossing1B,plane,fts.momentum().mag(),gp1,gm1,s1);
287  LogDebug("AnalyticalPropagator") << "after second attempt, check1, s1,gp1: "
288  << check1 << " , "
289  << s1 << " , " << gp1 << "\n";
290 
291  xGen = planeCrossing1B.position(s1+tolerance);
292  pGen = planeCrossing1B.direction(s1+tolerance);
293  }
294  */
295 
296  if(!check1){
297  LogDebug("AnalyticalPropagator") << "failed also second attempt. No idea what to do, then bailout" << "\n";
298  }
299 
300 
301  pGen *= gm1.mag()/pGen.mag();
302  GlobalPoint gp2(xGen);
303  GlobalVector gm2(pGen);
304  double s2 = 0;
305  double rho2 = rho1;
306  HelixPlaneCrossing::PositionType helixPos2(gp2);
307  HelixPlaneCrossing::DirectionType helixDir2(gm2);
308  HelixArbitraryPlaneCrossing planeCrossing2(helixPos2,helixDir2,rho2,propagationDirection());
309 
310  bool check2 = propagateWithHelixCrossing(planeCrossing2,plane,gm2.mag(),gp2,gm2,s2);
311 
312  if(!check2){
313  x = gp1;
314  p = gm1;
315  s = s1;
316  return check1;
317  }
318 
319  if(!check1){
320  edm::LogError("AnalyticalPropagator") << "LOGIC ERROR: I should not have entered here!" << "\n";
321  return false;
322  }
323 
324 
325  LogDebug("AnalyticalPropagator") << "check2, s2, gp2: "
326  << check2 << " , "
327  << s2 << " , " << gp2 << "\n";
328 
329 
330  double dist1 = (plane.position()-gp1).perp();
331  double dist2 = (plane.position()-gp2).perp();
332 
333 
334  LogDebug("AnalyticalPropagator") << "propDir, dist1, dist2: "
335  << propagationDirection() << " , "
336  << dist1 << " , "
337  << dist2 << "\n";
338 
339  //If there are two solutions, the one which is the closest to the module's center is chosen
340  if(dist1<2*dist2){
341  x = gp1;
342  p = gm1;
343  s = s1;
344  return check1;
345  }else if(dist2<2*dist1){
346  x = gp2;
347  p = gm2;
348  s = s1+s2+tolerance;
349  return check2;
350  }else{
351  if(fabs(s1)<fabs(s2)){
352  x = gp1;
353  p = gm1;
354  s = s1;
355  return check1;
356  }else{
357  x = gp2;
358  p = gm2;
359  s = s1+s2+tolerance;
360  return check2;
361  }
362  }
363 
364  //-------- END of ugly piece of code ---------------
365  }
366 
367  }
368 }
369 
370 bool
372  const GlobalVector& p0,
373  const Plane& plane,
374  GlobalPoint& x, double& s) const {
375  //
376  // Instantiate auxiliary object for finding intersection.
377  // Frame-independant point and vector are created explicitely to
378  // avoid confusing gcc (refuses to compile with temporary objects
379  // in the constructor).
380  //
383  StraightLinePlaneCrossing planeCrossing(pos,dir,propagationDirection());
384  //
385  // get solution
386  //
387  std::pair<bool,double> propResult = planeCrossing.pathLength(plane);
388  if ( !propResult.first ) return false;
389  s = propResult.second;
390  // point (reconverted to GlobalPoint)
391  x = GlobalPoint(planeCrossing.position(s));
392  //
393  return true;
394 }
395 
396 bool
398  const GlobalVector& p0,
399  const Cylinder& cylinder,
400  GlobalPoint& x, double& s) const {
401  //
402  // Instantiate auxiliary object for finding intersection.
403  // Frame-independant point and vector are created explicitely to
404  // avoid confusing gcc (refuses to compile with temporary objects
405  // in the constructor).
406  //
407  StraightLineBarrelCylinderCrossing cylCrossing(x0,p0,propagationDirection());
408  //
409  // get solution
410  //
411  std::pair<bool,double> propResult = cylCrossing.pathLength(cylinder);
412  if ( !propResult.first ) return false;
413  s = propResult.second;
414  // point (reconverted to GlobalPoint)
415  x = cylCrossing.position(s);
416  //
417  return true;
418 }
419 
420 bool
422  const Plane& plane,
423  const float pmag,
424  GlobalPoint& x,
425  GlobalVector& p,
426  double& s) const {
427  // get solution
428  std::pair<bool,double> propResult = planeCrossing.pathLength(plane);
429  if ( !propResult.first ) return false;
430  s = propResult.second;
431  // point (reconverted to GlobalPoint)
432  HelixPlaneCrossing::PositionType xGen = planeCrossing.position(s);
433  x = GlobalPoint(xGen);
434  // direction (reconverted to GlobalVector, renormalised)
435  HelixPlaneCrossing::DirectionType pGen = planeCrossing.direction(s);
436  pGen *= pmag/pGen.mag();
437  p = GlobalVector(pGen);
438  //
439  return true;
440 }
#define LogDebug(id)
T mag2() const
Definition: PV3DBase.h:65
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
bool propagateParametersOnPlane(const FreeTrajectoryState &fts, const Plane &plane, GlobalPoint &x, GlobalVector &p, double &s) const dso_internal
parameter propagation to plane (returns position, momentum and path length)
virtual ReferenceCountingPointer< TangentPlane > tangentPlane(const GlobalPoint &) const
tangent plane to surface from global point
Definition: Cylinder.cc:23
std::pair< TrajectoryStateOnSurface, double > propagatedStateWithPath(const FreeTrajectoryState &fts, const Surface &surface, const GlobalTrajectoryParameters &gtp, const double &s) const dso_internal
propagation of errors (if needed) and generation of a new TSOS
const GlobalTrajectoryParameters & parameters() const
const AlgebraicMatrix55 & jacobian() const
virtual DirectionType direction(double s) const =0
Definition: DDAxes.h:10
Geom::Phi< T > phi() const
Definition: PV3DBase.h:68
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
T y() const
Definition: PV3DBase.h:62
float localZ(const GlobalPoint &gp) const
Fast access to distance from plane for a point.
Definition: Plane.h:52
virtual DirectionType direction(double s) const
T perp2() const
Definition: PV3DBase.h:70
TrackCharge charge() const
virtual PositionType position(double s) const
const CurvilinearTrajectoryError & curvilinearError() const
Definition: Plane.h:17
tuple s2
Definition: indexGen.py:106
T mag() const
Definition: PV3DBase.h:66
std::pair< TrajectoryStateOnSurface, double > propagateWithPath(const FreeTrajectoryState &fts, const Plane &plane) const
propagation to plane with path length
Common base class.
bool propagateParametersOnCylinder(const FreeTrajectoryState &fts, const Cylinder &cylinder, GlobalPoint &x, GlobalVector &p, double &s) const dso_internal
parameter propagation to cylinder (returns position, momentum and path length)
GlobalVector momentum() const
std::pair< bool, double > pathLength(const Plane &plane) const
PositionType position(double s) const
Vector3DBase unit() const
Definition: Vector3DBase.h:57
GlobalPoint position() const
virtual std::pair< bool, double > pathLength(const Plane &)=0
double transverseCurvature() const
virtual PositionType position(double s) const =0
T perp() const
Magnitude of transverse component.
const AlgebraicSymMatrix55 & matrix() const
bool propagateWithHelixCrossing(HelixPlaneCrossing &, const Plane &, const float, GlobalPoint &, GlobalVector &, double &s) const dso_internal
helix parameter propagation to a plane using HelixPlaneCrossing
std::pair< TrajectoryStateOnSurface, double > TsosWP
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:61
const PositionType & position() const
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepStd< double, 5, 5 > > AlgebraicMatrix55
bool propagateWithLineCrossing(const GlobalPoint &, const GlobalVector &, const Plane &, GlobalPoint &, double &) const dso_internal
straight line parameter propagation to a plane
Global3DVector GlobalVector
Definition: GlobalVector.h:10