CMS 3D CMS Logo

Public Member Functions | Protected Member Functions | Protected Attributes

CaloDetIdAssociator Class Reference

#include <CaloDetIdAssociator.h>

Inheritance diagram for CaloDetIdAssociator:
DetIdAssociator EcalDetIdAssociator HcalDetIdAssociator HODetIdAssociator PreshowerDetIdAssociator

List of all members.

Public Member Functions

 CaloDetIdAssociator ()
 CaloDetIdAssociator (const int nPhi, const int nEta, const double etaBinSize)
 CaloDetIdAssociator (const edm::ParameterSet &pSet)
virtual const GeomDetgetGeomDet (const DetId &id) const
virtual const char * name () const
virtual void setGeometry (const DetIdAssociatorRecord &iRecord)
virtual void setGeometry (const CaloGeometry *ptr)

Protected Member Functions

virtual void check_setup () const
virtual bool crossedElement (const GlobalPoint &, const GlobalPoint &, const DetId &id, const double tolerance=-1, const SteppingHelixStateInfo *=0) const
virtual std::pair
< const_iterator,
const_iterator
getDetIdPoints (const DetId &id) const
virtual GlobalPoint getPosition (const DetId &id) const
virtual const std::vector
< DetId > & 
getValidDetIds (unsigned int subDetectorIndex) const
virtual bool insideElement (const GlobalPoint &point, const DetId &id) const

Protected Attributes

std::vector< GlobalPointdummy_
const CaloGeometrygeometry_

Detailed Description

Definition at line 33 of file CaloDetIdAssociator.h.


Constructor & Destructor Documentation

CaloDetIdAssociator::CaloDetIdAssociator ( ) [inline]

Definition at line 35 of file CaloDetIdAssociator.h.

:DetIdAssociator(72, 70 ,0.087),geometry_(0){};
CaloDetIdAssociator::CaloDetIdAssociator ( const int  nPhi,
const int  nEta,
const double  etaBinSize 
) [inline]

Definition at line 36 of file CaloDetIdAssociator.h.

     :DetIdAssociator(nPhi, nEta, etaBinSize),geometry_(0){};
CaloDetIdAssociator::CaloDetIdAssociator ( const edm::ParameterSet pSet) [inline]

Definition at line 39 of file CaloDetIdAssociator.h.

     :DetIdAssociator(pSet.getParameter<int>("nPhi"),pSet.getParameter<int>("nEta"),pSet.getParameter<double>("etaBinSize")),geometry_(0){};

Member Function Documentation

void CaloDetIdAssociator::check_setup ( ) const [protected, virtual]

Reimplemented from DetIdAssociator.

Definition at line 194 of file CaloDetIdAssociator.cc.

References Exception, and geometry_.

{
  DetIdAssociator::check_setup();
  if (geometry_==0) throw cms::Exception("CaloGeometry is not set");
}
bool CaloDetIdAssociator::crossedElement ( const GlobalPoint point1,
const GlobalPoint point2,
const DetId id,
const double  tolerance = -1,
const SteppingHelixStateInfo initialState = 0 
) const [protected, virtual]

Reimplemented from DetIdAssociator.

Definition at line 3 of file CaloDetIdAssociator.cc.

References newFWLiteAna::build, Exception, getDetIdPoints(), SteppingHelixStateInfo::getStateOnSurface(), h, TrajectoryStateOnSurface::hasError(), i, TrajectoryStateOnSurface::isValid(), j, TrajectoryStateOnSurface::localError(), M_PI, mag(), max(), AlCaHLTBitMon_ParallelJobs::p, p1, p2, PV3DBase< T, PVType, FrameType >::phi(), phi, LocalTrajectoryError::positionError(), funct::pow(), dttmaxenums::r32, LocalError::rotate(), LocalError::scale(), mathSSE::sqrt(), swap(), csvLumiCalc::unit, PV3DBase< T, PVType, FrameType >::x(), LocalError::xx(), LocalError::xy(), xy(), PV3DBase< T, PVType, FrameType >::y(), LocalError::yy(), and PV3DBase< T, PVType, FrameType >::z().

