TkDetLayers/src/Phase2EndcapRing.cc

0001 #include "Phase2EndcapRing.h"
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 #include "TrackingTools/DetLayers/interface/DetLayerException.h"
0006 #include "TrackingTools/DetLayers/interface/MeasurementEstimator.h"
0007 #include "TrackingTools/GeomPropagators/interface/HelixForwardPlaneCrossing.h"
0008 #include "TrackingTools/DetLayers/interface/rangesIntersect.h"
0009 #include "TrackingTools/DetLayers/interface/ForwardRingDiskBuilderFromDet.h"
0010
0011 #include "LayerCrossingSide.h"
0012 #include "DetGroupMerger.h"
0013 #include "CompatibleDetToGroupAdder.h"
0014
0015 #include "TkDetUtil.h"
0016 #include "DataFormats/GeometryVector/interface/VectorUtil.h"
0017
0018 #include "DetGroupElementZLess.h"
0019
0020 using namespace std;
0021
0022 typedef GeometricSearchDet::DetWithState DetWithState;
0023
0024 Phase2EndcapRing::Phase2EndcapRing(vector<const GeomDet*>& innerDets,
0025                                    vector<const GeomDet*>& outerDets,
0026                                    const vector<const GeomDet*>& innerDetBrothers,
0027                                    const vector<const GeomDet*>& outerDetBrothers)
0028     : GeometricSearchDet(true),
0029       theFrontDets(innerDets.begin(), innerDets.end()),
0030       theBackDets(outerDets.begin(), outerDets.end()),
0031       theFrontDetBrothers(innerDetBrothers.begin(), innerDetBrothers.end()),
0032       theBackDetBrothers(outerDetBrothers.begin(), outerDetBrothers.end()) {
0033   theDets.assign(theFrontDets.begin(), theFrontDets.end());
0034   theDets.insert(theDets.end(), theBackDets.begin(), theBackDets.end());
0035   theDets.insert(theDets.end(), theFrontDetBrothers.begin(), theFrontDetBrothers.end());
0036   theDets.insert(theDets.end(), theBackDetBrothers.begin(), theBackDetBrothers.end());
0037
0038   // the dets should be already phi-ordered. TO BE CHECKED
0039   //sort( theFrontDets.begin(), theFrontDets.end(), DetLessPhi() );
0040   //sort( theBackDets.begin(), theBackDets.end(), DetLessPhi() );
0041
0042   theDisk = ForwardRingDiskBuilderFromDet()(theDets);
0043
0044   theFrontDisk = ForwardRingDiskBuilderFromDet()(theFrontDets);
0045   theBackDisk = ForwardRingDiskBuilderFromDet()(theBackDets);
0046
0047   theFrontBinFinder = BinFinderType(theFrontDets.front()->surface().position().phi(), theFrontDets.size());
0048   theBackBinFinder = BinFinderType(theBackDets.front()->surface().position().phi(), theBackDets.size());
0049
0050 #ifdef EDM_ML_DEBUG
0051   LogDebug("TkDetLayers") << "DEBUG INFO for Phase2EndcapRing";
0052   for (vector<const GeomDet*>::const_iterator it = theFrontDets.begin(); it != theFrontDets.end(); it++) {
0053     LogDebug("TkDetLayers") << "frontDet detId,phi,z,r: " << (*it)->geographicalId().rawId() << " , "
0054                             << (*it)->surface().position().phi() << " , " << (*it)->surface().position().z() << " , "
0055                             << (*it)->surface().position().perp();
0056   }
0057
0058   if (!theFrontDetBrothers.empty()) {
0059     for (vector<const GeomDet*>::const_iterator it = theFrontDetBrothers.begin(); it != theFrontDetBrothers.end();
0060          it++) {
0061       LogDebug("TkDetLayers") << "frontDet brothers detId,phi,z,r: " << (*it)->geographicalId().rawId() << " , "
0062                               << (*it)->surface().position().phi() << " , " << (*it)->surface().position().z() << " , "
0063                               << (*it)->surface().position().perp();
0064     }
0065   }
0066
0067   for (vector<const GeomDet*>::const_iterator it = theBackDets.begin(); it != theBackDets.end(); it++) {
0068     LogDebug("TkDetLayers") << "backDet detId,phi,z,r: " << (*it)->geographicalId().rawId() << " , "
0069                             << (*it)->surface().position().phi() << " , " << (*it)->surface().position().z() << " , "
0070                             << (*it)->surface().position().perp();
0071   }
0072
0073   if (!theBackDetBrothers.empty()) {
0074     for (vector<const GeomDet*>::const_iterator it = theBackDetBrothers.begin(); it != theBackDetBrothers.end(); it++) {
0075       LogDebug("TkDetLayers") << "backDet brothers detId,phi,z,r: " << (*it)->geographicalId().rawId() << " , "
0076                               << (*it)->surface().position().phi() << " , " << (*it)->surface().position().z() << " , "
0077                               << (*it)->surface().position().perp();
0078     }
0079   }
0080 #endif
0081 }
0082
0083 Phase2EndcapRing::~Phase2EndcapRing() {}
0084
0085 const vector<const GeometricSearchDet*>& Phase2EndcapRing::components() const {
0086   throw DetLayerException("Phase2EndcapRing doesn't have GeometricSearchDet components");
0087 }
0088
0089 pair<bool, TrajectoryStateOnSurface> Phase2EndcapRing::compatible(const TrajectoryStateOnSurface&,
0090                                                                   const Propagator&,
0091                                                                   const MeasurementEstimator&) const {
0092   edm::LogError("TkDetLayers") << "temporary dummy implementation of Phase2EndcapRing::compatible()!!";
0093   return pair<bool, TrajectoryStateOnSurface>();
0094 }
0095
0096 void Phase2EndcapRing::groupedCompatibleDetsV(const TrajectoryStateOnSurface& tsos,
0097                                               const Propagator& prop,
0098                                               const MeasurementEstimator& est,
0099                                               std::vector<DetGroup>& result) const {
0100   SubLayerCrossings crossings;
0101   crossings = computeCrossings(tsos, prop.propagationDirection());
0102   if (!crossings.isValid())
0103     return;
0104
0105   std::vector<DetGroup> closestResult;
0106   std::vector<DetGroup> closestBrotherResult;
0107   addClosest(tsos, prop, est, crossings.closest(), closestResult, closestBrotherResult);
0108   if (closestResult.empty())
0109     return;
0110
0111   DetGroupElement closestGel(closestResult.front().front());
0112   int crossingSide = LayerCrossingSide().endcapSide(closestGel.trajectoryState(), prop);
0113   float phiWindow = tkDetUtil::computeWindowSize(closestGel.det(), closestGel.trajectoryState(), est);
0114   searchNeighbors(tsos, prop, est, crossings.closest(), phiWindow, closestResult, closestBrotherResult, false);
0115
0116   vector<DetGroup> closestCompleteResult;
0117   DetGroupMerger::orderAndMergeTwoLevels(
0118       std::move(closestResult), std::move(closestBrotherResult), closestCompleteResult, 0, crossingSide);
0119
0120   vector<DetGroup> nextResult;
0121   vector<DetGroup> nextBrotherResult;
0122   searchNeighbors(tsos, prop, est, crossings.other(), phiWindow, nextResult, nextBrotherResult, true);
0123
0124   vector<DetGroup> nextCompleteResult;
0125   DetGroupMerger::orderAndMergeTwoLevels(
0126       std::move(nextResult), std::move(nextBrotherResult), nextCompleteResult, 0, crossingSide);
0127
0128   DetGroupMerger::orderAndMergeTwoLevels(
0129       std::move(closestCompleteResult), std::move(nextCompleteResult), result, crossings.closestIndex(), crossingSide);
0130
0131   //due to propagator problems, when we add single pt sub modules, we should order them in z (endcap)
0132   if (!theFrontDetBrothers.empty() && !theBackDetBrothers.empty())
0133     sort(result.begin(), result.end(), DetGroupElementZLess());
0134
0135 #ifdef EDM_ML_DEBUG
0136   LogTrace("TkDetLayers") << "Number of groups : " << result.size() << std::endl;
0137   for (auto& grp : result) {
0138     if (grp.empty())
0139       continue;
0140     LogTrace("TkDetLayers") << "New group in Phase2EndcapRing made by : " << std::endl;
0141     for (auto const& det : grp) {
0142       LogTrace("TkDetLayers") << " geom det at r: " << det.det()->position().perp()
0143                               << " id:" << det.det()->geographicalId().rawId()
0144                               << " tsos at:" << det.trajectoryState().globalPosition() << std::endl;
0145     }
0146   }
0147 #endif
0148 }
0149
0150 SubLayerCrossings Phase2EndcapRing::computeCrossings(const TrajectoryStateOnSurface& startingState,
0151                                                      PropagationDirection propDir) const {
0152   auto rho = startingState.transverseCurvature();
0153
0154   HelixPlaneCrossing::PositionType startPos(startingState.globalPosition());
0155   HelixPlaneCrossing::DirectionType startDir(startingState.globalMomentum());
0156   HelixForwardPlaneCrossing crossing(startPos, startDir, rho, propDir);
0157
0158   pair<bool, double> frontPath = crossing.pathLength(*theFrontDisk);
0159   if (!frontPath.first)
0160     return SubLayerCrossings();
0161
0162   pair<bool, double> backPath = crossing.pathLength(*theBackDisk);
0163   if (!backPath.first)
0164     return SubLayerCrossings();
0165
0166   GlobalPoint gFrontPoint(crossing.position(frontPath.second));
0167   GlobalPoint gBackPoint(crossing.position(backPath.second));
0168
0169   int frontIndex = theFrontBinFinder.binIndex(gFrontPoint.barePhi());
0170   SubLayerCrossing frontSLC(0, frontIndex, gFrontPoint);
0171
0172   int backIndex = theBackBinFinder.binIndex(gBackPoint.barePhi());
0173   SubLayerCrossing backSLC(1, backIndex, gBackPoint);
0174
0175   // 0ss: frontDisk has index=0, backDisk has index=1
0176   float frontDist = std::abs(Geom::deltaPhi(gFrontPoint.barePhi(), theFrontDets[frontIndex]->surface().phi()));
0177   float backDist = std::abs(Geom::deltaPhi(gBackPoint.barePhi(), theBackDets[backIndex]->surface().phi()));
0178
0179   if (frontDist < backDist) {
0180     return SubLayerCrossings(frontSLC, backSLC, 0);
0181   } else {
0182     return SubLayerCrossings(backSLC, frontSLC, 1);
0183   }
0184 }
0185
0186 bool Phase2EndcapRing::addClosest(const TrajectoryStateOnSurface& tsos,
0187                                   const Propagator& prop,
0188                                   const MeasurementEstimator& est,
0189                                   const SubLayerCrossing& crossing,
0190                                   vector<DetGroup>& result,
0191                                   vector<DetGroup>& brotherresult) const {
0192   const vector<const GeomDet*>& sub(subLayer(crossing.subLayerIndex()));
0193   const GeomDet* det(sub[crossing.closestDetIndex()]);
0194   bool firstgroup = CompatibleDetToGroupAdder::add(*det, tsos, prop, est, result);
0195   if (theFrontDetBrothers.empty() && theBackDetBrothers.empty())
0196     return firstgroup;
0197   // it assumes that the closestDetIndex is ok also for the brother detectors: the crossing is NOT recomputed
0198   const vector<const GeomDet*>& subBrothers(subLayerBrothers(crossing.subLayerIndex()));
0199   const GeomDet* detBrother(subBrothers[crossing.closestDetIndex()]);
0200   bool brothergroup = CompatibleDetToGroupAdder::add(*detBrother, tsos, prop, est, brotherresult);
0201   return firstgroup || brothergroup;
0202 }
0203
0204 void Phase2EndcapRing::searchNeighbors(const TrajectoryStateOnSurface& tsos,
0205                                        const Propagator& prop,
0206                                        const MeasurementEstimator& est,
0207                                        const SubLayerCrossing& crossing,
0208                                        float window,
0209                                        vector<DetGroup>& result,
0210                                        vector<DetGroup>& brotherresult,
0211                                        bool checkClosest) const {
0212   const GlobalPoint& gCrossingPos = crossing.position();
0213
0214   const vector<const GeomDet*>& sLayer(subLayer(crossing.subLayerIndex()));
0215   // It assumes that what is ok for the front modules in the pt modules is ok also for the back module
0216   const vector<const GeomDet*>& sBrotherLayer(subLayerBrothers(crossing.subLayerIndex()));
0217
0218   int closestIndex = crossing.closestDetIndex();
0219   int negStartIndex = closestIndex - 1;
0220   int posStartIndex = closestIndex + 1;
0221
0222   if (checkClosest) {  // must decide if the closest is on the neg or pos side
0223     if (Geom::phiLess(gCrossingPos.barePhi(), sLayer[closestIndex]->surface().phi())) {
0224       posStartIndex = closestIndex;
0225     } else {
0226       negStartIndex = closestIndex;
0227     }
0228   }
0229
0230   const BinFinderType& binFinder = (crossing.subLayerIndex() == 0 ? theFrontBinFinder : theBackBinFinder);
0231
0232   typedef CompatibleDetToGroupAdder Adder;
0233   int half = sLayer.size() / 2;  // to check if dets are called twice....
0234   for (int idet = negStartIndex; idet >= negStartIndex - half; idet--) {
0235     const GeomDet& neighborDet = *sLayer[binFinder.binIndex(idet)];
0236     if (!tkDetUtil::overlapInPhi(gCrossingPos, neighborDet, window))
0237       break;
0238     if (!Adder::add(neighborDet, tsos, prop, est, result))
0239       break;
0240     if (theFrontDetBrothers.empty() && theBackDetBrothers.empty())
0241       break;
0242     // If the two above checks are passed also the brother module will be added with no further checks
0243     const GeomDet& neighborBrotherDet = *sBrotherLayer[binFinder.binIndex(idet)];
0244     Adder::add(neighborBrotherDet, tsos, prop, est, brotherresult);
0245     // maybe also add shallow crossing angle test here???
0246   }
0247   for (int idet = posStartIndex; idet < posStartIndex + half; idet++) {
0248     const GeomDet& neighborDet = *sLayer[binFinder.binIndex(idet)];
0249     if (!tkDetUtil::overlapInPhi(gCrossingPos, neighborDet, window))
0250       break;
0251     if (!Adder::add(neighborDet, tsos, prop, est, result))
0252       break;
0253     if (theFrontDetBrothers.empty() && theBackDetBrothers.empty())
0254       break;
0255     // If the two above checks are passed also the brother module will be added with no further checks
0256     const GeomDet& neighborBrotherDet = *sBrotherLayer[binFinder.binIndex(idet)];
0257     Adder::add(neighborBrotherDet, tsos, prop, est, brotherresult);
0258     // maybe also add shallow crossing angle test here???
0259   }
0260 }