AnalyticalPropagator.cc

Go to the documentation of this file.

#include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"

#include "DataFormats/GeometrySurface/interface/Cylinder.h"
#include "DataFormats/GeometrySurface/interface/TangentPlane.h"
#include "DataFormats/GeometrySurface/interface/Plane.h"

#include "MagneticField/Engine/interface/MagneticField.h"

#include "TrackingTools/GeomPropagators/interface/PropagationExceptions.h"
#include "TrackingTools/GeomPropagators/interface/StraightLinePlaneCrossing.h"
#include "TrackingTools/GeomPropagators/interface/StraightLineBarrelCylinderCrossing.h"
#include "TrackingTools/GeomPropagators/interface/OptimalHelixPlaneCrossing.h"
#include "TrackingTools/GeomPropagators/interface/HelixBarrelCylinderCrossing.h"
#include "TrackingTools/AnalyticalJacobians/interface/AnalyticalCurvilinearJacobian.h"
#include "TrackingTools/GeomPropagators/interface/PropagationDirectionFromPath.h"
#include "TrackingTools/TrajectoryState/interface/SurfaceSideDefinition.h"
#include "TrackingTools/GeomPropagators/interface/PropagationExceptions.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/Likely.h"

#include <cmath>

using namespace SurfaceSideDefinition;

std::pair<TrajectoryStateOnSurface,double>
AnalyticalPropagator::propagateWithPath(const FreeTrajectoryState& fts, 
                    const Plane& plane) const
{
  // check curvature
  float rho = fts.transverseCurvature();
  
  // propagate parameters
  GlobalPoint x;
  GlobalVector p;
  double s;
  
  // check if already on plane
  if LIKELY (plane.localZclamped(fts.position()) !=0)  {
      // propagate
      bool parametersOK = this->propagateParametersOnPlane(fts, plane, x, p, s);
      // check status and deltaPhi limit
      float dphi2 = float(s)*rho;
      dphi2 = dphi2*dphi2*fts.momentum().perp2();
      if UNLIKELY( !parametersOK || dphi2>theMaxDPhi2*fts.momentum().mag2() )  return TsosWP(TrajectoryStateOnSurface(),0.);
    }
  else {
    LogDebug("AnalyticalPropagator")<<"not going anywhere. Already on surface.\n"
                    <<"plane.localZ(fts.position()): "<<plane.localZ(fts.position())<<"\n"
                    <<"plane.position().mag(): "<<plane.position().mag() <<"\n"
                    <<"plane.posPrec: "<<plane.posPrec();
    x = fts.position();
    p = fts.momentum();
    s = 0;
  }
  //
  // Compute propagated state and check change in curvature
  //
  GlobalTrajectoryParameters gtp(x,p,fts.charge(),theField);
  if UNLIKELY(std::abs(gtp.transverseCurvature()-rho)>theMaxDBzRatio*std::abs(rho) ) 
    return TsosWP(TrajectoryStateOnSurface(),0.);
  //
  // construct TrajectoryStateOnSurface
  //
  return propagatedStateWithPath(fts,plane,gtp,s);
}


std::pair<TrajectoryStateOnSurface,double>
AnalyticalPropagator::propagateWithPath(const FreeTrajectoryState& fts, 
                    const Cylinder& cylinder) const
{
  // check curvature
  auto rho = fts.transverseCurvature();

  // propagate parameters
  GlobalPoint x;
  GlobalVector p;
  double s = 0;

  bool parametersOK = this->propagateParametersOnCylinder(fts, cylinder, x, p, s);
  // check status and deltaPhi limit
  float dphi2 = s*rho;
  dphi2 = dphi2*dphi2*fts.momentum().perp2();
  if UNLIKELY( !parametersOK || dphi2>theMaxDPhi2*fts.momentum().mag2() )  return TsosWP(TrajectoryStateOnSurface(),0.);
  //
  // Compute propagated state and check change in curvature
  //
  GlobalTrajectoryParameters gtp(x,p,fts.charge(),theField);
  if UNLIKELY( std::abs(gtp.transverseCurvature()-rho)>theMaxDBzRatio*std::abs(rho) ) 
    return TsosWP(TrajectoryStateOnSurface(),0.);
  //
  // create result TSOS on TangentPlane (local parameters & errors are better defined)
  //

  //try {
    ConstReferenceCountingPointer<TangentPlane> plane(cylinder.tangentPlane(x));  // need to be here until tsos is created!
    return propagatedStateWithPath(fts,*plane,gtp,s);
  /*
  } catch(...) {
    std::cout << "wrong tangent to cylinder " << x 
              << " pos, rad " << cylinder.position() << " " << cylinder.radius()
              << std::endl;
    return TsosWP(TrajectoryStateOnSurface(),0.);
  }
  */
}