{
   const std::pair<const_iterator,const_iterator>& points = getDetIdPoints(id);
   // fast check
   bool xLess(false), xIn(false), xMore(false);
   bool yLess(false), yIn(false), yMore(false);
   bool zLess(false), zIn(false), zMore(false);
   double xMin(point1.x()), xMax(point2.x());
   double yMin(point1.y()), yMax(point2.y());
   double zMin(point1.z()), zMax(point2.z());
   if ( xMin>xMax ) std::swap(xMin,xMax);
   if ( yMin>yMax ) std::swap(yMin,yMax);
   if ( zMin>zMax ) std::swap(zMin,zMax);
   for ( std::vector<GlobalPoint>::const_iterator it = points.first;
         it != points.second; ++it ){
     if ( it->x()<xMin ){
       xLess = true;
     } else {
       if ( it->x()>xMax )
         xMore = true;
       else
         xIn = true;
     }
     if ( it->y()<yMin ){
       yLess = true;
     } else {
       if ( it->y()>yMax )
         yMore = true;
       else
         yIn = true;
     }
     if ( it->z()<zMin ){
       zLess = true;
     } else {
       if ( it->z()>zMax )
         zMore = true;
       else
         zIn = true;
     }
   }
   if ( ( (xLess && !xIn && !xMore) || (!xLess && !xIn && xMore) ) ||
        ( (yLess && !yIn && !yMore) || (!yLess && !yIn && yMore) ) ||
        ( (zLess && !zIn && !zMore) || (!zLess && !zIn && zMore) ) ) return false;

   // Define plane normal to the trajectory direction at the first point
   GlobalVector vector = (point2-point1).unit();
   float r21 = 0;
   float r22 = vector.z()/sqrt(1-pow(vector.x(),2));
   float r23 = -vector.y()/sqrt(1-pow(vector.x(),2));
   float r31 = vector.x();
   float r32 = vector.y();
   float r33 = vector.z();
   float r11 = r22*r33-r23*r32;
   float r12 = r23*r31;
   float r13 = -r22*r31;
   
   Surface::RotationType rotation(r11, r12, r13,
                                  r21, r22, r23,
                                  r31, r32, r33);
   Plane::PlanePointer plane = Plane::build(point1, rotation);
   double absoluteTolerance = -1;
   if ( toleranceInSigmas>0 && initialState ){
      TrajectoryStateOnSurface tsos = initialState->getStateOnSurface(*plane);
      if ( tsos.isValid() and tsos.hasError()) {
         LocalError localErr = tsos.localError().positionError();
         localErr.scale(toleranceInSigmas); 
         float xx = localErr.xx();
         float xy = localErr.xy();
         float yy = localErr.yy();

         float denom = yy - xx;
         float phi = 0., phi_temp=0.;
         if(xy == 0 && denom==0) phi = M_PI_4;
         else phi = 0.5 * atan2(2.*xy,denom); // angle of MAJOR axis
         phi_temp = phi;
         // Unrotate the error ellipse to get the semimajor and minor axes. Then place points on
         // the endpoints of semiminor an seminajor axes on original(rotated) error ellipse.
         LocalError rotErr = localErr.rotate(-phi); // xy covariance of rotErr should be zero
         float semi1 = sqrt(rotErr.xx());
         float semi2 = sqrt(rotErr.yy());
         absoluteTolerance = std::max(semi1,semi2);
      }
   }
   
   // distance between the points.
   double trajectorySegmentLength = (point2-point1).mag();

   // we need to find the angle that covers all points. 
   // if it's bigger than 180 degree, we are inside
   // otherwise we are outside, i.e. the volume is not crossed
   bool allBehind = true;
   bool allTooFar = true;
   std::vector<GlobalPoint>::const_iterator p = points.first;
   if ( p == points.second ) {
      edm::LogWarning("TrackAssociator") << "calo geometry for element " << id.rawId() << "is empty. Ignored"; 
      return false; 
   }
   LocalPoint localPoint = plane->toLocal(*p);
   double minPhi = localPoint.phi();
   double maxPhi = localPoint.phi();
   if ( localPoint.z() < 0 ) 
     allTooFar = false;
   else {  
      allBehind = false;
      if ( localPoint.z() < trajectorySegmentLength )  allTooFar = false;
   }
   ++p;
   for (; p!=points.second; ++p){
      localPoint = plane->toLocal(*p);
      double localPhi = localPoint.phi();
      if ( localPoint.z() < 0 ) 
        allTooFar = false;
      else {  
         allBehind = false;
         if ( localPoint.z() < trajectorySegmentLength )  allTooFar = false;
      }
      if ( localPhi >= minPhi && localPhi <= maxPhi ) continue;
      if ( localPhi+2*M_PI >= minPhi && localPhi+2*M_PI <= maxPhi ) continue;
      if ( localPhi-2*M_PI >= minPhi && localPhi-2*M_PI <= maxPhi ) continue;
      // find the closest limit
      if ( localPhi > maxPhi ){
         double delta1 = fabs(localPhi-maxPhi);
         double delta2 = fabs(localPhi-2*M_PI-minPhi);
         if ( delta1 < delta2 )
           maxPhi = localPhi;
         else
           minPhi = localPhi-2*M_PI;
         continue;
      }
      if ( localPhi < minPhi ){
         double delta1 = fabs(localPhi-minPhi);
         double delta2 = fabs(localPhi+2*M_PI-maxPhi);
         if ( delta1 < delta2 )
           minPhi = localPhi;
         else
           maxPhi = localPhi+2*M_PI;
         continue;
      }
      cms::Exception("FatalError") << "Algorithm logic error - this should never happen. Problems with trajectory-volume matching.";
   }
   if ( allBehind ) return false;
   if ( allTooFar ) return false;
   if ( fabs(maxPhi-minPhi)>M_PI ) return true;
   
   // now if the tolerance is positive, check how far we are 
   // from the closest line segment
   if (absoluteTolerance < 0 ) return false;
   double distanceToClosestLineSegment = 1e9;
   std::vector<GlobalPoint>::const_iterator i,j;
   for ( i = points.first; i != points.second; ++i )
     for ( j = i+1; j != points.second; ++j )
       {
          LocalPoint p1(plane->toLocal(*i));
          LocalPoint p2(plane->toLocal(*j));
          // now we deal with high school level math to get
          // the triangle paramaters
          double side1squared = p1.perp2();
          double side2squared = p2.perp2();
          double side3squared = (p2.x()-p1.x())*(p2.x()-p1.x()) + (p2.y()-p1.y())*(p2.y()-p1.y());
          double area = fabs(p1.x()*p2.y()-p2.x()*p1.y())/2;
          // all triangle angles must be smaller than 90 degree
          // otherwise the projection is out of the line segment
          if ( side1squared + side2squared > side3squared &&
               side2squared + side3squared > side1squared &&
               side1squared + side3squared > side1squared )
            {
               double h(2*area/sqrt(side3squared));
               if ( h < distanceToClosestLineSegment ) distanceToClosestLineSegment = h;
            }
          else
            {
               if ( sqrt(side1squared) < distanceToClosestLineSegment ) distanceToClosestLineSegment = sqrt(side1squared);
               if ( sqrt(side2squared) < distanceToClosestLineSegment ) distanceToClosestLineSegment = sqrt(side2squared);
            }
       }
   if ( distanceToClosestLineSegment < absoluteTolerance ) return true;
   return false;
}
std::pair< DetIdAssociator::const_iterator, DetIdAssociator::const_iterator > CaloDetIdAssociator::getDetIdPoints ( const DetId id) const [protected, virtual]

