TSCPBuilderNoMaterial.cc

Go to the documentation of this file.

#include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
#include "DataFormats/GeometrySurface/interface/Surface.h" 
#include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h" 
#include "TrackingTools/AnalyticalJacobians/interface/AnalyticalCurvilinearJacobian.h"
#include "DataFormats/GeometrySurface/interface/BoundPlane.h"
#include "TrackingTools/GeomPropagators/interface/HelixBarrelPlaneCrossingByCircle.h"
#include "TrackingTools/TrajectoryParametrization/interface/TrajectoryStateExceptions.h"
#include "MagneticField/Engine/interface/MagneticField.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

TrajectoryStateClosestToPoint 
TSCPBuilderNoMaterial::operator() (const FTS& originalFTS, 
    const GlobalPoint& referencePoint) const
{
  if (positionEqual(referencePoint, originalFTS.position()))
    return constructTSCP(originalFTS, referencePoint);

  // Now do the propagation or whatever...

  PairBoolFTS newStatePair =
    createFTSatTransverseImpactPoint(originalFTS, referencePoint);
  if (newStatePair.first) {
    return constructTSCP(newStatePair.second, referencePoint);
  } else {
    return TrajectoryStateClosestToPoint();
  }  
}

TrajectoryStateClosestToPoint 
TSCPBuilderNoMaterial::operator() (const TSOS& originalTSOS, 
    const GlobalPoint& referencePoint) const
{
  if (positionEqual(referencePoint, originalTSOS.globalPosition()))
    return constructTSCP(*originalTSOS.freeState(), referencePoint);

  // Now do the propagation
  
  PairBoolFTS newStatePair =
    createFTSatTransverseImpactPoint(*originalTSOS.freeState(), referencePoint);
  if (newStatePair.first) {
    return constructTSCP(newStatePair.second, referencePoint);
  } else {
    return TrajectoryStateClosestToPoint();
  }
}

TSCPBuilderNoMaterial::PairBoolFTS 
TSCPBuilderNoMaterial::createFTSatTransverseImpactPoint(
    const FTS& originalFTS, const GlobalPoint& referencePoint) const 
{
  //
  // Straight line approximation? |rho|<1.e-10 equivalent to ~ 1um 
  // difference in transversal position at 10m.
  //
  if( fabs(originalFTS.transverseCurvature())<1.e-10 ) {
    return createFTSatTransverseImpactPointNeutral(originalFTS, referencePoint);
  } else {
    return createFTSatTransverseImpactPointCharged(originalFTS, referencePoint);
  }
}

