TkDetLayers/src/TkDetUtil.h

0001 #ifndef TkDetLayers_TkDetUtil_h
0002 #define TkDetLayers_TkDetUtil_h
0003
0004 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0005 #include "TrackingTools/DetLayers/interface/MeasurementEstimator.h"
0006 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0007 #include "TrackingTools/DetLayers/interface/rangesIntersect.h"
0008 #include "DataFormats/GeometryVector/interface/VectorUtil.h"
0009 #include "DataFormats/GeometrySurface/interface/BoundDisk.h"
0010 #include "TrackingTools/DetLayers/interface/RingedForwardLayer.h"
0011 #include "DataFormats/GeometrySurface/interface/SimpleDiskBounds.h"
0012 #include "TrackingTools/DetLayers/interface/DetLayerException.h"
0013
0014 #include "TrackingTools/GeomPropagators/interface/HelixForwardPlaneCrossing.h"
0015
0016 #include "DataFormats/GeometrySurface/interface/SimpleDiskBounds.h"
0017
0018 #include "DetGroupMerger.h"
0019
0020 class GeomDet;
0021 class Plane;
0022 class TrajectoryStateOnSurface;
0023
0024 #pragma GCC visibility push(hidden)
0025
0026 namespace tkDetUtil {
0027
0028   struct RingPar {
0029     float theRingR, thetaRingMin, thetaRingMax;
0030   };
0031
0032   inline bool overlapInPhi(float phi, const GeomDet& det, float phiWindow) {
0033     std::pair<float, float> phiRange(phi - phiWindow, phi + phiWindow);
0034     return rangesIntersect(phiRange, det.surface().phiSpan(), [](auto x, auto y) { return Geom::phiLess(x, y); });
0035   }
0036
0037   inline bool overlapInPhi(GlobalPoint crossPoint, const GeomDet& det, float phiWindow) {
0038     return overlapInPhi(crossPoint.barePhi(), det, phiWindow);
0039   }
0040
0041   float computeWindowSize(const GeomDet* det, const TrajectoryStateOnSurface& tsos, const MeasurementEstimator& est);
0042
0043   float calculatePhiWindow(const MeasurementEstimator::Local2DVector& maxDistance,
0044                            const TrajectoryStateOnSurface& ts,
0045                            const Plane& plane);
0046
0047   float computeYdirWindowSize(const GeomDet* det,
0048                               const TrajectoryStateOnSurface& tsos,
0049                               const MeasurementEstimator& est);
0050
0051   std::array<int, 3> findThreeClosest(const std::vector<RingPar>& ringParams,
0052                                       const std::vector<GlobalPoint>& ringCrossing,
0053                                       const int ringSize);
0054
0055   bool overlapInR(const TrajectoryStateOnSurface& tsos, int index, double ymax, const std::vector<RingPar>& ringParams);
0056
0057   RingPar fillRingParametersFromDisk(const BoundDisk& ringDisk);
0058
0059   template <class T>
0060   std::array<int, 3> ringIndicesByCrossingProximity(const TrajectoryStateOnSurface& startingState,
0061                                                     const Propagator& prop,
0062                                                     const int ringSize,
0063                                                     const T& diskComponents,
0064                                                     const std::vector<RingPar>& ringParams) {
0065     typedef HelixForwardPlaneCrossing Crossing;
0066     typedef MeasurementEstimator::Local2DVector Local2DVector;
0067
0068     HelixPlaneCrossing::PositionType startPos(startingState.globalPosition());
0069     HelixPlaneCrossing::DirectionType startDir(startingState.globalMomentum());
0070     PropagationDirection propDir(prop.propagationDirection());
0071     float rho(startingState.transverseCurvature());
0072
0073     // calculate the crossings with the ring surfaces
0074     // rings are assumed to be sorted in R !
0075
0076     Crossing myXing(startPos, startDir, rho, propDir);
0077
0078     std::vector<GlobalPoint> ringCrossings;
0079     ringCrossings.reserve(ringSize);
0080
0081     for (int i = 0; i < ringSize; i++) {
0082       const BoundDisk& theRing = static_cast<const BoundDisk&>(diskComponents[i]->surface());
0083       std::pair<bool, double> pathlen = myXing.pathLength(theRing);
0084       if (pathlen.first) {
0085         ringCrossings.push_back(GlobalPoint(myXing.position(pathlen.second)));
0086       } else {
0087         // TO FIX.... perhaps there is something smarter to do
0088         ringCrossings.push_back(GlobalPoint(0., 0., 0.));
0089       }
0090     }
0091
0092     //find three closest rings to the crossing
0093
0094     std::array<int, 3> closests = findThreeClosest(ringParams, ringCrossings, ringSize);
0095
0096     return closests;
0097   }
0098
0099   template <class T>
0100   void groupedCompatibleDetsV(const TrajectoryStateOnSurface& startingState,
0101                               const Propagator& prop,
0102                               const MeasurementEstimator& est,
0103                               std::vector<DetGroup>& result,
0104                               const int ringSize,
0105                               const std::vector<const T*>& diskComponents,
0106                               const std::vector<RingPar>& ringParams) {
0107     std::array<int, 3> const& ringIndices =
0108         ringIndicesByCrossingProximity(startingState, prop, ringSize, diskComponents, ringParams);
0109     if (ringIndices[0] == -1 || ringIndices[1] == -1 || ringIndices[2] == -1) {
0110       edm::LogError("TkDetLayers") << "TkRingedForwardLayer::groupedCompatibleDets : error in CrossingProximity";
0111       return;
0112     }
0113
0114     //order is: rings in front = 0; rings in back = 1
0115     //rings should be already ordered in r
0116     //if the layer has 1 ring, it does not matter
0117     //FIXME: to be optimized once the geometry is stable
0118     std::vector<int> ringOrder(ringSize);
0119     std::fill(ringOrder.begin(), ringOrder.end(), 1);
0120     if (ringSize > 1) {
0121       if (std::abs(diskComponents[0]->position().z()) < std::abs(diskComponents[1]->position().z())) {
0122         for (int i = 0; i < ringSize; i++) {
0123           if (i % 2 == 0)
0124             ringOrder[i] = 0;
0125         }
0126       } else if (std::abs(diskComponents[0]->position().z()) > std::abs(diskComponents[1]->position().z())) {
0127         std::fill(ringOrder.begin(), ringOrder.end(), 0);
0128         for (int i = 0; i < ringSize; i++) {
0129           if (i % 2 == 0)
0130             ringOrder[i] = 1;
0131         }
0132       } else {
0133         throw DetLayerException("Rings in Endcap Layer have same z position, no idea how to order them!");
0134       }
0135     }
0136
0137     auto index = [&ringIndices, &ringOrder](int i) { return ringOrder[ringIndices[i]]; };
0138
0139     std::vector<DetGroup> closestResult;
0140     diskComponents[ringIndices[0]]->groupedCompatibleDetsV(startingState, prop, est, closestResult);
0141     // if the closest is empty, use the next one and exit: inherited from TID !
0142     if (closestResult.empty()) {
0143       diskComponents[ringIndices[1]]->groupedCompatibleDetsV(startingState, prop, est, result);
0144       return;
0145     }
0146
0147     DetGroupElement closestGel(closestResult.front().front());
0148     float rWindow = computeYdirWindowSize(closestGel.det(), closestGel.trajectoryState(), est);
0149
0150     // check if next ring and next next ring are found and if there is overlap
0151
0152     bool ring1ok =
0153         ringIndices[1] != -1 && overlapInR(closestGel.trajectoryState(), ringIndices[1], rWindow, ringParams);
0154     bool ring2ok =
0155         ringIndices[2] != -1 && overlapInR(closestGel.trajectoryState(), ringIndices[2], rWindow, ringParams);
0156
0157     // look for the two rings in the same plane (are they only two?)
0158
0159     // determine if we are propagating from in to out (0) or from out to in (1)
0160
0161     int direction = 0;
0162     if (startingState.globalPosition().z() * startingState.globalMomentum().z() > 0) {
0163       if (prop.propagationDirection() == alongMomentum)
0164         direction = 0;
0165       else
0166         direction = 1;
0167     } else {
0168       if (prop.propagationDirection() == alongMomentum)
0169         direction = 1;
0170       else
0171         direction = 0;
0172     }
0173
0174     if ((index(0) == index(1)) && (index(0) == index(2))) {
0175       edm::LogInfo("AllRingsInOnePlane") << " All rings: " << ringIndices[0] << " " << ringIndices[1] << " "
0176                                          << ringIndices[2] << " in one plane. Only the first two will be considered";
0177       ring2ok = false;
0178     }
0179
0180     if (index(0) == index(1)) {
0181       if (ring1ok) {
0182         std::vector<DetGroup> ring1res;
0183         diskComponents[ringIndices[1]]->groupedCompatibleDetsV(startingState, prop, est, ring1res);
0184         DetGroupMerger::addSameLevel(std::move(ring1res), closestResult);
0185       }
0186       if (ring2ok) {
0187         std::vector<DetGroup> ring2res;
0188         diskComponents[ringIndices[2]]->groupedCompatibleDetsV(startingState, prop, est, ring2res);
0189         DetGroupMerger::orderAndMergeTwoLevels(
0190             std::move(closestResult), std::move(ring2res), result, index(0), direction);
0191         return;
0192       } else {
0193         result.swap(closestResult);
0194         return;
0195       }
0196     } else if (index(0) == index(2)) {
0197       if (ring2ok) {
0198         std::vector<DetGroup> ring2res;
0199         diskComponents[ringIndices[2]]->groupedCompatibleDetsV(startingState, prop, est, ring2res);
0200         DetGroupMerger::addSameLevel(std::move(ring2res), closestResult);
0201       }
0202       if (ring1ok) {
0203         std::vector<DetGroup> ring1res;
0204         diskComponents[ringIndices[1]]->groupedCompatibleDetsV(startingState, prop, est, ring1res);
0205         DetGroupMerger::orderAndMergeTwoLevels(
0206             std::move(closestResult), std::move(ring1res), result, index(0), direction);
0207         return;
0208       } else {
0209         result.swap(closestResult);
0210         return;
0211       }
0212     } else {
0213       std::vector<DetGroup> ring12res;
0214       if (ring1ok) {
0215         std::vector<DetGroup> ring1res;
0216         diskComponents[ringIndices[1]]->groupedCompatibleDetsV(startingState, prop, est, ring1res);
0217         ring12res.swap(ring1res);
0218       }
0219       if (ring2ok) {
0220         std::vector<DetGroup> ring2res;
0221         diskComponents[ringIndices[2]]->groupedCompatibleDetsV(startingState, prop, est, ring2res);
0222         DetGroupMerger::addSameLevel(std::move(ring2res), ring12res);
0223       }
0224       if (!ring12res.empty()) {
0225         DetGroupMerger::orderAndMergeTwoLevels(
0226             std::move(closestResult), std::move(ring12res), result, index(0), direction);
0227         return;
0228       } else {
0229         result.swap(closestResult);
0230         return;
0231       }
0232     }
0233   }
0234
0235   template <class T>
0236   BoundDisk* computeDisk(const std::vector<const T*>& structures) {
0237     float theRmin = structures.front()->specificSurface().innerRadius();
0238     float theRmax = structures.front()->specificSurface().outerRadius();
0239     float theZmin = structures.front()->position().z() - structures.front()->surface().bounds().thickness() / 2;
0240     float theZmax = structures.front()->position().z() + structures.front()->surface().bounds().thickness() / 2;
0241
0242     for (typename std::vector<const T*>::const_iterator i = structures.begin(); i != structures.end(); i++) {
0243       float rmin = (**i).specificSurface().innerRadius();
0244       float rmax = (**i).specificSurface().outerRadius();
0245       float zmin = (**i).position().z() - (**i).surface().bounds().thickness() / 2.;
0246       float zmax = (**i).position().z() + (**i).surface().bounds().thickness() / 2.;
0247       theRmin = std::min(theRmin, rmin);
0248       theRmax = std::max(theRmax, rmax);
0249       theZmin = std::min(theZmin, zmin);
0250       theZmax = std::max(theZmax, zmax);
0251     }
0252
0253     float zPos = (theZmax + theZmin) / 2.;
0254     Plane::PositionType pos(0., 0., zPos);
0255     Plane::RotationType rot;
0256
0257     return new BoundDisk(pos, rot, new SimpleDiskBounds(theRmin, theRmax, theZmin - zPos, theZmax - zPos));
0258   }
0259
0260 }  // namespace tkDetUtil
0261
0262 #pragma GCC visibility pop
0263 #endif  // TkDetLayers_TkDetUtil_h