tkDetUtil Namespace Reference

Classes
struct	RingPar

Functions
float	calculatePhiWindow (const MeasurementEstimator::Local2DVector &imaxDistance, const TrajectoryStateOnSurface &ts, const Plane &plane)
template<class T >
BoundDisk *	computeDisk (const std::vector< const T *> &structures)
float	computeWindowSize (const GeomDet *det, const TrajectoryStateOnSurface &tsos, const MeasurementEstimator &est)
float	computeYdirWindowSize (const GeomDet *det, const TrajectoryStateOnSurface &tsos, const MeasurementEstimator &est)
RingPar	fillRingParametersFromDisk (const BoundDisk &ringDisk)
std::array< int, 3 >	findThreeClosest (const std::vector< RingPar > &ringParams, const std::vector< GlobalPoint > &ringCrossing, const int ringSize)
template<class T >
void	groupedCompatibleDetsV (const TrajectoryStateOnSurface &startingState, const Propagator &prop, const MeasurementEstimator &est, std::vector< DetGroup > &result, const int ringSize, const std::vector< const T *> &diskComponents, const std::vector< RingPar > &ringParams)
bool	overlapInPhi (float phi, const GeomDet &det, float phiWindow)
bool	overlapInPhi (GlobalPoint crossPoint, const GeomDet &det, float phiWindow)
bool	overlapInR (const TrajectoryStateOnSurface &tsos, int index, double ymax, const std::vector< RingPar > &ringParams)
template<class T >
std::array< int, 3 >	ringIndicesByCrossingProximity (const TrajectoryStateOnSurface &startingState, const Propagator &prop, const int ringSize, const T &diskComponents, const std::vector< RingPar > &ringParams)

Function Documentation

calculatePhiWindow()

float tkDetUtil::calculatePhiWindow	(	const MeasurementEstimator::Local2DVector &	imaxDistance,
		const TrajectoryStateOnSurface &	ts,
		const Plane &	plane
	)

Definition at line 16 of file TkDetUtil.cc.

References funct::abs(), PV3DBase< T, PVType, FrameType >::barePhi(), ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr(), f, mps_fire::i, LIKELY, TrajectoryStateOnSurface::localPosition(), M_PI, SiStripPI::max, particleFlowClusterHGC_cfi::maxDistance, SiStripPI::min, Plane::normalVector(), Geom::phiLess(), phimax, phimin, trackingPOGFilters_cfi::phiWindow, Geom::pi(), mathSSE::sqrt(), Surface::toGlobal(), tolerance, GloballyPositioned< T >::toLocal(), PV2DBase< T, PVType, FrameType >::x(), PV2DBase< T, PVType, FrameType >::y(), testProducerWithPsetDescEmpty_cfi::y1, and PV3DBase< T, PVType, FrameType >::z().

Referenced by computeWindowSize().

                                               {
    MeasurementEstimator::Local2DVector maxDistance(std::abs(imaxDistance.x()), std::abs(imaxDistance.y()));

    constexpr float tolerance = 1.e-6;
    LocalPoint start = ts.localPosition();
    //     std::cout << "plane z " << plane.normalVector() << std::endl;
    float dphi = 0;
    if LIKELY (std::abs(1.f - std::abs(plane.normalVector().z())) < tolerance) {
      auto ori = plane.toLocal(GlobalPoint(0., 0., 0.));
      auto xc = std::abs(start.x() - ori.x());
      auto yc = std::abs(start.y() - ori.y());

      if (yc < maxDistance.y() && xc < maxDistance.x())
        return M_PI;

      auto hori = yc > maxDistance.y();  // quadrant 1 (&2), otherwiase quadrant 1&4
      auto y0 = hori ? yc + std::copysign(maxDistance.y(), xc - maxDistance.x()) : xc - maxDistance.x();
      auto x0 = hori ? xc - maxDistance.x() : -yc - maxDistance.y();
      auto y1 = hori ? yc - maxDistance.y() : xc - maxDistance.x();
      auto x1 = hori ? xc + maxDistance.x() : -yc + maxDistance.y();

      auto sp = (x0 * x1 + y0 * y1) / std::sqrt((x0 * x0 + y0 * y0) * (x1 * x1 + y1 * y1));
      sp = std::min(std::max(sp, -1.f), 1.f);
      dphi = std::acos(sp);

      return dphi;
    }

    // generic algo
    float corners[] = {plane.toGlobal(LocalPoint(start.x() + maxDistance.x(), start.y() + maxDistance.y())).barePhi(),
                       plane.toGlobal(LocalPoint(start.x() - maxDistance.x(), start.y() + maxDistance.y())).barePhi(),
                       plane.toGlobal(LocalPoint(start.x() - maxDistance.x(), start.y() - maxDistance.y())).barePhi(),
                       plane.toGlobal(LocalPoint(start.x() + maxDistance.x(), start.y() - maxDistance.y())).barePhi()};

    float phimin = corners[0];
    float phimax = phimin;
    for (int i = 1; i < 4; i++) {
      float cPhi = corners[i];
      if (Geom::phiLess(cPhi, phimin)) {
        phimin = cPhi;
      }
      if (Geom::phiLess(phimax, cPhi)) {
        phimax = cPhi;
      }
    }
    float phiWindow = phimax - phimin;
    if (phiWindow < 0.) {
      phiWindow += 2. * Geom::pi();
    }
    // std::cout << "phiWindow " << phiWindow << ' ' << dphi << ' ' << dphi-phiWindow  << std::endl;
    return phiWindow;
  }