Implements DetIdAssociator.

Definition at line 211 of file CaloDetIdAssociator.cc.

References dummy_, geometry_, CaloCellGeometry::getCorners(), CaloSubdetectorGeometry::getGeometry(), CaloGeometry::getSubdetectorGeometry(), and LogDebug.

Referenced by crossedElement().

{
  const CaloSubdetectorGeometry* subDetGeom = geometry_->getSubdetectorGeometry(id);
  if(! subDetGeom){
    LogDebug("TrackAssociator") << "Cannot find sub-detector geometry for " << id.rawId() <<"\n";
    return std::pair<const_iterator,const_iterator>(dummy_.end(),dummy_.end());
  }
  const CaloCellGeometry* cellGeom = subDetGeom->getGeometry(id);
  if(! cellGeom) {
    LogDebug("TrackAssociator") << "Cannot find CaloCell geometry for " << id.rawId() <<"\n";
    return std::pair<const_iterator,const_iterator>(dummy_.end(),dummy_.end());
  } 
  const CaloCellGeometry::CornersVec& cor (cellGeom->getCorners() ) ; 
  return std::pair<const_iterator,const_iterator>( cor.begin(), cor.end() ) ;
}
virtual const GeomDet* CaloDetIdAssociator::getGeomDet ( const DetId id) const [inline, virtual]

Implements DetIdAssociator.

Definition at line 46 of file CaloDetIdAssociator.h.

{ return 0; };
GlobalPoint CaloDetIdAssociator::getPosition ( const DetId id) const [protected, virtual]
const std::vector< DetId > & CaloDetIdAssociator::getValidDetIds ( unsigned int  subDetectorIndex) const [protected, virtual]

Implements DetIdAssociator.

Reimplemented in EcalDetIdAssociator, HcalDetIdAssociator, HODetIdAssociator, and PreshowerDetIdAssociator.

Definition at line 204 of file CaloDetIdAssociator.cc.

References DetId::Calo, Exception, geometry_, and CaloGeometry::getValidDetIds().

{
  if ( subDectorIndex!=0 ) cms::Exception("FatalError") << "Calo sub-dectors are all handle as one sub-system, but subDetectorIndex is not zero.\n";
  return geometry_->getValidDetIds(DetId::Calo, 1);
}
virtual bool CaloDetIdAssociator::insideElement ( const GlobalPoint point,
const DetId id 
) const [inline, protected, virtual]
virtual const char* CaloDetIdAssociator::name ( void  ) const [inline, virtual]

Implements DetIdAssociator.

Reimplemented in EcalDetIdAssociator, HcalDetIdAssociator, HODetIdAssociator, and PreshowerDetIdAssociator.

Definition at line 48 of file CaloDetIdAssociator.h.

{ return "CaloTowers"; }
virtual void CaloDetIdAssociator::setGeometry ( const CaloGeometry ptr) [inline, virtual]

Definition at line 42 of file CaloDetIdAssociator.h.

References geometry_.

Referenced by setGeometry().

{ geometry_ = ptr; };
void CaloDetIdAssociator::setGeometry ( const DetIdAssociatorRecord iRecord) [virtual]

Member Data Documentation

std::vector<GlobalPoint> CaloDetIdAssociator::dummy_ [protected]

Definition at line 69 of file CaloDetIdAssociator.h.

Referenced by getDetIdPoints().