SpecialSeedGenerators/src/CosmicTrackingRegion.cc

0001 #include "RecoTracker/SpecialSeedGenerators/interface/CosmicTrackingRegion.h"
0002 #include "TrackingTools/KalmanUpdators/interface/EtaPhiMeasurementEstimator.h"
0003 #include "RecoTracker/Record/interface/CkfComponentsRecord.h"
0004 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
0005
0006 #include "TrackingTools/MeasurementDet/interface/LayerMeasurements.h"
0007 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
0008 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0009 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0010 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0011 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
0012 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0013
0014 #include "DataFormats/GeometrySurface/interface/BoundPlane.h"
0015
0016 namespace {
0017   template <class T>
0018   T sqr(T t) {
0019     return t * t;
0020   }
0021 }  // namespace
0022
0023 using namespace std;
0024
0025 void CosmicTrackingRegion::checkTracks(reco::TrackCollection const& tracks, std::vector<bool>& mask) const {
0026   const math::XYZPoint regOrigin(origin().x(), origin().y(), origin().z());
0027
0028   assert(mask.size() == tracks.size());
0029   int i = -1;
0030   for (auto const& track : tracks) {
0031     i++;
0032     if (mask[i])
0033       continue;
0034
0035     if (track.pt() < ptMin()) {
0036       continue;
0037     }
0038     if (std::abs(track.dxy(regOrigin)) > originRBound()) {
0039       continue;
0040     }
0041     if (std::abs(track.dz(regOrigin)) > originZBound()) {
0042       continue;
0043     }
0044
0045     mask[i] = true;
0046   }
0047 }
0048
0049 TrackingRegion::Hits CosmicTrackingRegion::hits(const SeedingLayerSetsHits::SeedingLayer& layer) const {
0050   TrackingRegion::Hits result;
0051   hits_(layer, result);
0052   return result;
0053 }
0054
0055 template <typename T>
0056 void CosmicTrackingRegion::hits_(const T& layer, TrackingRegion::Hits& result) const {
0057   //get and name collections
0058   //++++++++++++++++++++++++
0059
0060   //detector layer
0061   const DetLayer* detLayer = layer.detLayer();
0062   LogDebug("CosmicTrackingRegion") << "Looking at hits on subdet/layer " << layer.name();
0063   EtaPhiMeasurementEstimator est(0.3, 0.3);
0064
0065   //region
0066   const GlobalPoint vtx = origin();
0067   GlobalVector dir = direction();
0068   LogDebug("CosmicTrackingRegion") << "The initial region characteristics are:"
0069                                    << "\n"
0070                                    << " Origin    = " << origin() << "\n"
0071                                    << " Direction = " << direction() << "\n"
0072                                    << " Eta = " << origin().eta() << "\n"
0073                                    << " Phi = " << origin().phi();
0074
0075   //trajectory state on surface
0076   float phi = dir.phi();
0077   Surface::RotationType rot(sin(phi), -cos(phi), 0, 0, 0, -1, cos(phi), sin(phi), 0);
0078
0079   Plane::PlanePointer surface = Plane::build(vtx, rot);
0080   FreeTrajectoryState fts(GlobalTrajectoryParameters(vtx, dir, 1, theMagneticField_));
0081   TrajectoryStateOnSurface tsos(fts, *surface);
0082   LogDebug("CosmicTrackingRegion") << "The state used to find measurement with the measurement tracker is:\n" << tsos;
0083
0084   //propagator
0085   AnalyticalPropagator prop(theMagneticField_, alongMomentum);
0086
0087   //propagation verification (debug)
0088   //++++++++++++++++++++++++++++++++
0089
0090   //creation of the state
0091   TrajectoryStateOnSurface stateOnLayer = prop.propagate(*tsos.freeState(), detLayer->surface());
0092
0093   //verification of the state
0094   if (stateOnLayer.isValid()) {
0095     LogDebug("CosmicTrackingRegion") << "The initial state propagates to the layer surface: \n"
0096                                      << stateOnLayer << "R   = " << stateOnLayer.globalPosition().perp() << "\n"
0097                                      << "Eta = " << stateOnLayer.globalPosition().eta() << "\n"
0098                                      << "Phi = " << stateOnLayer.globalPosition().phi();
0099
0100   } else {
0101     LogDebug("CosmicTrackingRegion") << "The initial state does not propagate to the layer surface.";
0102   }
0103
0104   //number of compatible dets
0105   typedef DetLayer::DetWithState DetWithState;
0106   vector<DetWithState> compatDets = detLayer->compatibleDets(tsos, prop, est);
0107   LogDebug("CosmicTrackingRegion") << "Compatible dets = " << compatDets.size();
0108
0109   //get hits
0110   //++++++++
0111
0112   //measurement tracker (find hits)
0113   LayerMeasurements lm(theMeasurementTracker_->measurementTracker(), *theMeasurementTracker_);
0114   vector<TrajectoryMeasurement> meas = lm.measurements(*detLayer, tsos, prop, est);
0115   LogDebug("CosmicTrackingRegion") << "Number of Trajectory measurements = " << meas.size()
0116                                    << " but the last one is always an invalid hit, by construction.";
0117
0118   //trajectory measurement
0119
0120   // std::cout <<"CRegion b " << cache.size() << std::endl;
0121
0122   // waiting for a migration at LayerMeasurements level and at seed builder level
0123   for (auto const& im : meas) {
0124     if (!im.recHit()->isValid())
0125       continue;
0126     assert(!trackerHitRTTI::isUndef(*im.recHit()->hit()));
0127     auto ptrHit = (BaseTrackerRecHit*)(im.recHit()->hit()->clone());
0128     cache.emplace_back(ptrHit);
0129     result.emplace_back(ptrHit);
0130   }
0131
0132   // std::cout <<"CRegion a " << cache.size() << std::endl;
0133 }