TkDetLayers/src/Phase2OTBarrelRod.cc

0001 #include "Phase2OTBarrelRod.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 DetGroupElementPerpLess {
0020   public:
0021     bool operator()(DetGroup a, DetGroup b) const {
0022       return (a.front().det()->position().perp() < b.front().det()->position().perp());
0023     }
0024   };
0025 }  // namespace
0026
0027 Phase2OTBarrelRod::Phase2OTBarrelRod(vector<const GeomDet*>& innerDets,
0028                                      vector<const GeomDet*>& outerDets,
0029                                      const vector<const GeomDet*>& innerDetBrothers,
0030                                      const vector<const GeomDet*>& outerDetBrothers)
0031     : DetRod(true),
0032       theInnerDets(innerDets),
0033       theOuterDets(outerDets),
0034       theInnerDetBrothers(innerDetBrothers),
0035       theOuterDetBrothers(outerDetBrothers) {
0036   theDets.assign(theInnerDets.begin(), theInnerDets.end());
0037   theDets.insert(theDets.end(), theOuterDets.begin(), theOuterDets.end());
0038   theDets.insert(theDets.end(), theInnerDetBrothers.begin(), theInnerDetBrothers.end());
0039   theDets.insert(theDets.end(), theOuterDetBrothers.begin(), theOuterDetBrothers.end());
0040
0041   RodPlaneBuilderFromDet planeBuilder;
0042   setPlane(planeBuilder(theDets));
0043   theInnerPlane = planeBuilder(theInnerDets);
0044   theOuterPlane = planeBuilder(theOuterDets);
0045
0046   // modules be already sorted in z
0047
0048   theInnerBinFinder = BinFinderType(theInnerDets.begin(), theInnerDets.end());
0049   theOuterBinFinder = BinFinderType(theOuterDets.begin(), theOuterDets.end());
0050
0051 #ifdef EDM_ML_DEBUG
0052   LogDebug("TkDetLayers") << "==== DEBUG Phase2OTBarrelRod =====";
0053   for (vector<const GeomDet*>::const_iterator i = theInnerDets.begin(); i != theInnerDets.end(); i++) {
0054     LogDebug("TkDetLayers") << "inner Phase2OTBarrelRod's Det pos z,perp,eta,phi: " << (**i).position().z() << " , "
0055                             << (**i).position().perp() << " , " << (**i).position().eta() << " , "
0056                             << (**i).position().phi();
0057   }
0058
0059   for (vector<const GeomDet*>::const_iterator i = theInnerDetBrothers.begin(); i != theInnerDetBrothers.end(); i++) {
0060     LogDebug("TkDetLayers") << "inner Phase2OTBarrelRod's Det Brother pos z,perp,eta,phi: " << (**i).position().z()
0061                             << " , " << (**i).position().perp() << " , " << (**i).position().eta() << " , "
0062                             << (**i).position().phi();
0063   }
0064   if (theInnerDetBrothers.empty() && theOuterDetBrothers.empty())
0065     LogDebug("TkDetLayers") << "====       with stacks       =====";
0066   if (!theInnerDetBrothers.empty() && !theOuterDetBrothers.empty())
0067     LogDebug("TkDetLayers") << "====     without stacks      =====";
0068
0069   for (vector<const GeomDet*>::const_iterator i = theOuterDets.begin(); i != theOuterDets.end(); i++) {
0070     LogDebug("TkDetLayers") << "outer Phase2OTBarrelRod's Det pos z,perp,eta,phi: " << (**i).position().z() << " , "
0071                             << (**i).position().perp() << " , " << (**i).position().eta() << " , "
0072                             << (**i).position().phi();
0073   }
0074
0075   for (vector<const GeomDet*>::const_iterator i = theOuterDetBrothers.begin(); i != theOuterDetBrothers.end(); i++) {
0076     LogDebug("TkDetLayers") << "outer Phase2OTBarrelRod's Det Brother pos z,perp,eta,phi: " << (**i).position().z()
0077                             << " , " << (**i).position().perp() << " , " << (**i).position().eta() << " , "
0078                             << (**i).position().phi();
0079   }
0080   LogDebug("TkDetLayers") << "==== end DEBUG Phase2OTBarrelRod =====";
0081 #endif
0082 }
0083
0084 Phase2OTBarrelRod::~Phase2OTBarrelRod() {}
0085
0086 const vector<const GeometricSearchDet*>& Phase2OTBarrelRod::components() const {
0087   throw DetLayerException("Phase2OTBarrelRod doesn't have GeometricSearchDet components");
0088 }
0089
0090 pair<bool, TrajectoryStateOnSurface> Phase2OTBarrelRod::compatible(const TrajectoryStateOnSurface& ts,
0091                                                                    const Propagator&,
0092                                                                    const MeasurementEstimator&) const {
0093   edm::LogError("TkDetLayers") << "temporary dummy implementation of Phase2OTBarrelRod::compatible()!!";
0094   return pair<bool, TrajectoryStateOnSurface>();
0095 }
0096
0097 void Phase2OTBarrelRod::groupedCompatibleDetsV(const TrajectoryStateOnSurface& tsos,
0098                                                const Propagator& prop,
0099                                                const MeasurementEstimator& est,
0100                                                std::vector<DetGroup>& result) const {
0101   SubLayerCrossings crossings;
0102   crossings = computeCrossings(tsos, prop.propagationDirection());
0103   if (!crossings.isValid())
0104     return;
0105
0106   std::vector<DetGroup> closestResult;
0107   std::vector<DetGroup> closestBrotherResult;
0108   addClosest(tsos, prop, est, crossings.closest(), closestResult, closestBrotherResult);
0109   if (closestResult.empty()) {
0110     std::vector<DetGroup> nextResult;
0111     std::vector<DetGroup> nextBrotherResult;
0112     addClosest(tsos, prop, est, crossings.other(), nextResult, nextBrotherResult);
0113     if (nextResult.empty())
0114       return;
0115
0116     DetGroupElement nextGel(nextResult.front().front());
0117     int crossingSide = LayerCrossingSide().barrelSide(nextGel.trajectoryState(), prop);
0118     std::vector<DetGroup> closestCompleteResult;
0119     DetGroupMerger::orderAndMergeTwoLevels(
0120         std::move(closestResult), std::move(closestBrotherResult), closestCompleteResult, 0, crossingSide);
0121     std::vector<DetGroup> nextCompleteResult;
0122     DetGroupMerger::orderAndMergeTwoLevels(
0123         std::move(nextResult), std::move(nextBrotherResult), nextCompleteResult, 0, crossingSide);
0124
0125     DetGroupMerger::orderAndMergeTwoLevels(
0126         std::move(closestCompleteResult), std::move(nextCompleteResult), result, crossings.closestIndex(), crossingSide);
0127   } else {
0128     DetGroupElement closestGel(closestResult.front().front());
0129     int crossingSide = LayerCrossingSide().barrelSide(closestGel.trajectoryState(), prop);
0130     float window = computeWindowSize(closestGel.det(), closestGel.trajectoryState(), est);
0131
0132     searchNeighbors(tsos, prop, est, crossings.closest(), window, closestResult, closestBrotherResult, false);
0133
0134     std::vector<DetGroup> closestCompleteResult;
0135     DetGroupMerger::orderAndMergeTwoLevels(
0136         std::move(closestResult), std::move(closestBrotherResult), closestCompleteResult, 0, crossingSide);
0137
0138     std::vector<DetGroup> nextResult;
0139     std::vector<DetGroup> nextBrotherResult;
0140     searchNeighbors(tsos, prop, est, crossings.other(), window, nextResult, nextBrotherResult, true);
0141
0142     std::vector<DetGroup> nextCompleteResult;
0143     DetGroupMerger::orderAndMergeTwoLevels(
0144         std::move(nextResult), std::move(nextBrotherResult), nextCompleteResult, 0, crossingSide);
0145
0146     DetGroupMerger::orderAndMergeTwoLevels(
0147         std::move(closestCompleteResult), std::move(nextCompleteResult), result, crossings.closestIndex(), crossingSide);
0148   }
0149
0150   //due to propagator problems, when we add single pt sub modules, we should order them in r (barrel)
0151   sort(result.begin(), result.end(), DetGroupElementPerpLess());
0152   for (auto& grp : result) {
0153     if (grp.empty())
0154       continue;
0155     LogTrace("TkDetLayers") << "New group in Phase2OTBarrelRod made by : ";
0156     for (auto const& det : grp) {
0157       LogTrace("TkDetLayers") << " geom det at r: " << det.det()->position().perp()
0158                               << " id:" << det.det()->geographicalId().rawId()
0159                               << " tsos at:" << det.trajectoryState().globalPosition();
0160     }
0161   }
0162 }
0163
0164 SubLayerCrossings Phase2OTBarrelRod::computeCrossings(const TrajectoryStateOnSurface& startingState,
0165                                                       PropagationDirection propDir) const {
0166   GlobalPoint startPos(startingState.globalPosition());
0167   GlobalVector startDir(startingState.globalMomentum());
0168   double rho(startingState.transverseCurvature());
0169
0170   HelixBarrelPlaneCrossingByCircle crossing(startPos, startDir, rho, propDir);
0171
0172   std::pair<bool, double> outerPath = crossing.pathLength(*theOuterPlane);
0173   if (!outerPath.first)
0174     return SubLayerCrossings();
0175   GlobalPoint gOuterPoint(crossing.position(outerPath.second));
0176
0177   std::pair<bool, double> innerPath = crossing.pathLength(*theInnerPlane);
0178   if (!innerPath.first)
0179     return SubLayerCrossings();
0180   GlobalPoint gInnerPoint(crossing.position(innerPath.second));
0181
0182   int innerIndex = theInnerBinFinder.binIndex(gInnerPoint.z());
0183   float innerDist = std::abs(theInnerBinFinder.binPosition(innerIndex) - gInnerPoint.z());
0184   SubLayerCrossing innerSLC(0, innerIndex, gInnerPoint);
0185
0186   int outerIndex = theOuterBinFinder.binIndex(gOuterPoint.z());
0187   float outerDist = std::abs(theOuterBinFinder.binPosition(outerIndex) - gOuterPoint.z());
0188   SubLayerCrossing outerSLC(1, outerIndex, gOuterPoint);
0189
0190   if (innerDist < outerDist) {
0191     return SubLayerCrossings(innerSLC, outerSLC, 0);
0192   } else {
0193     return SubLayerCrossings(outerSLC, innerSLC, 1);
0194   }
0195 }
0196
0197 bool Phase2OTBarrelRod::addClosest(const TrajectoryStateOnSurface& tsos,
0198                                    const Propagator& prop,
0199                                    const MeasurementEstimator& est,
0200                                    const SubLayerCrossing& crossing,
0201                                    vector<DetGroup>& result,
0202                                    vector<DetGroup>& brotherresult) const {
0203   const vector<const GeomDet*>& sRod(subRod(crossing.subLayerIndex()));
0204   bool firstgroup = CompatibleDetToGroupAdder::add(*sRod[crossing.closestDetIndex()], tsos, prop, est, result);
0205   if (theInnerDetBrothers.empty() && theOuterDetBrothers.empty())
0206     return firstgroup;
0207
0208   // it assumes that the closestDetIndex is ok also for the brother detectors: the crossing is NOT recomputed
0209   const vector<const GeomDet*>& sRodBrothers(subRodBrothers(crossing.subLayerIndex()));
0210   bool brothergroup =
0211       CompatibleDetToGroupAdder::add(*sRodBrothers[crossing.closestDetIndex()], tsos, prop, est, brotherresult);
0212
0213   return firstgroup || brothergroup;
0214 }
0215
0216 float Phase2OTBarrelRod::computeWindowSize(const GeomDet* det,
0217                                            const TrajectoryStateOnSurface& tsos,
0218                                            const MeasurementEstimator& est) const {
0219   return est.maximalLocalDisplacement(tsos, det->surface()).y();
0220 }
0221
0222 namespace {
0223
0224   inline bool overlap(const GlobalPoint& crossPoint, const GeomDet& det, float window) {
0225     // check if the z window around TSOS overlaps with the detector theDet (with a 1% margin added)
0226
0227     //   const float tolerance = 0.1;
0228     constexpr float relativeMargin = 1.01;
0229
0230     LocalPoint localCrossPoint(det.surface().toLocal(crossPoint));
0231     //   if (std::abs(localCrossPoint.z()) > tolerance) {
0232     //     edm::LogInfo(TkDetLayers) << "TOBRod::overlap calculation assumes point on surface, but it is off by "
0233     //   << localCrossPoint.z() ;
0234     //   }
0235
0236     float localY = localCrossPoint.y();
0237     float detHalfLength = 0.5f * det.surface().bounds().length();
0238
0239     //   edm::LogInfo(TkDetLayers) << "TOBRod::overlap: Det at " << det.position() << " hit at " << localY
0240     //        << " Window " << window << " halflength "  << detHalfLength ;
0241
0242     return (std::abs(localY) - window) < relativeMargin * detHalfLength;
0243   }
0244
0245 }  // namespace
0246
0247 void Phase2OTBarrelRod::searchNeighbors(const TrajectoryStateOnSurface& tsos,
0248                                         const Propagator& prop,
0249                                         const MeasurementEstimator& est,
0250                                         const SubLayerCrossing& crossing,
0251                                         float window,
0252                                         vector<DetGroup>& result,
0253                                         vector<DetGroup>& brotherresult,
0254                                         bool checkClosest) const {
0255   const GlobalPoint& gCrossingPos = crossing.position();
0256
0257   const vector<const GeomDet*>& sRod(subRod(crossing.subLayerIndex()));
0258   const vector<const GeomDet*>& sBrotherRod(subRodBrothers(crossing.subLayerIndex()));
0259
0260   int closestIndex = crossing.closestDetIndex();
0261   int negStartIndex = closestIndex - 1;
0262   int posStartIndex = closestIndex + 1;
0263
0264   if (checkClosest) {  // must decide if the closest is on the neg or pos side
0265     if (gCrossingPos.z() < sRod[closestIndex]->surface().position().z()) {
0266       posStartIndex = closestIndex;
0267     } else {
0268       negStartIndex = closestIndex;
0269     }
0270   }
0271
0272   typedef CompatibleDetToGroupAdder Adder;
0273   for (int idet = negStartIndex; idet >= 0; idet--) {
0274     if (!overlap(gCrossingPos, *sRod[idet], window))
0275       break;
0276     if (!Adder::add(*sRod[idet], tsos, prop, est, result))
0277       break;
0278     if (theInnerDetBrothers.empty() && theOuterDetBrothers.empty())
0279       break;
0280     // If the two above checks are passed also the brother module will be added with no further checks
0281     Adder::add(*sBrotherRod[idet], tsos, prop, est, brotherresult);
0282   }
0283   for (int idet = posStartIndex; idet < static_cast<int>(sRod.size()); idet++) {
0284     if (!overlap(gCrossingPos, *sRod[idet], window))
0285       break;
0286     if (!Adder::add(*sRod[idet], tsos, prop, est, result))
0287       break;
0288     if (theInnerDetBrothers.empty() && theOuterDetBrothers.empty())
0289       break;
0290     // If the two above checks are passed also the brother module will be added with no further checks
0291     Adder::add(*sBrotherRod[idet], tsos, prop, est, brotherresult);
0292   }
0293 }