TSCPBuilderNoMaterial::PairBoolFTS 
TSCPBuilderNoMaterial::createFTSatTransverseImpactPointCharged(
    const FTS& originalFTS, const GlobalPoint& referencePoint) const 
{

  GlobalVector pvecOrig = originalFTS.momentum();
  GlobalPoint xvecOrig = originalFTS.position();
  double kappa = originalFTS.transverseCurvature();
  double pxOrig = pvecOrig.x();
  double pyOrig = pvecOrig.y();
  double pzOrig = pvecOrig.z();
  double xOrig = xvecOrig.x();
  double yOrig = xvecOrig.y();
  double zOrig = xvecOrig.z();

//  double fac = 1./originalFTS.charge()/MagneticField::inInverseGeV(referencePoint).z();
  double fac = 1./originalFTS.charge()/
    (originalFTS.parameters().magneticField().inInverseGeV(referencePoint).z());
  GlobalVectorDouble xOrig2Centre = GlobalVectorDouble(fac * pyOrig, -fac * pxOrig, 0.);
  GlobalVectorDouble xOrigProj = GlobalVectorDouble(xOrig, yOrig, 0.);
  GlobalVectorDouble xRefProj = GlobalVectorDouble(referencePoint.x(), referencePoint.y(), 0.);
  GlobalVectorDouble deltax = xRefProj-xOrigProj-xOrig2Centre;
  GlobalVectorDouble ndeltax = deltax.unit();

  PropagationDirection direction = anyDirection;
  Surface::PositionType pos(referencePoint);
  // Need to define plane with orientation as the
  // ImpactPointSurface
  GlobalVector X(ndeltax.x(), ndeltax.y(), ndeltax.z());
  GlobalVector Y(0.,0.,1.);
  Surface::RotationType rot(X,Y);
  BoundPlane* plane = new BoundPlane(pos,rot);
  // Using Teddy's HelixBarrelPlaneCrossingByCircle for general barrel planes. 
  // A large variety of other,
  // direct solutions turned out to be not so stable numerically.
  HelixBarrelPlaneCrossingByCircle 
    planeCrossing(HelixPlaneCrossing::PositionType(xOrig, yOrig, zOrig),
          HelixPlaneCrossing::DirectionType(pxOrig, pyOrig, pzOrig), 
          kappa, direction);
  std::pair<bool,double> propResult = planeCrossing.pathLength(*plane);
  if ( !propResult.first ) {
    edm::LogWarning ("TSCPBuilderNoMaterial") << "Propagation to perigee plane failed!";
    return PairBoolFTS(false, FreeTrajectoryState() );
  }
  double s = propResult.second;
  HelixPlaneCrossing::PositionType xGen = planeCrossing.position(s);
  GlobalPoint xPerigee = GlobalPoint(xGen.x(),xGen.y(),xGen.z());
  // direction (reconverted to GlobalVector, renormalised)
  HelixPlaneCrossing::DirectionType pGen = planeCrossing.direction(s);
  pGen *= pvecOrig.mag()/pGen.mag();
  GlobalVector pPerigee = GlobalVector(pGen.x(),pGen.y(),pGen.z());
  delete plane;
          
  if (originalFTS.hasError()) {
    const AlgebraicSymMatrix55 &errorMatrix = originalFTS.curvilinearError().matrix();
    AnalyticalCurvilinearJacobian curvilinJacobian(originalFTS.parameters(), xPerigee,
                           pPerigee, s);
    const AlgebraicMatrix55 &jacobian = curvilinJacobian.jacobian();
    CurvilinearTrajectoryError cte( ROOT::Math::Similarity(jacobian, errorMatrix) );
  
    return PairBoolFTS(true,
    FreeTrajectoryState(GlobalTrajectoryParameters(xPerigee, pPerigee, originalFTS.charge(), 
                &(originalFTS.parameters().magneticField())), cte) );
  } 
  else {
    return PairBoolFTS(true,
    FreeTrajectoryState(GlobalTrajectoryParameters(xPerigee, pPerigee, originalFTS.charge(), 
                &(originalFTS.parameters().magneticField()))) );
  }
  
}


TSCPBuilderNoMaterial::PairBoolFTS 
TSCPBuilderNoMaterial::createFTSatTransverseImpactPointNeutral(const FTS& originalFTS, 
    const GlobalPoint& referencePoint) const 
{

  GlobalPoint xvecOrig = originalFTS.position();
  double phi = originalFTS.momentum().phi();
  double theta = originalFTS.momentum().theta();
  double xOrig = xvecOrig.x();
  double yOrig = xvecOrig.y();
  double zOrig = xvecOrig.z();
  double xR = referencePoint.x();
  double yR = referencePoint.y();

  double s2D = (xR - xOrig)  * cos(phi) + (yR - yOrig)  * sin(phi);
  double s = s2D / sin(theta);
  double xGen = xOrig + s2D*cos(phi);
  double yGen = yOrig + s2D*sin(phi);
  double zGen = zOrig + s2D/tan(theta);
  GlobalPoint xPerigee = GlobalPoint(xGen, yGen, zGen);

  GlobalVector pPerigee = originalFTS.momentum();
          
  if (originalFTS.hasError()) {
    const AlgebraicSymMatrix55 &errorMatrix = originalFTS.curvilinearError().matrix();
    AnalyticalCurvilinearJacobian curvilinJacobian(originalFTS.parameters(), xPerigee,
                           pPerigee, s);
    const AlgebraicMatrix55 &jacobian = curvilinJacobian.jacobian();
    CurvilinearTrajectoryError cte( ROOT::Math::Similarity(jacobian, errorMatrix) );
  
    return PairBoolFTS(true,
    FreeTrajectoryState(GlobalTrajectoryParameters(xPerigee, pPerigee, originalFTS.charge(), 
                &(originalFTS.parameters().magneticField())), cte));
  } 
  else {
    return PairBoolFTS(true,
    FreeTrajectoryState(GlobalTrajectoryParameters(xPerigee, pPerigee, originalFTS.charge(),
                &(originalFTS.parameters().magneticField()))) );
  }
  
}