computeDisk()

template<class T >

BoundDisk* tkDetUtil::computeDisk

(

const std::vector< const T *> &

structures

)

Definition at line 236 of file TkDetUtil.h.

References mps_fire::i, SiStripPI::max, SiStripPI::min, makeMuonMisalignmentScenario::rot, SiStripMonitorCluster_cfi::zmax, and SiStripMonitorCluster_cfi::zmin.

Referenced by Phase2EndcapLayerDoubleDisk::computeDisk(), Phase2EndcapSubDisk::computeDisk(), Phase2EndcapLayer::computeDisk(), and TECLayer::TECLayer().

                                                                {
    float theRmin = structures.front()->specificSurface().innerRadius();
    float theRmax = structures.front()->specificSurface().outerRadius();
    float theZmin = structures.front()->position().z() - structures.front()->surface().bounds().thickness() / 2;
    float theZmax = structures.front()->position().z() + structures.front()->surface().bounds().thickness() / 2;

    for (typename std::vector<const T*>::const_iterator i = structures.begin(); i != structures.end(); i++) {
      float rmin = (**i).specificSurface().innerRadius();
      float rmax = (**i).specificSurface().outerRadius();
      float zmin = (**i).position().z() - (**i).surface().bounds().thickness() / 2.;
      float zmax = (**i).position().z() + (**i).surface().bounds().thickness() / 2.;
      theRmin = std::min(theRmin, rmin);
      theRmax = std::max(theRmax, rmax);
      theZmin = std::min(theZmin, zmin);
      theZmax = std::max(theZmax, zmax);
    }

    float zPos = (theZmax + theZmin) / 2.;
    Plane::PositionType pos(0., 0., zPos);
    Plane::RotationType rot;

    return new BoundDisk(pos, rot, new SimpleDiskBounds(theRmin, theRmax, theZmin - zPos, theZmax - zPos));
  }

computeWindowSize()

float tkDetUtil::computeWindowSize	(	const GeomDet *	det,
		const TrajectoryStateOnSurface &	tsos,
		const MeasurementEstimator &	est
	)

Definition at line 10 of file TkDetUtil.cc.

References calculatePhiWindow(), particleFlowClusterHGC_cfi::maxDistance, MeasurementEstimator::maximalLocalDisplacement(), and GeomDet::surface().

Referenced by TECLayer::groupedCompatibleDetsV(), TIDRing::groupedCompatibleDetsV(), CompositeTECWedge::groupedCompatibleDetsV(), Phase2EndcapSingleRing::groupedCompatibleDetsV(), and Phase2EndcapRing::groupedCompatibleDetsV().

                                                                                                                     {
    const Plane& startPlane = det->surface();
    auto maxDistance = est.maximalLocalDisplacement(tsos, startPlane);
    return std::copysign(calculatePhiWindow(maxDistance, tsos, startPlane), maxDistance.x());
  }

computeYdirWindowSize()

float tkDetUtil::computeYdirWindowSize	(	const GeomDet *	det,
		const TrajectoryStateOnSurface &	tsos,
		const MeasurementEstimator &	est
	)

Definition at line 71 of file TkDetUtil.cc.

References particleFlowClusterHGC_cfi::maxDistance, MeasurementEstimator::maximalLocalDisplacement(), and GeomDet::surface().

Referenced by Phase2EndcapSubDisk::computeWindowSize(), Phase2EndcapLayer::computeWindowSize(), and groupedCompatibleDetsV().

                                                               {
    const Plane& startPlane = det->surface();
    MeasurementEstimator::Local2DVector maxDistance = est.maximalLocalDisplacement(tsos, startPlane);
    return maxDistance.y();
  }

