/data/refman/pasoursint/CMSSW_4_1_8_patch13/src/Geometry/EcalAlgo/src/EcalEndcapGeometry.cc

Go to the documentation of this file.

#include "Geometry/CaloGeometry/interface/CaloGenericDetId.h"
#include "Geometry/EcalAlgo/interface/EcalEndcapGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
#include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
#include "FWCore/Utilities/interface/Exception.h"
#include <CLHEP/Geometry/Point3D.h>
#include <CLHEP/Geometry/Plane3D.h>

EcalEndcapGeometry::EcalEndcapGeometry() :
   _nnmods ( 316 ) ,
   _nncrys ( 25 ) ,
   m_borderMgr ( 0 ),
   m_borderPtrVec ( 0 ),
   m_avgZ ( -1 )
{
}

EcalEndcapGeometry::~EcalEndcapGeometry() 
{
   delete m_borderPtrVec ;
   delete m_borderMgr ;
}

unsigned int
EcalEndcapGeometry::alignmentTransformIndexLocal( const DetId& id )
{
   const CaloGenericDetId gid ( id ) ;

   assert( gid.isEE() ) ;
   unsigned int index ( EEDetId(id).ix()/51 + ( EEDetId(id).zside()<0 ? 0 : 2 ) ) ;

   return index ;
}

DetId 
EcalEndcapGeometry::detIdFromLocalAlignmentIndex( unsigned int iLoc )
{
   return EEDetId( 20 + 50*( iLoc%2 ), 50, 2*( iLoc/2 ) - 1 ) ;
}

unsigned int
EcalEndcapGeometry::alignmentTransformIndexGlobal( const DetId& id )
{
   return (unsigned int)DetId::Ecal - 1 ;
}