std::pair<TrajectoryStateOnSurface,double>
AnalyticalPropagator::propagatedStateWithPath (const FreeTrajectoryState& fts, 
                           const Surface& surface, 
                           const GlobalTrajectoryParameters& gtp, 
                           const double& s) const
{
  //
  // for forward propagation: state is before surface,
  // for backward propagation: state is after surface
  //
  SurfaceSide side = PropagationDirectionFromPath()(s,propagationDirection())==alongMomentum 
    ? beforeSurface : afterSurface;
  // 
  //
  // error propagation (if needed) and conversion to a TrajectoryStateOnSurface
  //
  if (fts.hasError()) {
    //
    // compute jacobian
    //
    AnalyticalCurvilinearJacobian analyticalJacobian(fts.parameters(), gtp.position(), gtp.momentum(), s);
    const AlgebraicMatrix55 &jacobian = analyticalJacobian.jacobian();
    // CurvilinearTrajectoryError cte(ROOT::Math::Similarity(jacobian, fts.curvilinearError().matrix()));
    return TsosWP(TrajectoryStateOnSurface(gtp,
                       ROOT::Math::Similarity(jacobian, fts.curvilinearError().matrix()),
                       surface,side),s);
  }
  else {
    //
    // return state without errors
    //
    return TsosWP(TrajectoryStateOnSurface(gtp,surface,side),s);
  }
}

bool AnalyticalPropagator::propagateParametersOnCylinder(
  const FreeTrajectoryState& fts, const Cylinder& cylinder, 
  GlobalPoint& x, GlobalVector& p, double& s) const
{

  GlobalPoint const & sp = cylinder.position();
  if UNLIKELY(sp.x()!=0. || sp.y()!=0.) {
    throw PropagationException("Cannot propagate to an arbitrary cylinder");
  }
  // preset output
  x = fts.position();
  p = fts.momentum();
  s = 0;
  // (transverse) curvature
  auto rho = fts.transverseCurvature();
  //
  // Straight line approximation? |rho|<1.e-10 equivalent to ~ 1um 
  // difference in transversal position at 10m.
  //
  if UNLIKELY( std::abs(rho)<1.e-10f )
    return propagateWithLineCrossing(fts.position(),p,cylinder,x,s);
  //
  // Helix case
  //
  // check for possible intersection
  constexpr float tolerance = 1.e-4; // 1 micron distance
  auto rdiff = x.perp() - cylinder.radius();
  if ( std::abs(rdiff) < tolerance )  return true;
  //
  // Instantiate HelixBarrelCylinderCrossing and get solutions
  //
  HelixBarrelCylinderCrossing cylinderCrossing(fts.position(),fts.momentum(),rho,
                           propagationDirection(),cylinder);
  if UNLIKELY( !cylinderCrossing.hasSolution() )  return false;
  // path length
  s = cylinderCrossing.pathLength();
  // point
  x = cylinderCrossing.position();
  // direction (renormalised)
  p = cylinderCrossing.direction().unit()*fts.momentum().mag();
  return true;
}
  
