NuclearSeedGenerator/src/SeedFromNuclearInteraction.cc

0001 #include <memory>
0002
0003 #include "RecoTracker/NuclearSeedGenerator/interface/SeedFromNuclearInteraction.h"
0004
0005 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0006 #include "MagneticField/Engine/interface/MagneticField.h"
0007
0008 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0009 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0010 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
0011
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013
0014 SeedFromNuclearInteraction::SeedFromNuclearInteraction(const Propagator* prop,
0015                                                        const TrackerGeometry* geom,
0016                                                        double iPtMin)
0017     : ptMin(iPtMin), thePropagator(prop), theTrackerGeom(geom) {
0018   isValid_ = true;
0019   initialTSOS_ = std::make_shared<TrajectoryStateOnSurface>();
0020   updatedTSOS_ = std::make_shared<TrajectoryStateOnSurface>();
0021   freeTS_ = std::make_shared<FreeTrajectoryState>();
0022 }
0023
0024 //----------------------------------------------------------------------
0025 void SeedFromNuclearInteraction::setMeasurements(const TSOS& inner_TSOS,
0026                                                  ConstRecHitPointer ihit,
0027                                                  ConstRecHitPointer ohit) {
0028   // delete pointer to TrackingRecHits
0029   theHits.clear();
0030
0031   // get the inner and outer transient TrackingRecHits
0032   innerHit_ = ihit;
0033   outerHit_ = ohit;
0034
0035   //theHits.push_back(  inner_TM.recHit() ); // put temporarily - TODO: remove this line
0036   theHits.push_back(outerHit_);
0037
0038   initialTSOS_ = std::make_shared<TrajectoryStateOnSurface>(inner_TSOS);
0039
0040   // calculate the initial FreeTrajectoryState.
0041   freeTS_.reset(stateWithError());
0042
0043   // check transverse momentum
0044   if (freeTS_->momentum().perp() < ptMin) {
0045     isValid_ = false;
0046   } else {
0047     // convert freeTS_ into a persistent TSOS on the outer surface
0048     isValid_ = construct();
0049   }
0050 }
0051 //----------------------------------------------------------------------
0052 void SeedFromNuclearInteraction::setMeasurements(TangentHelix& thePrimaryHelix,
0053                                                  const TSOS& inner_TSOS,
0054                                                  ConstRecHitPointer ihit,
0055                                                  ConstRecHitPointer ohit) {
0056   // delete pointer to TrackingRecHits
0057   theHits.clear();
0058
0059   // get the inner and outer transient TrackingRecHits
0060   innerHit_ = ihit;
0061   outerHit_ = ohit;
0062
0063   GlobalPoint innerPos =
0064       theTrackerGeom->idToDet(innerHit_->geographicalId())->surface().toGlobal(innerHit_->localPosition());
0065   GlobalPoint outerPos =
0066       theTrackerGeom->idToDet(outerHit_->geographicalId())->surface().toGlobal(outerHit_->localPosition());
0067
0068   TangentHelix helix(thePrimaryHelix, outerPos, innerPos);
0069
0070   theHits.push_back(innerHit_);
0071   theHits.push_back(outerHit_);
0072
0073   initialTSOS_ = std::make_shared<TrajectoryStateOnSurface>(inner_TSOS);
0074
0075   // calculate the initial FreeTrajectoryState from the inner and outer TM assuming that the helix equation is already known.
0076   freeTS_.reset(stateWithError(helix));
0077
0078   if (freeTS_->momentum().perp() < ptMin) {
0079     isValid_ = false;
0080   } else {
0081     // convert freeTS_ into a persistent TSOS on the outer surface
0082     isValid_ = construct();
0083   }
0084 }
0085 //----------------------------------------------------------------------
0086 FreeTrajectoryState* SeedFromNuclearInteraction::stateWithError() const {
0087   // Calculation of the helix assuming that the secondary track has the same direction
0088   // than the primary track and pass through the inner and outer hits.
0089   GlobalVector direction = initialTSOS_->globalDirection();
0090   GlobalPoint inner = initialTSOS_->globalPosition();
0091   TangentHelix helix(direction, inner, outerHitPosition());
0092
0093   return stateWithError(helix);
0094 }
0095 //----------------------------------------------------------------------
0096 FreeTrajectoryState* SeedFromNuclearInteraction::stateWithError(TangentHelix& helix) const {
0097   //   typedef TkRotation<float> Rotation;
0098
0099   GlobalVector dirAtVtx = helix.directionAtVertex();
0100   const MagneticField& mag = initialTSOS_->globalParameters().magneticField();
0101
0102   // Get the global parameters of the trajectory
0103   // we assume that the magnetic field at the vertex is equal to the magnetic field at the inner TM.
0104   GlobalTrajectoryParameters gtp(
0105       helix.vertexPoint(), dirAtVtx, helix.charge(mag.inTesla(helix.vertexPoint()).z()) / helix.rho(), 0, &mag);
0106
0107   // Error matrix in a frame where z is in the direction of the track at the vertex
0108   AlgebraicSymMatrix66 primaryError(initialTSOS_->cartesianError().matrix());
0109   double p_max = initialTSOS_->globalParameters().momentum().mag();
0110   AlgebraicMatrix33 rot = this->rotationMatrix(dirAtVtx);
0111
0112   AlgebraicMatrix66 globalRotation;
0113   globalRotation.Place_at(rot, 0, 0);
0114   globalRotation.Place_at(rot, 3, 3);
0115   AlgebraicSymMatrix66 primaryErrorInNewFrame = ROOT::Math::Similarity(globalRotation, primaryError);
0116
0117   AlgebraicSymMatrix66 secondaryErrorInNewFrame = AlgebraicMatrixID();
0118   double p_perp_max = 2;  // energy max of a secondary track emited perpendicularly to the
0119                           // primary track is +/- 2 GeV
0120   secondaryErrorInNewFrame(0, 0) = primaryErrorInNewFrame(0, 0) + helix.vertexError() * p_perp_max / p_max;
0121   secondaryErrorInNewFrame(1, 1) = primaryErrorInNewFrame(1, 1) + helix.vertexError() * p_perp_max / p_max;
0122   secondaryErrorInNewFrame(2, 2) = helix.vertexError() * helix.vertexError();
0123   secondaryErrorInNewFrame(3, 3) = p_perp_max * p_perp_max;
0124   secondaryErrorInNewFrame(4, 4) = p_perp_max * p_perp_max;
0125   secondaryErrorInNewFrame(5, 5) = p_max * p_max;
0126
0127   AlgebraicSymMatrix66 secondaryError = ROOT::Math::SimilarityT(globalRotation, secondaryErrorInNewFrame);
0128
0129   return new FreeTrajectoryState(gtp, CartesianTrajectoryError(secondaryError));
0130 }
0131
0132 //----------------------------------------------------------------------
0133 bool SeedFromNuclearInteraction::construct() {
0134   // loop on all hits in theHits
0135   KFUpdator theUpdator;
0136
0137   const TrackingRecHit* hit = nullptr;
0138
0139   LogDebug("NuclearSeedGenerator") << "Seed ** initial state " << freeTS_->cartesianError().matrix();
0140
0141   for (unsigned int iHit = 0; iHit < theHits.size(); iHit++) {
0142     hit = theHits[iHit]->hit();
0143     TrajectoryStateOnSurface state =
0144         (iHit == 0)
0145             ? thePropagator->propagate(*freeTS_, theTrackerGeom->idToDet(hit->geographicalId())->surface())
0146             : thePropagator->propagate(*updatedTSOS_, theTrackerGeom->idToDet(hit->geographicalId())->surface());
0147
0148     if (!state.isValid())
0149       return false;
0150
0151     const TransientTrackingRecHit::ConstRecHitPointer& tth = theHits[iHit];
0152     updatedTSOS_ = std::make_shared<TrajectoryStateOnSurface>(theUpdator.update(state, *tth));
0153   }
0154
0155   LogDebug("NuclearSeedGenerator") << "Seed ** updated state " << updatedTSOS_->cartesianError().matrix();
0156
0157   pTraj = trajectoryStateTransform::persistentState(*updatedTSOS_, outerHitDetId().rawId());
0158   return true;
0159 }
0160
0161 //----------------------------------------------------------------------
0162 edm::OwnVector<TrackingRecHit> SeedFromNuclearInteraction::hits() const {
0163   recHitContainer _hits;
0164   for (ConstRecHitContainer::const_iterator it = theHits.begin(); it != theHits.end(); it++) {
0165     _hits.push_back(it->get()->hit()->clone());
0166   }
0167   return _hits;
0168 }
0169 //----------------------------------------------------------------------
0170 AlgebraicMatrix33 SeedFromNuclearInteraction::rotationMatrix(const GlobalVector& perp) const {
0171   AlgebraicMatrix33 result;
0172
0173   // z axis coincides with perp
0174   GlobalVector zAxis = perp.unit();
0175
0176   // x axis has no global Z component
0177   GlobalVector xAxis;
0178   if (zAxis.x() != 0 || zAxis.y() != 0) {
0179     // precision is not an issue here, just protect against divizion by zero
0180     xAxis = GlobalVector(-zAxis.y(), zAxis.x(), 0).unit();
0181   } else {  // perp coincides with global Z
0182     xAxis = GlobalVector(1, 0, 0);
0183   }
0184
0185   // y axis obtained by cross product
0186   GlobalVector yAxis(zAxis.cross(xAxis));
0187
0188   result(0, 0) = xAxis.x();
0189   result(0, 1) = xAxis.y();
0190   result(0, 2) = xAxis.z();
0191   result(1, 0) = yAxis.x();
0192   result(1, 1) = yAxis.y();
0193   result(1, 2) = yAxis.z();
0194   result(2, 0) = zAxis.x();
0195   result(2, 1) = zAxis.y();
0196   result(2, 2) = zAxis.z();
0197   return result;
0198 }