TkDetLayers/src/TECLayer.cc

0001 #include "TECLayer.h"
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 #include "CompatibleDetToGroupAdder.h"
0006 #include "DetGroupMerger.h"
0007 #include "LayerCrossingSide.h"
0008
0009 #include "TrackingTools/DetLayers/interface/DetLayerException.h"
0010 #include "TrackingTools/DetLayers/interface/rangesIntersect.h"
0011 #include "TrackingTools/GeomPropagators/interface/HelixForwardPlaneCrossing.h"
0012 #include "DataFormats/GeometrySurface/interface/SimpleDiskBounds.h"
0013 #include "DataFormats/GeometryVector/interface/VectorUtil.h"
0014 #include "TkDetUtil.h"
0015
0016 using namespace std;
0017
0018 typedef GeometricSearchDet::DetWithState DetWithState;
0019
0020 namespace {
0021
0022   template <typename T>
0023   BoundDisk* computeDisk(vector<const T*>& petals) {
0024     // Attention: it is assumed that the petals do belong to one layer, and are all
0025     // of the same rmin/rmax extension !!
0026
0027     const BoundDiskSector& diskSector = static_cast<const BoundDiskSector&>(petals.front()->surface());
0028
0029     float rmin = diskSector.innerRadius();
0030     float rmax = diskSector.outerRadius();
0031
0032     float theZmax(petals.front()->position().z());
0033     float theZmin(theZmax);
0034     for (auto i = petals.begin(); i != petals.end(); i++) {
0035       float zmin = (**i).position().z() - (**i).surface().bounds().thickness() / 2.;
0036       float zmax = (**i).position().z() + (**i).surface().bounds().thickness() / 2.;
0037       theZmax = max(theZmax, zmax);
0038       theZmin = min(theZmin, zmin);
0039     }
0040
0041     float zPos = (theZmax + theZmin) / 2.;
0042     Surface::PositionType pos(0., 0., zPos);
0043     Surface::RotationType rot;
0044
0045     return new BoundDisk(pos, rot, new SimpleDiskBounds(rmin, rmax, theZmin - zPos, theZmax - zPos));
0046   }
0047
0048 }  // namespace
0049
0050 TECLayer::TECLayer(vector<const TECPetal*>& innerPetals, vector<const TECPetal*>& outerPetals)
0051     : ForwardDetLayer(true),
0052       theFrontComps(innerPetals.begin(), innerPetals.end()),
0053       theBackComps(outerPetals.begin(), outerPetals.end()) {
0054   theComps.assign(theFrontComps.begin(), theFrontComps.end());
0055   theComps.insert(theComps.end(), theBackComps.begin(), theBackComps.end());
0056
0057   for (vector<const GeometricSearchDet*>::const_iterator it = theComps.begin(); it != theComps.end(); it++) {
0058     theBasicComps.insert(theBasicComps.end(), (**it).basicComponents().begin(), (**it).basicComponents().end());
0059   }
0060
0061   //This should be no necessary. TO BE CHECKED
0062   //sort(theFrontPetals.begin(), theFrontPetals.end(), PetalLessPhi());
0063   //sort(theBackPetals.begin(), theBackPetals.end(), PetalLessPhi());
0064
0065   // building disk for front and back petals
0066   setSurface(computeDisk(theComps));
0067   theFrontDisk = computeDisk(theFrontComps);
0068   theBackDisk = computeDisk(theBackComps);
0069
0070   // set up the bin finders
0071   theFrontBinFinder = BinFinderPhi(theFrontComps.front()->position().phi(), theFrontComps.size());
0072   theBackBinFinder = BinFinderPhi(theBackComps.front()->position().phi(), theBackComps.size());
0073
0074   //--------- DEBUG INFO --------------
0075   LogDebug("TkDetLayers") << "DEBUG INFO for TECLayer"
0076                           << "\n"
0077                           << "TECLayer z,perp, innerRadius, outerR: " << this->position().z() << " , "
0078                           << this->position().perp() << " , " << this->specificSurface().innerRadius() << " , "
0079                           << this->specificSurface().outerRadius();
0080
0081   for (auto it = theFrontComps.begin(); it != theFrontComps.end(); it++) {
0082     LogDebug("TkDetLayers") << "frontPetal phi,z,r: " << (*it)->surface().position().phi() << " , "
0083                             << (*it)->surface().position().z() << " , " << (*it)->surface().position().perp();
0084   }
0085
0086   for (auto it = theBackComps.begin(); it != theBackComps.end(); it++) {
0087     LogDebug("TkDetLayers") << "backPetal phi,z,r: " << (*it)->surface().position().phi() << " , "
0088                             << (*it)->surface().position().z() << " , " << (*it)->surface().position().perp();
0089   }
0090   //-----------------------------------
0091 }
0092
0093 TECLayer::~TECLayer() {
0094   for (auto i = theComps.begin(); i != theComps.end(); i++) {
0095     delete *i;
0096   }
0097 }
0098
0099 void TECLayer::groupedCompatibleDetsV(const TrajectoryStateOnSurface& tsos,
0100                                       const Propagator& prop,
0101                                       const MeasurementEstimator& est,
0102                                       std::vector<DetGroup>& result) const {
0103   SubLayerCrossings crossings;
0104   crossings = computeCrossings(tsos, prop.propagationDirection());
0105   if (!crossings.isValid())
0106     return;
0107
0108   vector<DetGroup> closestResult;
0109   addClosest(tsos, prop, est, crossings.closest(), closestResult);
0110   LogDebug("TkDetLayers") << "in TECLayer, closestResult.size(): " << closestResult.size();
0111
0112   // this differs from other groupedCompatibleDets logic, which DON'T check next in such cases!!!
0113   if (closestResult.empty()) {
0114     vector<DetGroup> nextResult;
0115     addClosest(tsos, prop, est, crossings.other(), nextResult);
0116     LogDebug("TkDetLayers") << "in TECLayer, nextResult.size(): " << nextResult.size();
0117     if (nextResult.empty())
0118       return;
0119
0120     DetGroupElement nextGel(nextResult.front().front());
0121     int crossingSide = LayerCrossingSide::endcapSide(nextGel.trajectoryState(), prop);
0122     DetGroupMerger::orderAndMergeTwoLevels(
0123         std::move(closestResult), std::move(nextResult), result, crossings.closestIndex(), crossingSide);
0124   } else {
0125     DetGroupElement closestGel(closestResult.front().front());
0126     float phiWindow = tkDetUtil::computeWindowSize(closestGel.det(), closestGel.trajectoryState(), est);
0127     searchNeighbors(tsos, prop, est, crossings.closest(), phiWindow, closestResult, false);
0128     vector<DetGroup> nextResult;
0129     searchNeighbors(tsos, prop, est, crossings.other(), phiWindow, nextResult, true);
0130
0131     int crossingSide = LayerCrossingSide::endcapSide(closestGel.trajectoryState(), prop);
0132     DetGroupMerger::orderAndMergeTwoLevels(
0133         std::move(closestResult), std::move(nextResult), result, crossings.closestIndex(), crossingSide);
0134   }
0135 }
0136
0137 SubLayerCrossings TECLayer::computeCrossings(const TrajectoryStateOnSurface& startingState,
0138                                              PropagationDirection propDir) const {
0139   double rho(startingState.transverseCurvature());
0140
0141   HelixPlaneCrossing::PositionType startPos(startingState.globalPosition());
0142   HelixPlaneCrossing::DirectionType startDir(startingState.globalMomentum());
0143   HelixForwardPlaneCrossing crossing(startPos, startDir, rho, propDir);
0144
0145   pair<bool, double> frontPath = crossing.pathLength(*theFrontDisk);
0146   if (!frontPath.first)
0147     return SubLayerCrossings();
0148
0149   GlobalPoint gFrontPoint(crossing.position(frontPath.second));
0150
0151   LogDebug("TkDetLayers") << "in TECLayer,front crossing point: r,z,phi: (" << gFrontPoint.perp() << ","
0152                           << gFrontPoint.z() << "," << gFrontPoint.phi() << ")" << endl;
0153
0154   int frontIndex = theFrontBinFinder.binIndex(gFrontPoint.barePhi());
0155   SubLayerCrossing frontSLC(0, frontIndex, gFrontPoint);
0156
0157   pair<bool, double> backPath = crossing.pathLength(*theBackDisk);
0158   if (!backPath.first)
0159     return SubLayerCrossings();
0160
0161   GlobalPoint gBackPoint(crossing.position(backPath.second));
0162   LogDebug("TkDetLayers") << "in TECLayer,back crossing point: r,z,phi: (" << gBackPoint.perp() << ","
0163                           << gFrontPoint.z() << "," << gBackPoint.phi() << ")" << endl;
0164
0165   int backIndex = theBackBinFinder.binIndex(gBackPoint.barePhi());
0166   SubLayerCrossing backSLC(1, backIndex, gBackPoint);
0167
0168   // 0ss: frontDisk has index=0, backDisk has index=1
0169   float frontDist = std::abs(Geom::deltaPhi(gFrontPoint.barePhi(), theFrontComps[frontIndex]->surface().phi()));
0170   float backDist = std::abs(Geom::deltaPhi(gBackPoint.barePhi(), theBackComps[backIndex]->surface().phi()));
0171
0172   if (frontDist < backDist) {
0173     return SubLayerCrossings(frontSLC, backSLC, 0);
0174   } else {
0175     return SubLayerCrossings(backSLC, frontSLC, 1);
0176   }
0177 }
0178
0179 bool TECLayer::addClosest(const TrajectoryStateOnSurface& tsos,
0180                           const Propagator& prop,
0181                           const MeasurementEstimator& est,
0182                           const SubLayerCrossing& crossing,
0183                           vector<DetGroup>& result) const {
0184   const auto& sub(subLayer(crossing.subLayerIndex()));
0185   const auto det(sub[crossing.closestDetIndex()]);
0186
0187   LogDebug("TkDetLayers") << "in TECLayer, adding petal at r,z,phi: (" << det->position().perp() << ","
0188                           << det->position().z() << "," << det->position().phi() << ")" << endl;
0189
0190   return CompatibleDetToGroupAdder().add(*det, tsos, prop, est, result);
0191 }
0192
0193 namespace {
0194   inline bool overlap(float phi, const TECPetal& gsdet, float phiWin) {
0195     const BoundDiskSector& diskSector = gsdet.specificSurface();
0196     pair<float, float> phiRange(phi - phiWin, phi + phiWin);
0197     pair<float, float> petalPhiRange(diskSector.phi() - diskSector.phiHalfExtension(),
0198                                      diskSector.phi() + diskSector.phiHalfExtension());
0199
0200     return rangesIntersect(phiRange, petalPhiRange, [](auto x, auto y) { return Geom::phiLess(x, y); });
0201   }
0202
0203 }  // namespace
0204
0205 void TECLayer::searchNeighbors(const TrajectoryStateOnSurface& tsos,
0206                                const Propagator& prop,
0207                                const MeasurementEstimator& est,
0208                                const SubLayerCrossing& crossing,
0209                                float window,
0210                                vector<DetGroup>& result,
0211                                bool checkClosest) const {
0212   const GlobalPoint& gCrossingPos = crossing.position();
0213   auto gphi = gCrossingPos.barePhi();
0214
0215   const auto& sLayer(subLayer(crossing.subLayerIndex()));
0216
0217   int closestIndex = crossing.closestDetIndex();
0218   int negStartIndex = closestIndex - 1;
0219   int posStartIndex = closestIndex + 1;
0220
0221   if (checkClosest) {  // must decide if the closest is on the neg or pos side
0222     if (Geom::phiLess(gphi, sLayer[closestIndex]->surface().phi())) {
0223       posStartIndex = closestIndex;
0224     } else {
0225       negStartIndex = closestIndex;
0226     }
0227   }
0228
0229   const BinFinderPhi& binFinder = (crossing.subLayerIndex() == 0 ? theFrontBinFinder : theBackBinFinder);
0230
0231   typedef CompatibleDetToGroupAdder Adder;
0232   int half = sLayer.size() / 2;  // to check if dets are called twice....
0233   for (int idet = negStartIndex; idet >= negStartIndex - half; idet--) {
0234     const auto& neighborPetal = *sLayer[binFinder.binIndex(idet)];
0235     if (!overlap(gphi, neighborPetal, window))
0236       break;
0237     if (!Adder::add(neighborPetal, tsos, prop, est, result))
0238       break;
0239     // maybe also add shallow crossing angle test here???
0240   }
0241   for (int idet = posStartIndex; idet < posStartIndex + half; idet++) {
0242     const auto& neighborPetal = *sLayer[binFinder.binIndex(idet)];
0243     if (!overlap(gphi, neighborPetal, window))
0244       break;
0245     if (!Adder::add(neighborPetal, tsos, prop, est, result))
0246       break;
0247     // maybe also add shallow crossing angle test here???
0248   }
0249 }