bool 
AnalyticalPropagator::propagateParametersOnPlane(const FreeTrajectoryState& fts, 
                         const Plane& plane, 
                         GlobalPoint& x, 
                         GlobalVector& p, 
                         double& s) const
{
  // initialisation of position, momentum and path length
  x = fts.position();
  p = fts.momentum();
  s = 0;
  // (transverse) curvature
  auto rho = fts.transverseCurvature();
  //
  // Straight line approximation? |rho|<1.e-10 equivalent to ~ 1um 
  // difference in transversal position at 10m.
  //
  if UNLIKELY( std::abs(rho)<1.e-10f )
    return propagateWithLineCrossing(fts.position(),p,plane,x,s);
  //
  // Helix case 
  //

  //
  // Frame-independant point and vector are created explicitely to 
  // avoid confusing gcc (refuses to compile with temporary objects
  // in the constructor).
  //
  HelixPlaneCrossing::PositionType helixPos(x);
  HelixPlaneCrossing::DirectionType helixDir(p);
  if LIKELY(isOldPropagationType) {
      OptimalHelixPlaneCrossing planeCrossing(plane,helixPos,helixDir,rho,propagationDirection());
      return propagateWithHelixCrossing(*planeCrossing,plane,fts.momentum().mag(),x,p,s);
    }




  //--- Alternative implementation to be used for the propagation of the parameters  of looping
  //    particles that cross twice the (infinite) surface of the plane. It is not trivial to determine
  //    which of the two intersections has to be returned.
  
  //---- FIXME: WHAT FOLLOWS HAS TO BE REWRITTEN IN A CLEANER (AND CPU-OPTIMIZED) WAY ---------
  LogDebug("AnalyticalPropagator") << "In AnaliticalProp, calling HAPC " << "\n"
                   << "plane is centered in xyz: " 
                   << plane.position().x() << " , "
                   << plane.position().y() << " , " 
                   << plane.position().z() << "\n";
  
  
  GlobalPoint gp1 = fts.position();
  GlobalVector gm1 = fts.momentum();
  double s1 = 0;
  double rho1 = fts.transverseCurvature();
  HelixPlaneCrossing::PositionType helixPos1(gp1);
  HelixPlaneCrossing::DirectionType helixDir1(gm1);
  LogDebug("AnalyticalPropagator") << "gp1 before calling planeCrossing1: " << gp1 << "\n";
  OptimalHelixPlaneCrossing planeCrossing1(plane,helixPos1,helixDir1,rho1,propagationDirection());
  
  HelixPlaneCrossing::PositionType xGen;
  HelixPlaneCrossing::DirectionType pGen;
      
  double tolerance(0.0050);
  if(propagationDirection()==oppositeToMomentum)
    tolerance *=-1;
  
  bool check1 = propagateWithHelixCrossing(*planeCrossing1,plane,fts.momentum().mag(),gp1,gm1,s1);
  double dphi1 = fabs(fts.momentum().phi()-gm1.phi());
  LogDebug("AnalyticalPropagator") << "check1, s1, dphi, gp1: " 
                   << check1 << " , "
                   << s1 << " , " 
                   << dphi1 << " , "
                   << gp1 << "\n";
  
  //move forward a bit to avoid that the propagator doesn't propagate because the state is already on surface.
  //we want to go to the other point of intersection between the helix and the plane
  xGen = (*planeCrossing1).position(s1+tolerance);
  pGen = (*planeCrossing1).direction(s1+tolerance);
  
  /*
    if(!check1 || s1>170 ){
    //PropagationDirection newDir = (propagationDirection() == alongMomentum) ? oppositeToMomentum : alongMomentum;
    PropagationDirection newDir = anyDirection;
    HelixArbitraryPlaneCrossing  planeCrossing1B(helixPos1,helixDir1,rho1,newDir);
    check1 = propagateWithHelixCrossing(planeCrossing1B,plane,fts.momentum().mag(),gp1,gm1,s1);
    LogDebug("AnalyticalPropagator") << "after second attempt, check1, s1,gp1: "
    << check1 << " , "
    << s1 << " , " << gp1 << "\n";
    
    xGen = planeCrossing1B.position(s1+tolerance);
    pGen = planeCrossing1B.direction(s1+tolerance);
    }
      */
  
  if(!check1){
    LogDebug("AnalyticalPropagator") << "failed also second attempt. No idea what to do, then bailout" << "\n";
  }
  
  
  pGen *= gm1.mag()/pGen.mag();
  GlobalPoint gp2(xGen);
  GlobalVector gm2(pGen);
  double s2 = 0;
  double rho2 = rho1;
  HelixPlaneCrossing::PositionType helixPos2(gp2);
  HelixPlaneCrossing::DirectionType helixDir2(gm2);
  OptimalHelixPlaneCrossing planeCrossing2(plane,helixPos2,helixDir2,rho2,propagationDirection());
  
  bool check2 = propagateWithHelixCrossing(*planeCrossing2,plane,gm2.mag(),gp2,gm2,s2);
  
  if(!check2){
    x = gp1;
    p = gm1;
    s = s1;
    return check1;
  }
  
  if(!check1){
    edm::LogError("AnalyticalPropagator") << "LOGIC ERROR: I should not have entered here!" << "\n";
    return false;
  }
  
  
  LogDebug("AnalyticalPropagator")  << "check2, s2, gp2: " 
                    << check2 << " , "
                    << s2 << " , " << gp2 << "\n";
  
  
  double dist1 = (plane.position()-gp1).perp();
  double dist2 = (plane.position()-gp2).perp();
  
  
  LogDebug("AnalyticalPropagator") << "propDir, dist1, dist2: " 
                   << propagationDirection() << " , "
                   << dist1 << " , " 
                   << dist2 << "\n";
  
  //If there are two solutions, the one which is the closest to the module's center is chosen
  if(dist1<2*dist2){
    x = gp1;
    p = gm1;
    s = s1;
    return check1;
  }else if(dist2<2*dist1){
    x = gp2;
    p = gm2;
    s = s1+s2+tolerance;
    return check2;
  }else{
    if(fabs(s1)<fabs(s2)){
      x = gp1;
      p = gm1;
      s = s1;
      return check1;
    }else{
      x = gp2;
      p = gm2;
      s = s1+s2+tolerance;
      return check2;
    }   
  }
  
  //-------- END of ugly piece of code  ---------------
  
}

