PatternTools/src/TSCPBuilderNoMaterial.cc

0001 #include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
0002 #include "DataFormats/GeometrySurface/interface/Surface.h"
0003 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0004 #include "TrackingTools/AnalyticalJacobians/interface/AnalyticalCurvilinearJacobian.h"
0005 #include "DataFormats/GeometrySurface/interface/BoundPlane.h"
0006 #include "TrackingTools/GeomPropagators/interface/HelixBarrelPlaneCrossingByCircle.h"
0007 #include "TrackingTools/TrajectoryParametrization/interface/TrajectoryStateExceptions.h"
0008 #include "MagneticField/Engine/interface/MagneticField.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010
0011 TrajectoryStateClosestToPoint TSCPBuilderNoMaterial::operator()(const FTS& originalFTS,
0012                                                                 const GlobalPoint& referencePoint) const {
0013   if (positionEqual(referencePoint, originalFTS.position()))
0014     return constructTSCP(originalFTS, referencePoint);
0015
0016   // Now do the propagation or whatever...
0017
0018   PairBoolFTS newStatePair = createFTSatTransverseImpactPoint(originalFTS, referencePoint);
0019   if (newStatePair.first) {
0020     return constructTSCP(newStatePair.second, referencePoint);
0021   } else {
0022     return TrajectoryStateClosestToPoint();
0023   }
0024 }
0025
0026 TrajectoryStateClosestToPoint TSCPBuilderNoMaterial::operator()(const TSOS& originalTSOS,
0027                                                                 const GlobalPoint& referencePoint) const {
0028   if (positionEqual(referencePoint, originalTSOS.globalPosition()))
0029     return constructTSCP(*originalTSOS.freeState(), referencePoint);
0030
0031   // Now do the propagation
0032
0033   PairBoolFTS newStatePair = createFTSatTransverseImpactPoint(*originalTSOS.freeState(), referencePoint);
0034   if (newStatePair.first) {
0035     return constructTSCP(newStatePair.second, referencePoint);
0036   } else {
0037     return TrajectoryStateClosestToPoint();
0038   }
0039 }
0040
0041 TSCPBuilderNoMaterial::PairBoolFTS TSCPBuilderNoMaterial::createFTSatTransverseImpactPoint(
0042     const FTS& originalFTS, const GlobalPoint& referencePoint) const {
0043   //
0044   // Straight line approximation? |rho|<1.e-10 equivalent to ~ 1um
0045   // difference in transversal position at 10m.
0046   //
0047   if (fabs(originalFTS.transverseCurvature()) < 1.e-10) {
0048     return createFTSatTransverseImpactPointNeutral(originalFTS, referencePoint);
0049   } else {
0050     return createFTSatTransverseImpactPointCharged(originalFTS, referencePoint);
0051   }
0052 }
0053
0054 TSCPBuilderNoMaterial::PairBoolFTS TSCPBuilderNoMaterial::createFTSatTransverseImpactPointCharged(
0055     const FTS& originalFTS, const GlobalPoint& referencePoint) const {
0056   GlobalVector pvecOrig = originalFTS.momentum();
0057   GlobalPoint xvecOrig = originalFTS.position();
0058   double kappa = originalFTS.transverseCurvature();
0059   double pxOrig = pvecOrig.x();
0060   double pyOrig = pvecOrig.y();
0061   double pzOrig = pvecOrig.z();
0062   double xOrig = xvecOrig.x();
0063   double yOrig = xvecOrig.y();
0064   double zOrig = xvecOrig.z();
0065
0066   //  double fac = 1./originalFTS.charge()/MagneticField::inInverseGeV(referencePoint).z();
0067   double fac = 1. / originalFTS.charge() / (originalFTS.parameters().magneticField().inInverseGeV(referencePoint).z());
0068   GlobalVectorDouble xOrig2Centre = GlobalVectorDouble(fac * pyOrig, -fac * pxOrig, 0.);
0069   GlobalVectorDouble xOrigProj = GlobalVectorDouble(xOrig, yOrig, 0.);
0070   GlobalVectorDouble xRefProj = GlobalVectorDouble(referencePoint.x(), referencePoint.y(), 0.);
0071   GlobalVectorDouble deltax = xRefProj - xOrigProj - xOrig2Centre;
0072   GlobalVectorDouble ndeltax = deltax.unit();
0073
0074   PropagationDirection direction = anyDirection;
0075   Surface::PositionType pos(referencePoint);
0076   // Need to define plane with orientation as the
0077   // ImpactPointSurface
0078   GlobalVector X(ndeltax.x(), ndeltax.y(), ndeltax.z());
0079   GlobalVector Y(0., 0., 1.);
0080   Surface::RotationType rot(X, Y);
0081   BoundPlane* plane = new BoundPlane(pos, rot);
0082   // Using Teddy's HelixBarrelPlaneCrossingByCircle for general barrel planes.
0083   // A large variety of other,
0084   // direct solutions turned out to be not so stable numerically.
0085   HelixBarrelPlaneCrossingByCircle planeCrossing(HelixPlaneCrossing::PositionType(xOrig, yOrig, zOrig),
0086                                                  HelixPlaneCrossing::DirectionType(pxOrig, pyOrig, pzOrig),
0087                                                  kappa,
0088                                                  direction);
0089   std::pair<bool, double> propResult = planeCrossing.pathLength(*plane);
0090   if (!propResult.first) {
0091     edm::LogWarning("TSCPBuilderNoMaterial") << "Propagation to perigee plane failed!";
0092     return PairBoolFTS(false, FreeTrajectoryState());
0093   }
0094   double s = propResult.second;
0095   HelixPlaneCrossing::PositionType xGen = planeCrossing.position(s);
0096   GlobalPoint xPerigee = GlobalPoint(xGen.x(), xGen.y(), xGen.z());
0097   // direction (reconverted to GlobalVector, renormalised)
0098   HelixPlaneCrossing::DirectionType pGen = planeCrossing.direction(s);
0099   pGen *= pvecOrig.mag() / pGen.mag();
0100   GlobalVector pPerigee = GlobalVector(pGen.x(), pGen.y(), pGen.z());
0101   delete plane;
0102
0103   if (originalFTS.hasError()) {
0104     const AlgebraicSymMatrix55& errorMatrix = originalFTS.curvilinearError().matrix();
0105     AnalyticalCurvilinearJacobian curvilinJacobian(originalFTS.parameters(), xPerigee, pPerigee, s);
0106     const AlgebraicMatrix55& jacobian = curvilinJacobian.jacobian();
0107     CurvilinearTrajectoryError cte(ROOT::Math::Similarity(jacobian, errorMatrix));
0108
0109     return PairBoolFTS(
0110         true,
0111         FreeTrajectoryState(GlobalTrajectoryParameters(
0112                                 xPerigee, pPerigee, originalFTS.charge(), &(originalFTS.parameters().magneticField())),
0113                             cte));
0114   } else {
0115     return PairBoolFTS(true,
0116                        FreeTrajectoryState(GlobalTrajectoryParameters(
0117                            xPerigee, pPerigee, originalFTS.charge(), &(originalFTS.parameters().magneticField()))));
0118   }
0119 }
0120
0121 TSCPBuilderNoMaterial::PairBoolFTS TSCPBuilderNoMaterial::createFTSatTransverseImpactPointNeutral(
0122     const FTS& originalFTS, const GlobalPoint& referencePoint) const {
0123   GlobalPoint xvecOrig = originalFTS.position();
0124   double phi = originalFTS.momentum().phi();
0125   double theta = originalFTS.momentum().theta();
0126   double xOrig = xvecOrig.x();
0127   double yOrig = xvecOrig.y();
0128   double zOrig = xvecOrig.z();
0129   double xR = referencePoint.x();
0130   double yR = referencePoint.y();
0131
0132   double s2D = (xR - xOrig) * cos(phi) + (yR - yOrig) * sin(phi);
0133   double s = s2D / sin(theta);
0134   double xGen = xOrig + s2D * cos(phi);
0135   double yGen = yOrig + s2D * sin(phi);
0136   double zGen = zOrig + s2D / tan(theta);
0137   GlobalPoint xPerigee = GlobalPoint(xGen, yGen, zGen);
0138
0139   GlobalVector pPerigee = originalFTS.momentum();
0140
0141   if (originalFTS.hasError()) {
0142     const AlgebraicSymMatrix55& errorMatrix = originalFTS.curvilinearError().matrix();
0143     AnalyticalCurvilinearJacobian curvilinJacobian(originalFTS.parameters(), xPerigee, pPerigee, s);
0144     const AlgebraicMatrix55& jacobian = curvilinJacobian.jacobian();
0145     CurvilinearTrajectoryError cte(ROOT::Math::Similarity(jacobian, errorMatrix));
0146
0147     return PairBoolFTS(
0148         true,
0149         FreeTrajectoryState(GlobalTrajectoryParameters(
0150                                 xPerigee, pPerigee, originalFTS.charge(), &(originalFTS.parameters().magneticField())),
0151                             cte));
0152   } else {
0153     return PairBoolFTS(true,
0154                        FreeTrajectoryState(GlobalTrajectoryParameters(
0155                            xPerigee, pPerigee, originalFTS.charge(), &(originalFTS.parameters().magneticField()))));
0156   }
0157 }