void 
EcalEndcapGeometry::initializeParms()
{
  zeP=0.;
  zeN=0.;
  unsigned nP=0;
  unsigned nN=0;
  m_nref = 0 ;

  for( uint32_t i ( 0 ) ; i != cellGeometries().size() ; ++i )
  {
     if( 0 != cellGeometries()[i] )
     {
//      addCrystalToZGridmap(i->first,dynamic_cast<const TruncatedPyramid*>(i->second));
        const double z ( cellGeometries()[i]->getPosition().z() ) ;
        if(z>0.)
        {
           zeP+=z;
           ++nP;
        }
        else
        {
           zeN+=z;
           ++nN;
        }
        const EEDetId myId ( EEDetId::detIdFromDenseIndex(i) ) ;
        const unsigned int ix ( myId.ix() ) ;
        const unsigned int iy ( myId.iy() ) ;
        if( ix > m_nref ) m_nref = ix ;
        if( iy > m_nref ) m_nref = iy ;
     }
  }
  if( 0 < nP ) zeP/=(double)nP;
  if( 0 < nN ) zeN/=(double)nN;

  m_xlo[0] =  999 ;
  m_xhi[0] = -999 ;
  m_ylo[0] =  999 ;
  m_yhi[0] = -999 ;
  m_xlo[1] =  999 ;
  m_xhi[1] = -999 ;
  m_ylo[1] =  999 ;
  m_yhi[1] = -999 ;
  for( uint32_t i ( 0 ) ; i != cellGeometries().size() ; ++i )
  {
     const GlobalPoint p ( cellGeometries()[i]->getPosition()  ) ;
     const double z ( p.z() ) ;
     const double zz ( 0 > z ? zeN : zeP ) ;
     const double x ( p.x()*zz/z ) ;
     const double y ( p.y()*zz/z ) ;

     if( 0 > z && x < m_xlo[0] ) m_xlo[0] = x ;
     if( 0 < z && x < m_xlo[1] ) m_xlo[1] = x ;
     if( 0 > z && y < m_ylo[0] ) m_ylo[0] = y ;
     if( 0 < z && y < m_ylo[1] ) m_ylo[1] = y ;

     if( 0 > z && x > m_xhi[0] ) m_xhi[0] = x ;
     if( 0 < z && x > m_xhi[1] ) m_xhi[1] = x ;
     if( 0 > z && y > m_yhi[0] ) m_yhi[0] = y ;
     if( 0 < z && y > m_yhi[1] ) m_yhi[1] = y ;
  }

  m_xoff[0] = ( m_xhi[0] + m_xlo[0] )/2. ;
  m_xoff[1] = ( m_xhi[1] + m_xlo[1] )/2. ;
  m_yoff[0] = ( m_yhi[0] + m_ylo[0] )/2. ;
  m_yoff[1] = ( m_yhi[1] + m_ylo[1] )/2. ;

  m_del = ( m_xhi[0] - m_xlo[0] + m_xhi[1] - m_xlo[1] +
            m_yhi[0] - m_ylo[0] + m_yhi[1] - m_ylo[1]   ) ;

  m_wref = m_del/(4.*(m_nref-1)) ;

  m_xlo[0] -= m_wref/2 ;
  m_xlo[1] -= m_wref/2 ;
  m_xhi[0] += m_wref/2 ;
  m_xhi[1] += m_wref/2 ;

  m_ylo[0] -= m_wref/2 ;
  m_ylo[1] -= m_wref/2 ;
  m_yhi[0] += m_wref/2 ;
  m_yhi[1] += m_wref/2 ;

  m_del += m_wref ;
/*
  std::cout<<"zeP="<<zeP<<", zeN="<<zeN<<", nP="<<nP<<", nN="<<nN<<std::endl ;

  std::cout<<"xlo[0]="<<m_xlo[0]<<", xlo[1]="<<m_xlo[1]<<", xhi[0]="<<m_xhi[0]<<", xhi[1]="<<m_xhi[1]
           <<"\nylo[0]="<<m_ylo[0]<<", ylo[1]="<<m_ylo[1]<<", yhi[0]="<<m_yhi[0]<<", yhi[1]="<<m_yhi[1]<<std::endl ;

  std::cout<<"xoff[0]="<<m_xoff[0]<<", xoff[1]"<<m_xoff[1]<<", yoff[0]="<<m_yoff[0]<<", yoff[1]"<<m_yoff[1]<<std::endl ;

  std::cout<<"nref="<<m_nref<<", m_wref="<<m_wref<<std::endl ;
*/  
}


unsigned int 
EcalEndcapGeometry::xindex( double x,
                            double z ) const
{
   const double xlo ( 0 > z ? m_xlo[0]  : m_xlo[1]  ) ;
   const int i ( 1 + int( ( x - xlo )/m_wref ) ) ;

   return ( 1 > i ? 1 :
            ( m_nref < (unsigned int) i ? m_nref : (unsigned int) i ) ) ;

}

unsigned int 
EcalEndcapGeometry::yindex( double y,
                            double z  ) const
{
   const double ylo ( 0 > z ? m_ylo[0]  : m_ylo[1]  ) ;
   const int i ( 1 + int( ( y - ylo )/m_wref ) ) ;

   return ( 1 > i ? 1 :
            ( m_nref < (unsigned int) i ? m_nref : (unsigned int) i ) ) ;
}

EEDetId 
EcalEndcapGeometry::gId( float x, 
                         float y, 
                         float z ) const
{
   const double       fac ( fabs( ( 0 > z ? zeN : zeP )/z ) ) ;
   const unsigned int ix  ( xindex( x*fac, z ) ) ; 
   const unsigned int iy  ( yindex( y*fac, z ) ) ; 
   const unsigned int iz  ( z>0 ? 1 : -1 ) ;

   if( EEDetId::validDetId( ix, iy, iz ) ) 
   {
      return EEDetId( ix, iy, iz ) ; // first try is on target
   }
   else // try nearby coordinates, spiraling out from center
   {
      for( unsigned int i ( 1 ) ; i != 6 ; ++i )
      {
         for( unsigned int k ( 0 ) ; k != 8 ; ++k )
         {
            const int jx ( 0 == k || 4 == k || 5 == k ? +i :
                           ( 1 == k || 5 < k ? -i : 0 ) ) ;
            const int jy ( 2 == k || 4 == k || 6 == k ? +i :
                           ( 3 == k || 5 == k || 7 == k ? -i : 0 ) ) ;
            if( EEDetId::validDetId( ix + jx, iy + jy, iz ) ) 
            {
               return EEDetId( ix + jx, iy + jy, iz ) ;
            }
         }
      }
   }
   return EEDetId() ; // nowhere near any crystal
}