bool
AnalyticalPropagator::propagateWithLineCrossing (const GlobalPoint& x0, 
                         const GlobalVector& p0,
                         const Plane& plane,
                         GlobalPoint& x, double& s) const {
  //
  // Instantiate auxiliary object for finding intersection.
  // Frame-independant point and vector are created explicitely to 
  // avoid confusing gcc (refuses to compile with temporary objects
  // in the constructor).
  //
  StraightLinePlaneCrossing::PositionType pos(x0);
  StraightLinePlaneCrossing::DirectionType dir(p0);
  StraightLinePlaneCrossing planeCrossing(pos,dir,propagationDirection());
  //
  // get solution
  //
  std::pair<bool,double> propResult = planeCrossing.pathLength(plane);
  if ( !propResult.first )  return false;
  s = propResult.second;
  // point (reconverted to GlobalPoint)
  x = GlobalPoint(planeCrossing.position(s));
  //
  return true;
}

bool
AnalyticalPropagator::propagateWithLineCrossing (const GlobalPoint& x0, 
                         const GlobalVector& p0,
                         const Cylinder& cylinder,
                         GlobalPoint& x, double& s) const {
  //
  // Instantiate auxiliary object for finding intersection.
  // Frame-independant point and vector are created explicitely to 
  // avoid confusing gcc (refuses to compile with temporary objects
  // in the constructor).
  //
  StraightLineBarrelCylinderCrossing cylCrossing(x0,p0,propagationDirection());
  //
  // get solution
  //
  std::pair<bool,double> propResult = cylCrossing.pathLength(cylinder);
  if ( !propResult.first )  return false;
  s = propResult.second;
  // point (reconverted to GlobalPoint)
  x = cylCrossing.position(s);
  //
  return true;
}

bool
AnalyticalPropagator::propagateWithHelixCrossing (HelixPlaneCrossing& planeCrossing,
                          const Plane& plane,
                          const float pmag,
                          GlobalPoint& x,
                          GlobalVector& p,
                          double& s) const {
  // get solution
  std::pair<bool,double> propResult = planeCrossing.pathLength(plane);
  if UNLIKELY( !propResult.first )  return false;

  s = propResult.second;
  x = GlobalPoint(planeCrossing.position(s));
  // direction (reconverted to GlobalVector, renormalised)
  GlobalVector pGen =  GlobalVector(planeCrossing.direction(s));
  pGen *= pmag/pGen.mag();
  p = pGen;
  //
  return true;
}