Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:31:36 

0001 #include "CaloDetIdAssociator.h"
0002 #include "DataFormats/GeometrySurface/interface/Plane.h"
0003 bool CaloDetIdAssociator::crossedElement(const GlobalPoint& point1,
0004                                          const GlobalPoint& point2,
0005                                          const DetId& id,
0006                                          const double toleranceInSigmas,
0007                                          const SteppingHelixStateInfo* initialState) const {
0008   std::vector<GlobalPoint> pointBuffer;
0009   const std::pair<const_iterator, const_iterator>& points = getDetIdPoints(id, pointBuffer);
0010   // fast check
0011   bool xLess(false), xIn(false), xMore(false);
0012   bool yLess(false), yIn(false), yMore(false);
0013   bool zLess(false), zIn(false), zMore(false);
0014   double xMin(point1.x()), xMax(point2.x());
0015   double yMin(point1.y()), yMax(point2.y());
0016   double zMin(point1.z()), zMax(point2.z());
0017   if (xMin > xMax)
0018     std::swap(xMin, xMax);
0019   if (yMin > yMax)
0020     std::swap(yMin, yMax);
0021   if (zMin > zMax)
0022     std::swap(zMin, zMax);
0023   for (std::vector<GlobalPoint>::const_iterator it = points.first; it != points.second; ++it) {
0024     if (it->x() < xMin) {
0025       xLess = true;
0026     } else {
0027       if (it->x() > xMax)
0028         xMore = true;
0029       else
0030         xIn = true;
0031     }
0032     if (it->y() < yMin) {
0033       yLess = true;
0034     } else {
0035       if (it->y() > yMax)
0036         yMore = true;
0037       else
0038         yIn = true;
0039     }
0040     if (it->z() < zMin) {
0041       zLess = true;
0042     } else {
0043       if (it->z() > zMax)
0044         zMore = true;
0045       else
0046         zIn = true;
0047     }
0048   }
0049   if (((xLess && !xIn && !xMore) || (!xLess && !xIn && xMore)) ||
0050       ((yLess && !yIn && !yMore) || (!yLess && !yIn && yMore)) ||
0051       ((zLess && !zIn && !zMore) || (!zLess && !zIn && zMore)))
0052     return false;
0053 
0054   // Define plane normal to the trajectory direction at the first point
0055   GlobalVector vector = (point2 - point1).unit();
0056   float r21 = 0;
0057   float r22 = vector.z() / sqrt(1 - pow(vector.x(), 2));
0058   float r23 = -vector.y() / sqrt(1 - pow(vector.x(), 2));
0059   float r31 = vector.x();
0060   float r32 = vector.y();
0061   float r33 = vector.z();
0062   float r11 = r22 * r33 - r23 * r32;
0063   float r12 = r23 * r31;
0064   float r13 = -r22 * r31;
0065 
0066   Surface::RotationType rotation(r11, r12, r13, r21, r22, r23, r31, r32, r33);
0067   Plane::PlanePointer plane = Plane::build(point1, rotation);
0068   double absoluteTolerance = -1;
0069   if (toleranceInSigmas > 0 && initialState) {
0070     TrajectoryStateOnSurface tsos = initialState->getStateOnSurface(*plane);
0071     if (tsos.isValid() and tsos.hasError()) {
0072       LocalError localErr = tsos.localError().positionError();
0073       localErr.scale(toleranceInSigmas);
0074       float xx = localErr.xx();
0075       float xy = localErr.xy();
0076       float yy = localErr.yy();
0077 
0078       float denom = yy - xx;
0079       float phi = 0.;
0080       if (xy == 0 && denom == 0)
0081         phi = M_PI_4;
0082       else
0083         phi = 0.5 * atan2(2. * xy, denom);  // angle of MAJOR axis
0084       // Unrotate the error ellipse to get the semimajor and minor axes. Then place points on
0085       // the endpoints of semiminor an seminajor axes on original(rotated) error ellipse.
0086       LocalError rotErr = localErr.rotate(-phi);  // xy covariance of rotErr should be zero
0087       float semi1 = sqrt(rotErr.xx());
0088       float semi2 = sqrt(rotErr.yy());
0089       absoluteTolerance = std::max(semi1, semi2);
0090     }
0091   }
0092 
0093   // distance between the points.
0094   double trajectorySegmentLength = (point2 - point1).mag();
0095 
0096   // we need to find the angle that covers all points.
0097   // if it's bigger than 180 degree, we are inside
0098   // otherwise we are outside, i.e. the volume is not crossed
0099   bool allBehind = true;
0100   bool allTooFar = true;
0101   std::vector<GlobalPoint>::const_iterator p = points.first;
0102   if (p == points.second) {
0103     edm::LogWarning("TrackAssociator") << "calo geometry for element " << id.rawId() << "is empty. Ignored";
0104     return false;
0105   }
0106   LocalPoint localPoint = plane->toLocal(*p);
0107   double minPhi = localPoint.phi();
0108   double maxPhi = localPoint.phi();
0109   if (localPoint.z() < 0)
0110     allTooFar = false;
0111   else {
0112     allBehind = false;
0113     if (localPoint.z() < trajectorySegmentLength)
0114       allTooFar = false;
0115   }
0116   ++p;
0117   for (; p != points.second; ++p) {
0118     localPoint = plane->toLocal(*p);
0119     double localPhi = localPoint.phi();
0120     if (localPoint.z() < 0)
0121       allTooFar = false;
0122     else {
0123       allBehind = false;
0124       if (localPoint.z() < trajectorySegmentLength)
0125         allTooFar = false;
0126     }
0127     if (localPhi >= minPhi && localPhi <= maxPhi)
0128       continue;
0129     if (localPhi + 2 * M_PI >= minPhi && localPhi + 2 * M_PI <= maxPhi)
0130       continue;
0131     if (localPhi - 2 * M_PI >= minPhi && localPhi - 2 * M_PI <= maxPhi)
0132       continue;
0133     // find the closest limit
0134     if (localPhi > maxPhi) {
0135       double delta1 = fabs(localPhi - maxPhi);
0136       double delta2 = fabs(localPhi - 2 * M_PI - minPhi);
0137       if (delta1 < delta2)
0138         maxPhi = localPhi;
0139       else
0140         minPhi = localPhi - 2 * M_PI;
0141       continue;
0142     }
0143     if (localPhi < minPhi) {
0144       double delta1 = fabs(localPhi - minPhi);
0145       double delta2 = fabs(localPhi + 2 * M_PI - maxPhi);
0146       if (delta1 < delta2)
0147         minPhi = localPhi;
0148       else
0149         maxPhi = localPhi + 2 * M_PI;
0150       continue;
0151     }
0152     cms::Exception("FatalError")
0153         << "Algorithm logic error - this should never happen. Problems with trajectory-volume matching.";
0154   }
0155   if (allBehind)
0156     return false;
0157   if (allTooFar)
0158     return false;
0159   if (fabs(maxPhi - minPhi) > M_PI)
0160     return true;
0161 
0162   // now if the tolerance is positive, check how far we are
0163   // from the closest line segment
0164   if (absoluteTolerance < 0)
0165     return false;
0166   double distanceToClosestLineSegment = 1e9;
0167   std::vector<GlobalPoint>::const_iterator i, j;
0168   for (i = points.first; i != points.second; ++i)
0169     for (j = i + 1; j != points.second; ++j) {
0170       LocalPoint p1(plane->toLocal(*i));
0171       LocalPoint p2(plane->toLocal(*j));
0172       // now we deal with high school level math to get
0173       // the triangle paramaters
0174       double side1squared = p1.perp2();
0175       double side2squared = p2.perp2();
0176       double side3squared = (p2.x() - p1.x()) * (p2.x() - p1.x()) + (p2.y() - p1.y()) * (p2.y() - p1.y());
0177       double area = fabs(p1.x() * p2.y() - p2.x() * p1.y()) / 2;
0178       // all triangle angles must be smaller than 90 degree
0179       // otherwise the projection is out of the line segment
0180       if (side1squared + side2squared > side3squared && side2squared + side3squared > side1squared &&
0181           side1squared + side3squared > side1squared) {
0182         double h(2 * area / sqrt(side3squared));
0183         if (h < distanceToClosestLineSegment)
0184           distanceToClosestLineSegment = h;
0185       } else {
0186         if (sqrt(side1squared) < distanceToClosestLineSegment)
0187           distanceToClosestLineSegment = sqrt(side1squared);
0188         if (sqrt(side2squared) < distanceToClosestLineSegment)
0189           distanceToClosestLineSegment = sqrt(side2squared);
0190       }
0191     }
0192   if (distanceToClosestLineSegment < absoluteTolerance)
0193     return true;
0194   return false;
0195 }
0196 
0197 void CaloDetIdAssociator::check_setup() const {
0198   DetIdAssociator::check_setup();
0199   if (geometry_ == nullptr)
0200     throw cms::Exception("CaloGeometry is not set");
0201 }
0202 
0203 GlobalPoint CaloDetIdAssociator::getPosition(const DetId& id) const {
0204   return geometry_->getSubdetectorGeometry(id)->getGeometry(id)->getPosition();
0205 }
0206 
0207 void CaloDetIdAssociator::getValidDetIds(unsigned int subDectorIndex, std::vector<DetId>& detIds) const {
0208   if (subDectorIndex != 0)
0209     cms::Exception("FatalError")
0210         << "Calo sub-dectors are all handle as one sub-system, but subDetectorIndex is not zero.\n";
0211   detIds = geometry_->getValidDetIds(DetId::Calo, 1);
0212 }
0213 
0214 std::pair<DetIdAssociator::const_iterator, DetIdAssociator::const_iterator> CaloDetIdAssociator::getDetIdPoints(
0215     const DetId& id, std::vector<GlobalPoint>& points) const {
0216   const CaloSubdetectorGeometry* subDetGeom = geometry_->getSubdetectorGeometry(id);
0217   if (!subDetGeom) {
0218     LogDebug("TrackAssociator") << "Cannot find sub-detector geometry for " << id.rawId() << "\n";
0219     return std::pair<const_iterator, const_iterator>(dummy_.end(), dummy_.end());
0220   }
0221   auto cellGeom = subDetGeom->getGeometry(id);
0222   if (!cellGeom) {
0223     LogDebug("TrackAssociator") << "Cannot find CaloCell geometry for " << id.rawId() << "\n";
0224     return std::pair<const_iterator, const_iterator>(dummy_.end(), dummy_.end());
0225   }
0226   const CaloCellGeometry::CornersVec& cor(cellGeom->getCorners());
0227   return std::pair<const_iterator, const_iterator>(cor.begin(), cor.end());
0228 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!