// Get closest cell, etc...
DetId 
EcalEndcapGeometry::getClosestCell( const GlobalPoint& r ) const 
{
   try
   {
      EEDetId mycellID ( gId( r.x(), r.y(), r.z() ) ) ; // educated guess

      if( EEDetId::validDetId( mycellID.ix(), 
                               mycellID.iy(),
                               mycellID.zside() ) )
      {
         // now get points in convenient ordering

         HepGeom::Point3D<double>   A;
         HepGeom::Point3D<double>   B;
         HepGeom::Point3D<double>   C;
         HepGeom::Point3D<double>   point(r.x(),r.y(),r.z());
         // D.K. : equation of plane : AA*x+BB*y+CC*z+DD=0;
         // finding equation for each edge
         
         // ================================================================
         double x,y,z;
         unsigned offset=0;
         int zsign=1;
         //================================================================
         std::vector<double> SS;
      
         // compute the distance of the point with respect of the 4 crystal lateral planes

         
         if( 0 != getGeometry(mycellID) )
         {
            const GlobalPoint& myPosition=getGeometry(mycellID)->getPosition();
         
            x=myPosition.x();
            y=myPosition.y();
            z=myPosition.z();
         
            offset=0;
            // This will disappear when Andre has applied his fix
            zsign=1;
         
            if(z>0)
            {
               if(x>0&&y>0)
                  offset=1;
               else  if(x<0&&y>0)
                  offset=2;
               else if(x>0&&y<0)
                  offset=0;
               else if (x<0&&y<0)
                  offset=3;
               zsign=1;
            }
            else
            {
               if(x>0&&y>0)
                  offset=3;
               else if(x<0&&y>0)
                  offset=2;
               else if(x>0&&y<0)
                  offset=0;
               else if(x<0&&y<0)
                  offset=1;
               zsign=-1;
            }
            std::vector<GlobalPoint> corners;
            corners.clear();
            corners.resize(8);
            for(unsigned ic=0;ic<4;++ic)
            {
               corners[ic]=getGeometry(mycellID)->getCorners()[(unsigned)((zsign*ic+offset)%4)];
               corners[4+ic]=getGeometry(mycellID)->getCorners()[(unsigned)(4+(zsign*ic+offset)%4)];
            }
            
            for (short i=0; i < 4 ; ++i)
            {
               A = HepGeom::Point3D<double> (corners[i%4].x(),corners[i%4].y(),corners[i%4].z());
               B = HepGeom::Point3D<double> (corners[(i+1)%4].x(),corners[(i+1)%4].y(),corners[(i+1)%4].z());
               C = HepGeom::Point3D<double> (corners[4+(i+1)%4].x(),corners[4+(i+1)%4].y(),corners[4+(i+1)%4].z());
               HepGeom::Plane3D<double>  plane(A,B,C);
               plane.normalize();
               double distance = plane.distance(point);
               if (corners[0].z()<0.) distance=-distance;
               SS.push_back(distance);
            }
         
            // Only one move in necessary direction
         
            const bool yout ( 0 > SS[0]*SS[2] ) ;
            const bool xout ( 0 > SS[1]*SS[3] ) ;
         
            if( yout || xout )
            {
               const int ydel ( !yout ? 0 :  ( 0 < SS[0] ? -1 : 1 ) ) ;
               const int xdel ( !xout ? 0 :  ( 0 < SS[1] ? -1 : 1 ) ) ;
               const unsigned int ix ( mycellID.ix() + xdel ) ;
               const unsigned int iy ( mycellID.iy() + ydel ) ;
               const unsigned int iz ( mycellID.zside()     ) ;
               if( EEDetId::validDetId( ix, iy, iz ) ) 
                  mycellID = EEDetId( ix, iy, iz ) ;
            }
  
            return mycellID;
         }
         return DetId(0);
      }
   }
   catch ( cms::Exception &e ) 
   { 
      return DetId(0);
   }
   return DetId(0);
}