fillRingParametersFromDisk()

RingPar tkDetUtil::fillRingParametersFromDisk

(

const BoundDisk &

ringDisk

)

Definition at line 124 of file TkDetUtil.cc.

References funct::abs(), tkDetUtil::RingPar::theRingR, tkDetUtil::RingPar::thetaRingMax, and tkDetUtil::RingPar::thetaRingMin.

Referenced by Phase2EndcapSubDisk::fillRingPars(), and Phase2EndcapLayer::fillRingPars().

                                                                {
    float ringMinZ = std::abs(ringDisk.position().z()) - ringDisk.bounds().thickness() / 2.;
    float ringMaxZ = std::abs(ringDisk.position().z()) + ringDisk.bounds().thickness() / 2.;
    RingPar tempPar;
    tempPar.thetaRingMin = ringDisk.innerRadius() / ringMaxZ;
    tempPar.thetaRingMax = ringDisk.outerRadius() / ringMinZ;
    tempPar.theRingR = (ringDisk.innerRadius() + ringDisk.outerRadius()) / 2.;
    return tempPar;
  }

findThreeClosest()

std::array< int, 3 > tkDetUtil::findThreeClosest	(	const std::vector< RingPar > &	ringParams,
		const std::vector< GlobalPoint > &	ringCrossing,
		const int	ringSize
	)

Definition at line 79 of file TkDetUtil.cc.

References funct::abs(), mps_fire::i, perp(), and pfIsolationCalculator_cfi::ringSize.

Referenced by ringIndicesByCrossingProximity().

                                                          {
    std::array<int, 3> theBins = {{-1, -1, -1}};
    theBins[0] = 0;
    float initialR = ringParams[0].theRingR;
    float rDiff0 = std::abs(ringCrossing[0].perp() - initialR);
    float rDiff1 = -1.;
    float rDiff2 = -1.;
    for (int i = 1; i < ringSize; i++) {
      float ringR = ringParams[i].theRingR;
      float testDiff = std::abs(ringCrossing[i].perp() - ringR);
      if (testDiff < rDiff0) {
        rDiff2 = rDiff1;
        rDiff1 = rDiff0;
        rDiff0 = testDiff;
        theBins[2] = theBins[1];
        theBins[1] = theBins[0];
        theBins[0] = i;
      } else if (rDiff1 < 0 || testDiff < rDiff1) {
        rDiff2 = rDiff1;
        rDiff1 = testDiff;
        theBins[2] = theBins[1];
        theBins[1] = i;
      } else if (rDiff2 < 0 || testDiff < rDiff2) {
        rDiff2 = testDiff;
        theBins[2] = i;
      }
    }

    return theBins;
  }

groupedCompatibleDetsV()

template<class T >

void tkDetUtil::groupedCompatibleDetsV	(	const TrajectoryStateOnSurface &	startingState,
		const Propagator &	prop,
		const MeasurementEstimator &	est,
		std::vector< DetGroup > &	result,
		const int	ringSize,
		const std::vector< const T *> &	diskComponents,
		const std::vector< RingPar > &	ringParams
	)

Definition at line 100 of file TkDetUtil.h.

References funct::abs(), DetGroupMerger::addSameLevel(), alongMomentum, computeYdirWindowSize(), ntuplemaker::fill, TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalPosition(), mps_fire::i, eostools::move(), DetGroupMerger::orderAndMergeTwoLevels(), overlapInR(), position, Propagator::propagationDirection(), mps_fire::result, ringIndicesByCrossingProximity(), pfIsolationCalculator_cfi::ringSize, and PV3DBase< T, PVType, FrameType >::z().

