TkDetLayers/src/TOBRod.cc

0001 #include "TOBRod.h"
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 #include "TrackingTools/DetLayers/interface/RodPlaneBuilderFromDet.h"
0006 #include "TrackingTools/DetLayers/interface/DetLayerException.h"
0007 #include "TrackingTools/DetLayers/interface/MeasurementEstimator.h"
0008 #include "TrackingTools/GeomPropagators/interface/HelixBarrelPlaneCrossingByCircle.h"
0009
0010 #include "LayerCrossingSide.h"
0011 #include "DetGroupMerger.h"
0012 #include "CompatibleDetToGroupAdder.h"
0013
0014 using namespace std;
0015
0016 typedef GeometricSearchDet::DetWithState DetWithState;
0017
0018 namespace {
0019   class DetZLess {
0020   public:
0021     bool operator()(const GeomDet* a, const GeomDet* b) const { return (a->position().z() < b->position().z()); }
0022   };
0023 }  // namespace
0024
0025 TOBRod::TOBRod(vector<const GeomDet*>& innerDets, vector<const GeomDet*>& outerDets)
0026     : DetRod(true), theInnerDets(innerDets), theOuterDets(outerDets) {
0027   theDets.assign(theInnerDets.begin(), theInnerDets.end());
0028   theDets.insert(theDets.end(), theOuterDets.begin(), theOuterDets.end());
0029
0030   RodPlaneBuilderFromDet planeBuilder;
0031   setPlane(planeBuilder(theDets));
0032   theInnerPlane = planeBuilder(theInnerDets);
0033   theOuterPlane = planeBuilder(theOuterDets);
0034
0035   sort(theDets.begin(), theDets.end(), DetZLess());
0036   sort(theInnerDets.begin(), theInnerDets.end(), DetZLess());
0037   sort(theOuterDets.begin(), theOuterDets.end(), DetZLess());
0038   theInnerBinFinder = BinFinderType(theInnerDets.begin(), theInnerDets.end());
0039   theOuterBinFinder = BinFinderType(theOuterDets.begin(), theOuterDets.end());
0040
0041   LogDebug("TkDetLayers") << "==== DEBUG TOBRod =====";
0042   for (vector<const GeomDet*>::const_iterator i = theInnerDets.begin(); i != theInnerDets.end(); i++) {
0043     LogDebug("TkDetLayers") << "inner TOBRod's Det pos z,perp,eta,phi: " << (**i).position().z() << " , "
0044                             << (**i).position().perp() << " , " << (**i).position().eta() << " , "
0045                             << (**i).position().phi();
0046   }
0047
0048   for (vector<const GeomDet*>::const_iterator i = theOuterDets.begin(); i != theOuterDets.end(); i++) {
0049     LogDebug("TkDetLayers") << "outer TOBRod's Det pos z,perp,eta,phi: " << (**i).position().z() << " , "
0050                             << (**i).position().perp() << " , " << (**i).position().eta() << " , "
0051                             << (**i).position().phi();
0052   }
0053   LogDebug("TkDetLayers") << "==== end DEBUG TOBRod =====";
0054 }
0055
0056 TOBRod::~TOBRod() {}
0057
0058 const vector<const GeometricSearchDet*>& TOBRod::components() const {
0059   throw DetLayerException("TOBRod doesn't have GeometricSearchDet components");
0060 }
0061
0062 pair<bool, TrajectoryStateOnSurface> TOBRod::compatible(const TrajectoryStateOnSurface& ts,
0063                                                         const Propagator&,
0064                                                         const MeasurementEstimator&) const {
0065   edm::LogError("TkDetLayers") << "temporary dummy implementation of TOBRod::compatible()!!";
0066   return pair<bool, TrajectoryStateOnSurface>();
0067 }
0068
0069 void TOBRod::groupedCompatibleDetsV(const TrajectoryStateOnSurface& tsos,
0070                                     const Propagator& prop,
0071                                     const MeasurementEstimator& est,
0072                                     std::vector<DetGroup>& result) const {
0073   SubLayerCrossings crossings;
0074   crossings = computeCrossings(tsos, prop.propagationDirection());
0075   if (!crossings.isValid())
0076     return;
0077
0078   std::vector<DetGroup> closestResult;
0079   addClosest(tsos, prop, est, crossings.closest(), closestResult);
0080   if (closestResult.empty()) {
0081     std::vector<DetGroup> nextResult;
0082     addClosest(tsos, prop, est, crossings.other(), nextResult);
0083     if (nextResult.empty())
0084       return;
0085
0086     DetGroupElement nextGel(nextResult.front().front());
0087     int crossingSide = LayerCrossingSide().barrelSide(nextGel.trajectoryState(), prop);
0088     DetGroupMerger::orderAndMergeTwoLevels(
0089         std::move(closestResult), std::move(nextResult), result, crossings.closestIndex(), crossingSide);
0090   } else {
0091     DetGroupElement closestGel(closestResult.front().front());
0092     float window = computeWindowSize(closestGel.det(), closestGel.trajectoryState(), est);
0093
0094     searchNeighbors(tsos, prop, est, crossings.closest(), window, closestResult, false);
0095
0096     std::vector<DetGroup> nextResult;
0097     searchNeighbors(tsos, prop, est, crossings.other(), window, nextResult, true);
0098
0099     int crossingSide = LayerCrossingSide().barrelSide(closestGel.trajectoryState(), prop);
0100     DetGroupMerger::orderAndMergeTwoLevels(
0101         std::move(closestResult), std::move(nextResult), result, crossings.closestIndex(), crossingSide);
0102   }
0103 }
0104
0105 SubLayerCrossings TOBRod::computeCrossings(const TrajectoryStateOnSurface& startingState,
0106                                            PropagationDirection propDir) const {
0107   GlobalPoint startPos(startingState.globalPosition());
0108   GlobalVector startDir(startingState.globalMomentum());
0109   double rho(startingState.transverseCurvature());
0110
0111   HelixBarrelPlaneCrossingByCircle crossing(startPos, startDir, rho, propDir);
0112
0113   std::pair<bool, double> outerPath = crossing.pathLength(*theOuterPlane);
0114   if (!outerPath.first)
0115     return SubLayerCrossings();
0116   GlobalPoint gOuterPoint(crossing.position(outerPath.second));
0117
0118   std::pair<bool, double> innerPath = crossing.pathLength(*theInnerPlane);
0119   if (!innerPath.first)
0120     return SubLayerCrossings();
0121   GlobalPoint gInnerPoint(crossing.position(innerPath.second));
0122
0123   int innerIndex = theInnerBinFinder.binIndex(gInnerPoint.z());
0124   float innerDist = std::abs(theInnerBinFinder.binPosition(innerIndex) - gInnerPoint.z());
0125   SubLayerCrossing innerSLC(0, innerIndex, gInnerPoint);
0126
0127   int outerIndex = theOuterBinFinder.binIndex(gOuterPoint.z());
0128   float outerDist = std::abs(theOuterBinFinder.binPosition(outerIndex) - gOuterPoint.z());
0129   SubLayerCrossing outerSLC(1, outerIndex, gOuterPoint);
0130
0131   if (innerDist < outerDist) {
0132     return SubLayerCrossings(innerSLC, outerSLC, 0);
0133   } else {
0134     return SubLayerCrossings(outerSLC, innerSLC, 1);
0135   }
0136 }
0137
0138 bool TOBRod::addClosest(const TrajectoryStateOnSurface& tsos,
0139                         const Propagator& prop,
0140                         const MeasurementEstimator& est,
0141                         const SubLayerCrossing& crossing,
0142                         vector<DetGroup>& result) const {
0143   const vector<const GeomDet*>& sRod(subRod(crossing.subLayerIndex()));
0144   return CompatibleDetToGroupAdder::add(*sRod[crossing.closestDetIndex()], tsos, prop, est, result);
0145 }
0146
0147 float TOBRod::computeWindowSize(const GeomDet* det,
0148                                 const TrajectoryStateOnSurface& tsos,
0149                                 const MeasurementEstimator& est) const {
0150   return est.maximalLocalDisplacement(tsos, det->surface()).y();
0151 }
0152
0153 namespace {
0154
0155   inline bool overlap(const GlobalPoint& crossPoint, const GeomDet& det, float window) {
0156     // check if the z window around TSOS overlaps with the detector theDet (with a 1% margin added)
0157
0158     //   const float tolerance = 0.1;
0159     constexpr float relativeMargin = 1.01;
0160
0161     LocalPoint localCrossPoint(det.surface().toLocal(crossPoint));
0162     //   if (std::abs(localCrossPoint.z()) > tolerance) {
0163     //     edm::LogInfo(TkDetLayers) << "TOBRod::overlap calculation assumes point on surface, but it is off by "
0164     //   << localCrossPoint.z() ;
0165     //   }
0166
0167     float localY = localCrossPoint.y();
0168     float detHalfLength = 0.5f * det.surface().bounds().length();
0169
0170     //   edm::LogInfo(TkDetLayers) << "TOBRod::overlap: Det at " << det.position() << " hit at " << localY
0171     //        << " Window " << window << " halflength "  << detHalfLength ;
0172
0173     return (std::abs(localY) - window) < relativeMargin * detHalfLength;
0174   }
0175
0176 }  // namespace
0177
0178 void TOBRod::searchNeighbors(const TrajectoryStateOnSurface& tsos,
0179                              const Propagator& prop,
0180                              const MeasurementEstimator& est,
0181                              const SubLayerCrossing& crossing,
0182                              float window,
0183                              vector<DetGroup>& result,
0184                              bool checkClosest) const {
0185   const GlobalPoint& gCrossingPos = crossing.position();
0186
0187   const vector<const GeomDet*>& sRod(subRod(crossing.subLayerIndex()));
0188
0189   int closestIndex = crossing.closestDetIndex();
0190   int negStartIndex = closestIndex - 1;
0191   int posStartIndex = closestIndex + 1;
0192
0193   if (checkClosest) {  // must decide if the closest is on the neg or pos side
0194     if (gCrossingPos.z() < sRod[closestIndex]->surface().position().z()) {
0195       posStartIndex = closestIndex;
0196     } else {
0197       negStartIndex = closestIndex;
0198     }
0199   }
0200
0201   typedef CompatibleDetToGroupAdder Adder;
0202   for (int idet = negStartIndex; idet >= 0; idet--) {
0203     if (!overlap(gCrossingPos, *sRod[idet], window))
0204       break;
0205     if (!Adder::add(*sRod[idet], tsos, prop, est, result))
0206       break;
0207   }
0208   for (int idet = posStartIndex; idet < static_cast<int>(sRod.size()); idet++) {
0209     if (!overlap(gCrossingPos, *sRod[idet], window))
0210       break;
0211     if (!Adder::add(*sRod[idet], tsos, prop, est, result))
0212       break;
0213   }
0214 }