CaloSubdetectorGeometry::DetIdSet 
EcalEndcapGeometry::getCells( const GlobalPoint& r, 
                              double             dR ) const 
{
   CaloSubdetectorGeometry::DetIdSet dis ; // return object
   if( 0.000001 < dR )
   {
      if( dR > M_PI/2. ) // this code assumes small dR
      {
         dis = CaloSubdetectorGeometry::getCells( r, dR ) ; // base class version
      }
      else
      {
         const double dR2  ( dR*dR ) ;
         const double reta ( r.eta() ) ;
         const double rphi ( r.phi() ) ;
         const double rx   ( r.x() ) ;
         const double ry   ( r.y() ) ;
         const double rz   ( r.z() ) ;
         const double fac  ( fabs( zeP/rz ) ) ;
         const double xx   ( rx*fac ) ; // xyz at endcap z
         const double yy   ( ry*fac ) ; 
         const double zz   ( rz*fac ) ; 

         const double xang  ( atan( xx/zz ) ) ;
         const double lowX  ( zz>0 ? zz*tan( xang - dR ) : zz*tan( xang + dR ) ) ;
         const double highX ( zz>0 ? zz*tan( xang + dR ) : zz*tan( xang - dR ) ) ;
         const double yang  ( atan( yy/zz ) ) ;
         const double lowY  ( zz>0 ? zz*tan( yang - dR ) : zz*tan( yang + dR ) ) ;
         const double highY ( zz>0 ? zz*tan( yang + dR ) : zz*tan( yang - dR ) ) ;

         const double refxlo ( 0 > rz ? m_xlo[0] : m_xlo[1] ) ;
         const double refxhi ( 0 > rz ? m_xhi[0] : m_xhi[1] ) ;
         const double refylo ( 0 > rz ? m_ylo[0] : m_ylo[1] ) ;
         const double refyhi ( 0 > rz ? m_yhi[0] : m_yhi[1] ) ;

         if( lowX  <  refxhi &&   // proceed if any possible overlap with the endcap
             lowY  <  refyhi &&
             highX >  refxlo &&
             highY >  refylo    )
         {
            const int ix_ctr ( xindex( xx, rz ) ) ;
            const int iy_ctr ( yindex( yy, rz ) ) ;
            const int iz     ( rz>0 ? 1 : -1 ) ;
            
            const int ix_hi  ( ix_ctr + int( ( highX - xx )/m_wref ) + 2 ) ;
            const int ix_lo  ( ix_ctr - int( ( xx - lowX  )/m_wref ) - 2 ) ;
            
            const int iy_hi  ( iy_ctr + int( ( highY - yy )/m_wref ) + 2 ) ;
            const int iy_lo  ( iy_ctr - int( ( yy - lowY  )/m_wref ) - 2 ) ;
            
            for( int kx ( ix_lo ) ; kx <= ix_hi ; ++kx ) 
            {
               if( kx >  0      && 
                   kx <= (int) m_nref    )
               {
                  for( int ky ( iy_lo ) ; ky <= iy_hi ; ++ky ) 
                  {
                     if( ky >  0      && 
                         ky <= (int) m_nref    )
                     {
                        if( EEDetId::validDetId( kx, ky, iz ) ) // reject invalid ids
                        {
                           const EEDetId id ( kx, ky, iz ) ;
                           const CaloCellGeometry* cell ( getGeometry( id ) );
                           if( 0 != cell )
                           {
                              const GlobalPoint& p    ( cell->getPosition() ) ;
                              const double       eta  ( p.eta() ) ;
                              const double       phi  ( p.phi() ) ;
                              if( reco::deltaR2( eta, phi, reta, rphi ) < dR2 ) dis.insert( id ) ;
                           }
                        }
                     }
                  }
               }
            }
         }
      }
   }
   return dis;
}