Referenced by Phase2EndcapLayer::groupedCompatibleDetsV(), and Phase2EndcapSubDisk::groupedCompatibleDetsV().

                                                                    {
    std::array<int, 3> const& ringIndices =
        ringIndicesByCrossingProximity(startingState, prop, ringSize, diskComponents, ringParams);
    if (ringIndices[0] == -1 || ringIndices[1] == -1 || ringIndices[2] == -1) {
      edm::LogError("TkDetLayers") << "TkRingedForwardLayer::groupedCompatibleDets : error in CrossingProximity";
      return;
    }

    //order is: rings in front = 0; rings in back = 1
    //rings should be already ordered in r
    //if the layer has 1 ring, it does not matter
    //FIXME: to be optimized once the geometry is stable
    std::vector<int> ringOrder(ringSize);
    std::fill(ringOrder.begin(), ringOrder.end(), 1);
    if (ringSize > 1) {
      if (std::abs(diskComponents[0]->position().z()) < std::abs(diskComponents[1]->position().z())) {
        for (int i = 0; i < ringSize; i++) {
          if (i % 2 == 0)
            ringOrder[i] = 0;
        }
      } else if (std::abs(diskComponents[0]->position().z()) > std::abs(diskComponents[1]->position().z())) {
        std::fill(ringOrder.begin(), ringOrder.end(), 0);
        for (int i = 0; i < ringSize; i++) {
          if (i % 2 == 0)
            ringOrder[i] = 1;
        }
      } else {
        throw DetLayerException("Rings in Endcap Layer have same z position, no idea how to order them!");
      }
    }

    auto index = [&ringIndices, &ringOrder](int i) { return ringOrder[ringIndices[i]]; };

    std::vector<DetGroup> closestResult;
    diskComponents[ringIndices[0]]->groupedCompatibleDetsV(startingState, prop, est, closestResult);
    // if the closest is empty, use the next one and exit: inherited from TID !
    if (closestResult.empty()) {
      diskComponents[ringIndices[1]]->groupedCompatibleDetsV(startingState, prop, est, result);
      return;
    }

    DetGroupElement closestGel(closestResult.front().front());
    float rWindow = computeYdirWindowSize(closestGel.det(), closestGel.trajectoryState(), est);

    // check if next ring and next next ring are found and if there is overlap

    bool ring1ok =
        ringIndices[1] != -1 && overlapInR(closestGel.trajectoryState(), ringIndices[1], rWindow, ringParams);
    bool ring2ok =
        ringIndices[2] != -1 && overlapInR(closestGel.trajectoryState(), ringIndices[2], rWindow, ringParams);

    // look for the two rings in the same plane (are they only two?)

    // determine if we are propagating from in to out (0) or from out to in (1)

    int direction = 0;
    if (startingState.globalPosition().z() * startingState.globalMomentum().z() > 0) {
      if (prop.propagationDirection() == alongMomentum)
        direction = 0;
      else
        direction = 1;
    } else {
      if (prop.propagationDirection() == alongMomentum)
        direction = 1;
      else
        direction = 0;
    }

    if ((index(0) == index(1)) && (index(0) == index(2))) {
      edm::LogInfo("AllRingsInOnePlane") << " All rings: " << ringIndices[0] << " " << ringIndices[1] << " "
                                         << ringIndices[2] << " in one plane. Only the first two will be considered";
      ring2ok = false;
    }

    if (index(0) == index(1)) {
      if (ring1ok) {
        std::vector<DetGroup> ring1res;
        diskComponents[ringIndices[1]]->groupedCompatibleDetsV(startingState, prop, est, ring1res);
        DetGroupMerger::addSameLevel(std::move(ring1res), closestResult);
      }
      if (ring2ok) {
        std::vector<DetGroup> ring2res;
        diskComponents[ringIndices[2]]->groupedCompatibleDetsV(startingState, prop, est, ring2res);
        DetGroupMerger::orderAndMergeTwoLevels(
            std::move(closestResult), std::move(ring2res), result, index(0), direction);
        return;
      } else {
        result.swap(closestResult);
        return;
      }
    } else if (index(0) == index(2)) {
      if (ring2ok) {
        std::vector<DetGroup> ring2res;
        diskComponents[ringIndices[2]]->groupedCompatibleDetsV(startingState, prop, est, ring2res);
        DetGroupMerger::addSameLevel(std::move(ring2res), closestResult);
      }
      if (ring1ok) {
        std::vector<DetGroup> ring1res;
        diskComponents[ringIndices[1]]->groupedCompatibleDetsV(startingState, prop, est, ring1res);
        DetGroupMerger::orderAndMergeTwoLevels(
            std::move(closestResult), std::move(ring1res), result, index(0), direction);
        return;
      } else {
        result.swap(closestResult);
        return;
      }
    } else {
      std::vector<DetGroup> ring12res;
      if (ring1ok) {
        std::vector<DetGroup> ring1res;
        diskComponents[ringIndices[1]]->groupedCompatibleDetsV(startingState, prop, est, ring1res);
        ring12res.swap(ring1res);
      }
      if (ring2ok) {
        std::vector<DetGroup> ring2res;
        diskComponents[ringIndices[2]]->groupedCompatibleDetsV(startingState, prop, est, ring2res);
        DetGroupMerger::addSameLevel(std::move(ring2res), ring12res);
      }
      if (!ring12res.empty()) {
        DetGroupMerger::orderAndMergeTwoLevels(
            std::move(closestResult), std::move(ring12res), result, index(0), direction);
        return;
      } else {
        result.swap(closestResult);
        return;
      }
    }
  }

