PatternTools/src/TrajectoryExtrapolatorToLine.cc

0001
0002 #include "TrackingTools/PatternTools/interface/TrajectoryExtrapolatorToLine.h"
0003 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
0004 #include "DataFormats/GeometryCommonDetAlgo/interface/DeepCopyPointerByClone.h"
0005 #include "DataFormats/GeometrySurface/interface/BoundPlane.h"
0006 #include "DataFormats/GeometrySurface/interface/OpenBounds.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008
0009 TrajectoryStateOnSurface TrajectoryExtrapolatorToLine::extrapolate(const FreeTrajectoryState& fts,
0010                                                                    const Line& L,
0011                                                                    const Propagator& aPropagator) const {
0012   DeepCopyPointerByClone<Propagator> p(aPropagator.clone());
0013   p->setPropagationDirection(anyDirection);
0014
0015   FreeTrajectoryState fastFts(fts.parameters(), fts.curvilinearError());
0016   GlobalVector T1 = fastFts.momentum().unit();
0017   GlobalPoint T0 = fastFts.position();
0018   double distance = 9999999.9;
0019   double old_distance;
0020   bool refining = true;
0021
0022   LogDebug("TrajectoryExtrapolatorToLine") << "START REFINING";
0023   while (refining) {
0024     LogDebug("TrajectoryExtrapolatorToLine") << "Refining cycle...";
0025     // describe orientation of target surface on basis of track parameters
0026     Line T(T0, T1);
0027     GlobalPoint B = T.closerPointToLine(L);
0028     old_distance = distance;
0029
0030     //create surface
0031     GlobalPoint BB = B + 0.3 * (T0 - B);
0032     Surface::PositionType pos(BB);
0033     GlobalVector XX(T1.y(), -T1.x(), 0.);
0034     GlobalVector YY(T1.cross(XX));
0035     Surface::RotationType rot(XX, YY);
0036     ReferenceCountingPointer<Plane> surface = Plane::build(pos, rot);
0037     LogDebug("TrajectoryExtrapolatorToLine") << "Current plane position: " << surface->toGlobal(LocalPoint(0., 0., 0.));
0038     LogDebug("TrajectoryExtrapolatorToLine") << "Current plane normal: " << surface->toGlobal(LocalVector(0, 0, 1));
0039     LogDebug("TrajectoryExtrapolatorToLine") << "Current momentum:     " << T1;
0040
0041     // extrapolate fastFts to target surface
0042     TrajectoryStateOnSurface tsos = p->propagate(fastFts, *surface);
0043
0044     if (!tsos.isValid()) {
0045       LogDebug("TrajectoryExtrapolatorToLine") << "TETL - extrapolation failed";
0046       return tsos;
0047     } else {
0048       T0 = tsos.globalPosition();
0049       T1 = tsos.globalMomentum().unit();
0050       GlobalVector D = L.distance(T0);
0051       distance = D.mag();
0052       if (fabs(old_distance - distance) < 0.000001) {
0053         refining = false;
0054       }
0055       if (old_distance - distance < 0.) {
0056         refining = false;
0057         LogDebug("TrajectoryExtrapolatorToLine") << "TETL- stop to skip loops";
0058       }
0059     }
0060   }
0061   //
0062   // Now propagate with errors and (not for the moment) perform rotation
0063   //
0064   // Origin of local system: point of closest approach on the line
0065   // (w.r.t. to tangent to helix at last iteration)
0066   //
0067   Line T(T0, T1);
0068   GlobalPoint origin(L.closerPointToLine(T));
0069   //
0070   // Axes of local system:
0071   //   x from line to helix at closest approach
0072   //   z along the helix
0073   //   y to complete right-handed system
0074   //
0075   GlobalVector ZZ(T1.unit());
0076   GlobalVector YY(ZZ.cross(T0 - origin).unit());
0077   GlobalVector XX(YY.cross(ZZ));
0078   Surface::RotationType rot(XX, YY, ZZ);
0079   ReferenceCountingPointer<Plane> surface = Plane::build(origin, rot);
0080   TrajectoryStateOnSurface tsos = p->propagate(fts, *surface);
0081
0082   return tsos;
0083 }