const EcalEndcapGeometry::OrderedListOfEBDetId*
EcalEndcapGeometry::getClosestBarrelCells( EEDetId id ) const
{
   OrderedListOfEBDetId* ptr ( 0 ) ;
   if( 0 != id.rawId() &&
       0 != getGeometry( id ) )
   {
      const float phi ( 370. +
                        getGeometry( id )->getPosition().phi().degrees() );
      const int iPhi ( 1 + int(phi)%360 ) ;
      const int iz ( id.zside() ) ;
      if( 0 == m_borderMgr )
      {
         m_borderMgr = new EZMgrFL<EBDetId>( 720*9, 9 ) ;
      }
      if( 0 == m_borderPtrVec )
      {
         m_borderPtrVec = new VecOrdListEBDetIdPtr() ;
         m_borderPtrVec->reserve( 720 ) ;
         for( unsigned int i ( 0 ) ; i != 720 ; ++i )
         {
            const int kz     ( 360>i ? -1 : 1 ) ;
            const int iEta   ( kz*85 ) ;
            const int iEtam1 ( kz*84 ) ;
            const int iEtam2 ( kz*83 ) ;
            const int jPhi   ( i%360 + 1 ) ;
            OrderedListOfEBDetId& olist ( *new OrderedListOfEBDetId( m_borderMgr ) );
            olist[0]=EBDetId( iEta  ,        jPhi     ) ;
            olist[1]=EBDetId( iEta  , myPhi( jPhi+1 ) ) ;
            olist[2]=EBDetId( iEta  , myPhi( jPhi-1 ) ) ;
            olist[3]=EBDetId( iEtam1,        jPhi     ) ;
            olist[4]=EBDetId( iEtam1, myPhi( jPhi+1 ) ) ;
            olist[5]=EBDetId( iEtam1, myPhi( jPhi-1 ) ) ;
            olist[6]=EBDetId( iEta  , myPhi( jPhi+2 ) ) ;
            olist[7]=EBDetId( iEta  , myPhi( jPhi-2 ) ) ;
            olist[8]=EBDetId( iEtam2,        jPhi     ) ;
            m_borderPtrVec->push_back( &olist ) ;
         }
      }
      ptr = (*m_borderPtrVec)[ ( iPhi - 1 ) + ( 0>iz ? 0 : 360 ) ] ;
   }
   return ptr ;
}

std::vector<HepGeom::Point3D<double> > 
EcalEndcapGeometry::localCorners( const double* pv,
                                  unsigned int  i,
                                  HepGeom::Point3D<double> &   ref )
{
   return ( TruncatedPyramid::localCorners( pv, ref ) ) ;
}

CaloCellGeometry* 
EcalEndcapGeometry::newCell( const GlobalPoint& f1 ,
                             const GlobalPoint& f2 ,
                             const GlobalPoint& f3 ,
                             CaloCellGeometry::CornersMgr* mgr,
                             const double*      parm ,
                             const DetId&       detId   ) 
{
   return ( new TruncatedPyramid( mgr, f1, f2, f3, parm ) ) ;
}


double 
EcalEndcapGeometry::avgAbsZFrontFaceCenter() const
{
   if( 0 > m_avgZ )
   {
      double sum ( 0 ) ;
      const CaloSubdetectorGeometry::CellCont& cells ( cellGeometries() ) ;
      for( unsigned int i ( 0 ) ; i != cells.size() ; ++i )
      {
         sum += fabs( cells[i]->getPosition().z() ) ;
      }
      m_avgZ = sum/cells.size() ;
   }
   return m_avgZ ;
}