overlapInPhi() [1/2]

bool tkDetUtil::overlapInPhi ( float  phi, const GeomDet &  det, float  phiWindow  )

inline

Definition at line 32 of file TkDetUtil.h.

References Geom::phiLess(), Surface::phiSpan(), trackingPOGFilters_cfi::phiWindow, rangesIntersect(), GeomDet::surface(), and x.

Referenced by overlapInPhi(), CompositeTECWedge::searchNeighbors(), Phase2EndcapSingleRing::searchNeighbors(), TIDRing::searchNeighbors(), and Phase2EndcapRing::searchNeighbors().

                                                                           {
    std::pair<float, float> phiRange(phi - phiWindow, phi + phiWindow);
    return rangesIntersect(phiRange, det.surface().phiSpan(), [](auto x, auto y) { return Geom::phiLess(x, y); });
  }

overlapInPhi() [2/2]

bool tkDetUtil::overlapInPhi ( GlobalPoint  crossPoint, const GeomDet &  det, float  phiWindow  )

inline

Definition at line 37 of file TkDetUtil.h.

References PV3DBase< T, PVType, FrameType >::barePhi(), overlapInPhi(), and trackingPOGFilters_cfi::phiWindow.

                                                                                        {
    return overlapInPhi(crossPoint.barePhi(), det, phiWindow);
  }

overlapInR()

bool tkDetUtil::overlapInR	(	const TrajectoryStateOnSurface &	tsos,
		int	index,
		double	ymax,
		const std::vector< RingPar > &	ringParams
	)

Definition at line 112 of file TkDetUtil.cc.

References funct::abs(), f, TrajectoryStateOnSurface::globalPosition(), SiStripPI::max, PV3DBase< T, PVType, FrameType >::perp(), L1TOccupancyClient_cfi::ymax, and PV3DBase< T, PVType, FrameType >::z().

Referenced by groupedCompatibleDetsV(), Phase2EndcapSubDisk::overlapInR(), and Phase2EndcapLayer::overlapInR().

                                                                                                                      {
    // assume "fixed theta window", i.e. margin in local y = r is changing linearly with z
    float tsRadius = tsos.globalPosition().perp();
    float thetamin =
        (std::max(0., tsRadius - ymax)) / (std::abs(tsos.globalPosition().z()) + 10.f);  // add 10 cm contingency
    float thetamax = (tsRadius + ymax) / (std::abs(tsos.globalPosition().z()) - 10.f);

    // do the theta regions overlap ?

    return !(thetamin > ringParams[index].thetaRingMax || ringParams[index].thetaRingMin > thetamax);
  }

ringIndicesByCrossingProximity()

template<class T >

std::array<int, 3> tkDetUtil::ringIndicesByCrossingProximity	(	const TrajectoryStateOnSurface &	startingState,
		const Propagator &	prop,
		const int	ringSize,
		const T &	diskComponents,
		const std::vector< RingPar > &	ringParams
	)

Definition at line 60 of file TkDetUtil.h.

References findThreeClosest(), TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalPosition(), mps_fire::i, Propagator::propagationDirection(), pfIsolationCalculator_cfi::ringSize, and TrajectoryStateOnSurface::transverseCurvature().

Referenced by groupedCompatibleDetsV().

                                                                                          {
    typedef HelixForwardPlaneCrossing Crossing;
    typedef MeasurementEstimator::Local2DVector Local2DVector;

    HelixPlaneCrossing::PositionType startPos(startingState.globalPosition());
    HelixPlaneCrossing::DirectionType startDir(startingState.globalMomentum());
    PropagationDirection propDir(prop.propagationDirection());
    float rho(startingState.transverseCurvature());

    // calculate the crossings with the ring surfaces
    // rings are assumed to be sorted in R !

    Crossing myXing(startPos, startDir, rho, propDir);

    std::vector<GlobalPoint> ringCrossings;
    ringCrossings.reserve(ringSize);

    for (int i = 0; i < ringSize; i++) {
      const BoundDisk& theRing = static_cast<const BoundDisk&>(diskComponents[i]->surface());
      std::pair<bool, double> pathlen = myXing.pathLength(theRing);
      if (pathlen.first) {
        ringCrossings.push_back(GlobalPoint(myXing.position(pathlen.second)));
      } else {
        // TO FIX.... perhaps there is something smarter to do
        ringCrossings.push_back(GlobalPoint(0., 0., 0.));
      }
    }

    //find three closest rings to the crossing

    std::array<int, 3> closests = findThreeClosest(ringParams, ringCrossings, ringSize);

    return closests;
  }