TkDetLayers/src/TkDetUtil.cc

0001 #include "TkDetUtil.h"
0002
0003 #include "TrackingTools/DetLayers/interface/MeasurementEstimator.h"
0004 #include "DataFormats/GeometrySurface/interface/Plane.h"
0005 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0006 #include "DataFormats/GeometryVector/interface/VectorUtil.h"
0007
0008 namespace tkDetUtil {
0009
0010   float computeWindowSize(const GeomDet* det, const TrajectoryStateOnSurface& tsos, const MeasurementEstimator& est) {
0011     const Plane& startPlane = det->surface();
0012     auto maxDistance = est.maximalLocalDisplacement(tsos, startPlane);
0013     return std::copysign(calculatePhiWindow(maxDistance, tsos, startPlane), maxDistance.x());
0014   }
0015
0016   float calculatePhiWindow(const MeasurementEstimator::Local2DVector& imaxDistance,
0017                            const TrajectoryStateOnSurface& ts,
0018                            const Plane& plane) {
0019     MeasurementEstimator::Local2DVector maxDistance(std::abs(imaxDistance.x()), std::abs(imaxDistance.y()));
0020
0021     constexpr float tolerance = 1.e-6;
0022     LocalPoint start = ts.localPosition();
0023     //     std::cout << "plane z " << plane.normalVector() << std::endl;
0024     float dphi = 0;
0025     if LIKELY (std::abs(1.f - std::abs(plane.normalVector().z())) < tolerance) {
0026       auto ori = plane.toLocal(GlobalPoint(0., 0., 0.));
0027       auto xc = std::abs(start.x() - ori.x());
0028       auto yc = std::abs(start.y() - ori.y());
0029
0030       if (yc < maxDistance.y() && xc < maxDistance.x())
0031         return M_PI;
0032
0033       auto hori = yc > maxDistance.y();  // quadrant 1 (&2), otherwiase quadrant 1&4
0034       auto y0 = hori ? yc + std::copysign(maxDistance.y(), xc - maxDistance.x()) : xc - maxDistance.x();
0035       auto x0 = hori ? xc - maxDistance.x() : -yc - maxDistance.y();
0036       auto y1 = hori ? yc - maxDistance.y() : xc - maxDistance.x();
0037       auto x1 = hori ? xc + maxDistance.x() : -yc + maxDistance.y();
0038
0039       auto sp = (x0 * x1 + y0 * y1) / std::sqrt((x0 * x0 + y0 * y0) * (x1 * x1 + y1 * y1));
0040       sp = std::min(std::max(sp, -1.f), 1.f);
0041       dphi = std::acos(sp);
0042
0043       return dphi;
0044     }
0045
0046     // generic algo
0047     float corners[] = {plane.toGlobal(LocalPoint(start.x() + maxDistance.x(), start.y() + maxDistance.y())).barePhi(),
0048                        plane.toGlobal(LocalPoint(start.x() - maxDistance.x(), start.y() + maxDistance.y())).barePhi(),
0049                        plane.toGlobal(LocalPoint(start.x() - maxDistance.x(), start.y() - maxDistance.y())).barePhi(),
0050                        plane.toGlobal(LocalPoint(start.x() + maxDistance.x(), start.y() - maxDistance.y())).barePhi()};
0051
0052     float phimin = corners[0];
0053     float phimax = phimin;
0054     for (int i = 1; i < 4; i++) {
0055       float cPhi = corners[i];
0056       if (Geom::phiLess(cPhi, phimin)) {
0057         phimin = cPhi;
0058       }
0059       if (Geom::phiLess(phimax, cPhi)) {
0060         phimax = cPhi;
0061       }
0062     }
0063     float phiWindow = phimax - phimin;
0064     if (phiWindow < 0.) {
0065       phiWindow += 2. * Geom::pi();
0066     }
0067     // std::cout << "phiWindow " << phiWindow << ' ' << dphi << ' ' << dphi-phiWindow  << std::endl;
0068     return phiWindow;
0069   }
0070
0071   float computeYdirWindowSize(const GeomDet* det,
0072                               const TrajectoryStateOnSurface& tsos,
0073                               const MeasurementEstimator& est) {
0074     const Plane& startPlane = det->surface();
0075     MeasurementEstimator::Local2DVector maxDistance = est.maximalLocalDisplacement(tsos, startPlane);
0076     return maxDistance.y();
0077   }
0078
0079   std::array<int, 3> findThreeClosest(const std::vector<RingPar>& ringParams,
0080                                       const std::vector<GlobalPoint>& ringCrossing,
0081                                       const int ringSize) {
0082     std::array<int, 3> theBins = {{-1, -1, -1}};
0083     theBins[0] = 0;
0084     float initialR = ringParams[0].theRingR;
0085     float rDiff0 = std::abs(ringCrossing[0].perp() - initialR);
0086     float rDiff1 = -1.;
0087     float rDiff2 = -1.;
0088     for (int i = 1; i < ringSize; i++) {
0089       float ringR = ringParams[i].theRingR;
0090       float testDiff = std::abs(ringCrossing[i].perp() - ringR);
0091       if (testDiff < rDiff0) {
0092         rDiff2 = rDiff1;
0093         rDiff1 = rDiff0;
0094         rDiff0 = testDiff;
0095         theBins[2] = theBins[1];
0096         theBins[1] = theBins[0];
0097         theBins[0] = i;
0098       } else if (rDiff1 < 0 || testDiff < rDiff1) {
0099         rDiff2 = rDiff1;
0100         rDiff1 = testDiff;
0101         theBins[2] = theBins[1];
0102         theBins[1] = i;
0103       } else if (rDiff2 < 0 || testDiff < rDiff2) {
0104         rDiff2 = testDiff;
0105         theBins[2] = i;
0106       }
0107     }
0108
0109     return theBins;
0110   }
0111
0112   bool overlapInR(const TrajectoryStateOnSurface& tsos, int index, double ymax, const std::vector<RingPar>& ringParams) {
0113     // assume "fixed theta window", i.e. margin in local y = r is changing linearly with z
0114     float tsRadius = tsos.globalPosition().perp();
0115     float thetamin =
0116         (std::max(0., tsRadius - ymax)) / (std::abs(tsos.globalPosition().z()) + 10.f);  // add 10 cm contingency
0117     float thetamax = (tsRadius + ymax) / (std::abs(tsos.globalPosition().z()) - 10.f);
0118
0119     // do the theta regions overlap ?
0120
0121     return !(thetamin > ringParams[index].thetaRingMax || ringParams[index].thetaRingMin > thetamax);
0122   }
0123
0124   RingPar fillRingParametersFromDisk(const BoundDisk& ringDisk) {
0125     float ringMinZ = std::abs(ringDisk.position().z()) - ringDisk.bounds().thickness() / 2.;
0126     float ringMaxZ = std::abs(ringDisk.position().z()) + ringDisk.bounds().thickness() / 2.;
0127     RingPar tempPar;
0128     tempPar.thetaRingMin = ringDisk.innerRadius() / ringMaxZ;
0129     tempPar.thetaRingMax = ringDisk.outerRadius() / ringMinZ;
0130     tempPar.theRingR = (ringDisk.innerRadius() + ringDisk.outerRadius()) / 2.;
0131     return tempPar;
0132   }
0133
0134 }  // namespace tkDetUtil