CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
RKPropagatorInS.cc
Go to the documentation of this file.
21 
23 
24 
25 std::pair<TrajectoryStateOnSurface,double>
27  const Plane& plane) const
28 {
30  if (!gp) return TsosWP(TrajectoryStateOnSurface(),0.);
31  else {
35  return errorprop( fts, plane, side, gp.parameters(),gp.s());
36  }
37 }
38 
39 std::pair< TrajectoryStateOnSurface, double>
41 {
43  if (!gp) return TsosWP(TrajectoryStateOnSurface(),0.);
44  else {
48  return errorprop( fts, cyl, side, gp.parameters(),gp.s());
49  }
50 
51 }
52 
55  const Plane& plane) const
56 {
57 
58  GlobalPoint gpos( ts.position());
59  GlobalVector gmom( ts.momentum());
60  double startZ = plane.localZ(gpos);
61  // (transverse) curvature
62  double rho = ts.transverseCurvature();
63  //
64  // Straight line approximation? |rho|<1.e-10 equivalent to ~ 1um
65  // difference in transversal position at 10m.
66  //
67  if( fabs(rho)<1.e-10 ) {
68  //
69  // Instantiate auxiliary object for finding intersection.
70  // Frame-independant point and vector are created explicitely to
71  // avoid confusing gcc (refuses to compile with temporary objects
72  // in the constructor).
73  //
74  LogDebug("RKPropagatorInS")<<" startZ = "<<startZ;
75 
76  if (fabs(startZ) < 1e-5){
77  LogDebug("RKPropagatorInS")<< "Propagation is not performed: state is already on final surface.";
78  GlobalTrajectoryParameters res( gpos, gmom, ts.charge(),
79  theVolume);
80  return GlobalParametersWithPath( res, 0.0);
81  }
82 
85  StraightLinePlaneCrossing planeCrossing(pos,dir, propagationDirection());
86  //
87  // get solution
88  //
89  std::pair<bool,double> propResult = planeCrossing.pathLength(plane);
90  if ( propResult.first && theVolume !=0) {
91  double s = propResult.second;
92  // point (reconverted to GlobalPoint)
93  GlobalPoint x (planeCrossing.position(s));
94  GlobalTrajectoryParameters res( x, gmom, ts.charge(), theVolume);
95  return GlobalParametersWithPath( res, s);
96  } else {
97  //do someting
98  LogDebug("RKPropagatorInS")<< "Straight line propgation to plane failed !!";
99  return GlobalParametersWithPath( );
100  }
101  }
102 
103 
104 #ifdef EDM_LM_DEBUG
105  if (theVolume != 0) {
106  LogDebug("RKPropagatorInS") << "RKPropagatorInS: starting prop to plane in volume with pos " << theVolume->position()
107  << " Z axis " << theVolume->toGlobal( LocalVector(0,0,1)) ;
108 
109  LogDebug("RKPropagatorInS") << "The starting position is " << ts.position() << " (global) "
110  << theVolume->toLocal(ts.position()) << " (local) " ;
111 
112  FrameChanger changer;
113  FrameChanger::PlanePtr localPlane = changer.transformPlane( plane, *theVolume);
114  LogDebug("RKPropagatorInS") << "The plane position is " << plane.position() << " (global) "
115  << localPlane->position() << " (local) " ;
116 
117  LogDebug("RKPropagatorInS") << "The initial distance to plane is " << plane.localZ( ts.position()) ;
118 
120  std::pair<bool,double> res3 = cross.pathLength(plane);
121  LogDebug("RKPropagatorInS") << "straight line distance " << res3.first << " " << res3.second ;
122  }
123 #endif
124 
126  typedef Solver::Vector RKVector;
127 
129  PathToPlane2Order pathLength( field, &field.frame());
130  CartesianLorentzForce deriv(field, ts.charge());
131 
132  RKCartesianDistance dist;
133  double eps = theTolerance;
134  Solver solver;
135  double stot = 0;
136  PropagationDirection currentDirection = propagationDirection();
137 
138  // in magVolume frame
139  RKVector start( CartesianStateAdaptor::rkstate( rkPosition(gpos), rkMomentum(gmom)));
140  int safeGuard = 0;
141  while (safeGuard++<100) {
142  CartesianStateAdaptor startState(start);
143 
144  std::pair<bool,double> path = pathLength( plane, startState.position(),
145  startState.momentum(),
146  (double) ts.charge(), currentDirection);
147  if (!path.first) {
148  LogDebug("RKPropagatorInS") << "RKPropagatorInS: Path length calculation to plane failed!"
149  << "...distance to plane " << plane.localZ( globalPosition(startState.position()))
150  << "...Local starting position in volume " << startState.position()
151  << "...Magnetic field " << field.inTesla( startState.position()) ;
152 
153 
154  return GlobalParametersWithPath();
155  }
156  else {
157  LogDebug("RKPropagatorInS") << "RKPropagatorInS: Path lenght to plane is " << path.second ;
158  }
159 
160  double sstep = path.second;
161  if ( std::abs(sstep) < eps) {
162  LogDebug("RKPropagatorInS") << "On-surface accuracy not reached, but pathLength calculation says we are there! "
163  << "path " << path.second << " distance to plane is " << startZ ;
164  GlobalTrajectoryParameters res( gtpFromVolumeLocal( startState, ts.charge()));
165  return GlobalParametersWithPath( res, stot);
166  }
167 
168  LogDebug("RKPropagatorInS") << "RKPropagatorInS: Solving for " << sstep
169  << " current distance to plane is " << startZ ;
170 
171  RKVector rkresult = solver( 0, start, sstep, deriv, dist, eps);
172  stot += sstep;
173  CartesianStateAdaptor cur( rkresult);
174  double remainingZ = plane.localZ( globalPosition(cur.position()));
175 
176  if ( fabs(remainingZ) < eps) {
177  LogDebug("RKPropagatorInS") << "On-surface accuracy reached! " << remainingZ ;
179  return GlobalParametersWithPath( res, stot);
180  }
181 
182  start = rkresult;
183 
184  if (remainingZ * startZ > 0) {
185  LogDebug("RKPropagatorInS") << "Accuracy not reached yet, trying in same direction again "
186  << remainingZ ;
187  }
188  else {
189  LogDebug("RKPropagatorInS") << "Accuracy not reached yet, trying in opposite direction "
190  << remainingZ ;
191  currentDirection = invertDirection( currentDirection);
192  }
193  startZ = remainingZ;
194  }
195  edm::LogError("FailedPropagation") << " too many iterations trying to reach plane ";
196  return GlobalParametersWithPath();
197 }
198 
201  const Cylinder& cyl) const
202 {
204  typedef Solver::Vector RKVector;
205 
206 
207  GlobalPoint sp = cyl.position();
208  if (sp.x()!=0. || sp.y()!=0.) {
209  throw PropagationException("Cannot propagate to an arbitrary cylinder");
210  }
211 
212  GlobalPoint gpos( ts.position());
213  GlobalVector gmom( ts.momentum());
214  LocalPoint pos(cyl.toLocal(gpos));
215  LocalVector mom(cyl.toLocal(gmom));
216  double startR = cyl.radius() - pos.perp();
217 
218  // LogDebug("RKPropagatorInS") << "RKPropagatorInS: starting from FTS " << ts ;
219 
220 
221  // (transverse) curvature
222  double rho = ts.transverseCurvature();
223  //
224  // Straight line approximation? |rho|<1.e-10 equivalent to ~ 1um
225  // difference in transversal position at 10m.
226  //
227  if( fabs(rho)<1.e-10 ) {
228  //
229  // Instantiate auxiliary object for finding intersection.
230  // Frame-independant point and vector are created explicitely to
231  // avoid confusing gcc (refuses to compile with temporary objects
232  // in the constructor).
233  //
234 
236 
237  //
238  // get solution
239  //
240  std::pair<bool,double> propResult = cylCrossing.pathLength(cyl);
241  if ( propResult.first && theVolume !=0) {
242  double s = propResult.second;
243  // point (reconverted to GlobalPoint)
244  GlobalPoint x (cylCrossing.position(s));
245  GlobalTrajectoryParameters res( x, gmom, ts.charge(), theVolume);
246  LogDebug("RKPropagatorInS") << "Straight line propagation to cylinder succeeded !!";
247  return GlobalParametersWithPath( res, s);
248  } else {
249  //do someting
250  edm::LogError("RKPropagatorInS") <<"Straight line propagation to cylinder failed !!";
251  return GlobalParametersWithPath( );
252  }
253  }
254 
255 
256  RKLocalFieldProvider field( fieldProvider(cyl));
257  // StraightLineCylinderCrossing pathLength( pos, mom, propagationDirection());
258  CartesianLorentzForce deriv(field, ts.charge());
259 
260  RKCartesianDistance dist;
261  double eps = theTolerance;
262  Solver solver;
263  double stot = 0;
264  PropagationDirection currentDirection = propagationDirection();
265 
266  RKVector start( CartesianStateAdaptor::rkstate( pos.basicVector(), mom.basicVector()));
267  int safeGuard = 0;
268  while (safeGuard++<100) {
269  CartesianStateAdaptor startState(start);
270  StraightLineCylinderCrossing pathLength( LocalPoint(startState.position()),
271  LocalVector(startState.momentum()),
272  currentDirection, eps);
273 
274  std::pair<bool,double> path = pathLength.pathLength( cyl);
275  if (!path.first) {
276  LogDebug("RKPropagatorInS") << "RKPropagatorInS: Path length calculation to cylinder failed!"
277  << "Radius " << cyl.radius() << " pos.perp() " << LocalPoint(startState.position()).perp() ;
278  return GlobalParametersWithPath();
279  }
280  else {
281  LogDebug("RKPropagatorInS") << "RKPropagatorInS: Path lenght to cylinder is " << path.second
282  << " from point (R,z) " << startState.position().perp()
283  << ", " << startState.position().z()
284  << " to R " << cyl.radius()
285  ;
286  }
287 
288  double sstep = path.second;
289  if ( std::abs(sstep) < eps) {
290  LogDebug("RKPropagatorInS") << "accuracy not reached, but pathLength calculation says we are there! "
291  << path.second ;
292 
294  startState.momentum(),
295  ts.charge(), cyl));
296  return GlobalParametersWithPath( res, stot);
297  }
298 
299  LogDebug("RKPropagatorInS") << "RKPropagatorInS: Solving for " << sstep
300  << " current distance to cylinder is " << startR ;
301 
302  RKVector rkresult = solver( 0, start, sstep, deriv, dist, eps);
303  stot += sstep;
304  CartesianStateAdaptor cur( rkresult);
305  double remainingR = cyl.radius() - cur.position().perp();
306  if ( fabs(remainingR) < eps) {
307  LogDebug("RKPropagatorInS") << "Accuracy reached! " << remainingR ;
309  cur.momentum(),
310  ts.charge(), cyl));
311  return GlobalParametersWithPath( res, stot);
312  }
313 
314  start = rkresult;
315  if (remainingR * startR > 0) {
316  LogDebug("RKPropagatorInS") << "Accuracy not reached yet, trying in same direction again "
317  << remainingR ;
318  }
319  else {
320  LogDebug("RKPropagatorInS") << "Accuracy not reached yet, trying in opposite direction "
321  << remainingR ;
322  currentDirection = invertDirection( currentDirection);
323  }
324  startR = remainingR;
325  }
326  edm::LogError("FailedPropagation") << " too many iterations trying to reach cylinder ";
327  return GlobalParametersWithPath();
328 }
329 
331 RKPropagatorInS::propagate(const FreeTrajectoryState& fts, const Plane& plane) const
332 {
333  return propagateWithPath( fts, plane).first;
334 }
335 
338 {
339  return propagateWithPath( fts, cyl).first;
340 }
341 
343 {
344  return new RKPropagatorInS(*this);
345 }
346 
348  const Basic3DVector<double>& lmom,
349  TrackCharge ch, const Surface& surf) const
350 {
351  return GlobalTrajectoryParameters( surf.toGlobal( LocalPoint( lpos)),
352  surf.toGlobal( LocalVector( lmom)), ch, theVolume);
353 }
354 
356 {
358 }
359 
361 {
362  return RKLocalFieldProvider( *theVolume, cyl);
363 }
364 
366 {
367  if (dir == anyDirection) return dir;
368  return ( dir == alongMomentum ? oppositeToMomentum : alongMomentum);
369 }
370 
372 {
373  if (theVolume != 0) return theVolume->toLocal( pos).basicVector();
374  else return pos.basicVector();
375 }
376 
378 {
379  if (theVolume != 0) return theVolume->toLocal( mom).basicVector();
380  else return mom.basicVector();
381 }
382 
384 {
385  if (theVolume != 0) return theVolume->toGlobal( LocalPoint(pos));
386  else return GlobalPoint(pos);
387 }
388 
390 
391 {
392  if (theVolume != 0) return theVolume->toGlobal( LocalVector(mom));
393  else return GlobalVector(mom);
394 }
395 
398  TrackCharge charge) const
399 {
401  globalMomentum(state.momentum()),
402  charge, theVolume);
403 }
#define LogDebug(id)
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
Definition: Surface.h:78
virtual Propagator * clone() const
Local3DVector LocalVector
Definition: LocalVector.h:12
GlobalTrajectoryParameters gtpFromVolumeLocal(const CartesianStateAdaptor &state, TrackCharge charge) const
virtual PropagationDirection propagationDirection() const
Definition: Propagator.h:143
ROOT::Math::Plane3D::Vector Vector
Definition: EcalHitMaker.cc:28
GlobalTrajectoryParameters gtpFromLocal(const Basic3DVector< double > &lpos, const Basic3DVector< double > &lmom, TrackCharge ch, const Surface &surf) const
GlobalParametersWithPath propagateParametersOnCylinder(const FreeTrajectoryState &ts, const Cylinder &cyl) const
T perp() const
Magnitude of transverse component.
Definition: DDAxes.h:10
RKLocalFieldProvider fieldProvider() const
std::pair< bool, double > pathLength(const Cylinder &cyl) const
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
Basic3DVector cross(const Basic3DVector &v) const
Vector product, or &quot;cross&quot; product, with a vector of same type.
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
#define abs(x)
Definition: mlp_lapack.h:159
PropagationDirection
const GlobalTrajectoryParameters & parameters() const
TrackCharge charge() const
double charge(const std::vector< uint8_t > &Ampls)
Definition: Plane.h:17
int path() const
Definition: HLTadd.h:3
int TrackCharge
Definition: TrackCharge.h:4
Basic3DVector< double > rkMomentum(const GlobalVector &mom) const
std::pair< TrajectoryStateOnSurface, double > TsosWP
Estimator of the distance between two state vectors, e.g. for convergence test.
Scalar radius() const
Radius of the cylinder.
Definition: Cylinder.h:55
T z() const
Cartesian z coordinate.
PropagationDirection invertDirection(PropagationDirection dir) const
LocalPoint toLocal(const GlobalPoint &gp) const
virtual std::pair< TrajectoryStateOnSurface, double > propagateWithPath(const FreeTrajectoryState &, const Plane &) const
Derivative calculation for the 6D cartesian case.
static RKSmallVector< double, 6 > rkstate(const Vector3D &pos, const Vector3D &mom)
const Vector3D & momentum() const
virtual TrajectoryStateOnSurface propagate(const FreeTrajectoryState &, const Plane &) const
const MagVolume * theVolume
T perp() const
Magnitude of transverse component.
Common base class.
GlobalVector momentum() const
std::pair< bool, double > pathLength(const Plane &plane) const
PositionType position(double s) const
GlobalPoint toGlobal(const LocalPoint &lp) const
GlobalPoint position() const
const Frame & frame() const
The reference frame in which the field is defined.
static Plane transformPlane(const Plane &plane, const GloballyPositioned< T > &frame)
Definition: FrameChanger.h:14
double transverseCurvature() const
char state
Definition: procUtils.cc:75
GlobalPoint globalPosition(const Basic3DVector< double > &pos) const
GlobalParametersWithPath propagateParametersOnPlane(const FreeTrajectoryState &ts, const Plane &plane) const
Local3DPoint LocalPoint
Definition: LocalPoint.h:11
RKPropagatorInS(const MagVolume &vol, PropagationDirection dir=alongMomentum, double tolerance=5.e-5)
string s
Definition: asciidump.py:422
dbl *** dir
Definition: mlp_gen.cc:35
const Vector3D & position() const
std::pair< bool, double > pathLength(const Cylinder &cyl) const
T x() const
Definition: PV3DBase.h:56
const BasicVectorType & basicVector() const
Definition: PV3DBase.h:54
const PositionType & position() const
GlobalVector globalMomentum(const Basic3DVector< double > &mom) const
Basic3DVector< double > rkPosition(const GlobalPoint &pos) const
Global3DVector GlobalVector
Definition: GlobalVector.h:10
Vector inTesla(const LocalPoint &lp) const
the argument lp is in the local frame specified in the constructor