00001 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00002 #include "Geometry/CaloGeometry/interface/CaloGenericDetId.h"
00003
00004 #include <Math/Transform3D.h>
00005 #include <Math/EulerAngles.h>
00006
00007 typedef CaloCellGeometry::Pt3D Pt3D ;
00008 typedef CaloCellGeometry::Pt3DVec Pt3DVec ;
00009 typedef CaloCellGeometry::Tr3D Tr3D ;
00010
00011 typedef CaloSubdetectorGeometry::CCGFloat CCGFloat ;
00012
00013 CaloSubdetectorGeometry::CaloSubdetectorGeometry() :
00014 m_parMgr ( 0 ) ,
00015 m_cmgr ( 0 ) ,
00016 m_sortedIds (false) ,
00017 m_deltaPhi ( 0 ) ,
00018 m_deltaEta ( 0 )
00019 {}
00020
00021
00022 CaloSubdetectorGeometry::~CaloSubdetectorGeometry()
00023 {
00024 delete m_cmgr ;
00025 delete m_parMgr ;
00026 delete m_deltaPhi ;
00027 delete m_deltaEta ;
00028 }
00029
00030 const std::vector<DetId>&
00031 CaloSubdetectorGeometry::getValidDetIds( DetId::Detector ,
00032 int ) const
00033 {
00034 if( !m_sortedIds )
00035 {
00036 m_sortedIds = true ;
00037 std::sort( m_validIds.begin(), m_validIds.end() ) ;
00038 }
00039 return m_validIds ;
00040 }
00041
00042 const CaloCellGeometry*
00043 CaloSubdetectorGeometry::getGeometry( const DetId& id ) const
00044 {
00045 return cellGeomPtr( CaloGenericDetId( id ).denseIndex() ) ;
00046 }
00047
00048 bool
00049 CaloSubdetectorGeometry::present( const DetId& id ) const
00050 {
00051 return ( 0 != getGeometry( id ) ) ;
00052 }
00053
00054 DetId
00055 CaloSubdetectorGeometry::getClosestCell( const GlobalPoint& r ) const
00056 {
00057 const CCGFloat eta ( r.eta() ) ;
00058 const CCGFloat phi ( r.phi() ) ;
00059 uint32_t index ( ~0 ) ;
00060 CCGFloat closest ( 1e9 ) ;
00061
00062 for( uint32_t i ( 0 ); i != m_validIds.size() ; ++i )
00063 {
00064 const CaloCellGeometry* cell ( getGeometry( m_validIds[ i ] ) ) ;
00065 if( 0 != cell )
00066 {
00067 const GlobalPoint& p ( cell->getPosition() ) ;
00068 const CCGFloat eta0 ( p.eta() ) ;
00069 const CCGFloat phi0 ( p.phi() ) ;
00070 const CCGFloat dR2 ( reco::deltaR2( eta0, phi0, eta, phi ) ) ;
00071 if( dR2 < closest )
00072 {
00073 closest = dR2 ;
00074 index = i ;
00075 }
00076 }
00077 }
00078 return ( closest > 0.9e9 ||
00079 (uint32_t)(~0) == index ? DetId(0) :
00080 m_validIds[index] ) ;
00081 }
00082
00083 CaloSubdetectorGeometry::DetIdSet
00084 CaloSubdetectorGeometry::getCells( const GlobalPoint& r,
00085 double dR ) const
00086 {
00087 const double dR2 ( dR*dR ) ;
00088 const double eta ( r.eta() ) ;
00089 const double phi ( r.phi() ) ;
00090
00091 DetIdSet dss;
00092
00093 if( 0.000001 < dR )
00094 {
00095 for( uint32_t i ( 0 ); i != m_validIds.size() ; ++i )
00096 {
00097 const CaloCellGeometry* cell ( getGeometry( m_validIds[i] ) ) ;
00098 if( 0 != cell )
00099 {
00100 const GlobalPoint& p ( cell->getPosition() ) ;
00101 const CCGFloat eta0 ( p.eta() ) ;
00102 if( fabs( eta - eta0 ) < dR )
00103 {
00104 const CCGFloat phi0 ( p.phi() ) ;
00105 CCGFloat delp ( fabs( phi - phi0 ) ) ;
00106 if( delp > M_PI ) delp = 2*M_PI - delp ;
00107 if( delp < dR )
00108 {
00109 const CCGFloat dist2 ( reco::deltaR2( eta0, phi0, eta, phi ) ) ;
00110 if( dist2 < dR2 ) dss.insert( m_validIds[i] ) ;
00111 }
00112 }
00113 }
00114 }
00115 }
00116 return dss;
00117 }
00118
00119 CaloSubdetectorGeometry::CellSet
00120 CaloSubdetectorGeometry::getCellSet( const GlobalPoint& r, double dR ) const {
00121
00122 DetIdSet ids = getCells(r, dR);
00123 CellSet cells; cells.reserve(ids.size());
00124 for ( auto id : ids) cells.push_back(getGeometry(id));
00125 return cells;
00126 }
00127
00128 void
00129 CaloSubdetectorGeometry::allocateCorners( CaloCellGeometry::CornersVec::size_type n )
00130 {
00131 assert( 0 == m_cmgr ) ;
00132 m_cmgr = new CaloCellGeometry::CornersMgr( n*( CaloCellGeometry::k_cornerSize ),
00133 CaloCellGeometry::k_cornerSize ) ;
00134
00135 m_validIds.reserve( n ) ;
00136 }
00137
00138 void
00139 CaloSubdetectorGeometry::allocatePar( ParVec::size_type n,
00140 unsigned int m )
00141 {
00142 assert( 0 == m_parMgr ) ;
00143 m_parMgr = new ParMgr( n*m, m ) ;
00144 }
00145
00146 void
00147 CaloSubdetectorGeometry::getSummary( CaloSubdetectorGeometry::TrVec& tVec ,
00148 CaloSubdetectorGeometry::IVec& iVec ,
00149 CaloSubdetectorGeometry::DimVec& dVec ,
00150 CaloSubdetectorGeometry::IVec& ) const
00151 {
00152 tVec.reserve( m_validIds.size()*numberOfTransformParms() ) ;
00153 iVec.reserve( numberOfShapes()==1 ? 1 : m_validIds.size() ) ;
00154 dVec.reserve( numberOfShapes()*numberOfParametersPerShape() ) ;
00155
00156 for( ParVecVec::const_iterator ivv ( parVecVec().begin() ) ; ivv != parVecVec().end() ; ++ivv )
00157 {
00158 const ParVec& pv ( *ivv ) ;
00159 for( ParVec::const_iterator iv ( pv.begin() ) ; iv != pv.end() ; ++iv )
00160 {
00161 dVec.push_back( *iv ) ;
00162 }
00163 }
00164
00165 for( uint32_t i ( 0 ) ; i != m_validIds.size() ; ++i )
00166 {
00167 Tr3D tr ;
00168 const CaloCellGeometry* ptr ( cellGeomPtr( i ) ) ;
00169 assert( 0 != ptr ) ;
00170 ptr->getTransform( tr, ( Pt3DVec* ) 0 ) ;
00171
00172 if( Tr3D() == tr )
00173 {
00174 const GlobalPoint& gp ( ptr->getPosition() ) ;
00175 tr = HepGeom::Translate3D( gp.x(), gp.y(), gp.z() ) ;
00176 }
00177
00178 const CLHEP::Hep3Vector tt ( tr.getTranslation() ) ;
00179 tVec.push_back( tt.x() ) ;
00180 tVec.push_back( tt.y() ) ;
00181 tVec.push_back( tt.z() ) ;
00182 if( 6 == numberOfTransformParms() )
00183 {
00184 const CLHEP::HepRotation rr ( tr.getRotation() ) ;
00185 const ROOT::Math::Transform3D rtr ( rr.xx(), rr.xy(), rr.xz(), tt.x(),
00186 rr.yx(), rr.yy(), rr.yz(), tt.y(),
00187 rr.zx(), rr.zy(), rr.zz(), tt.z() ) ;
00188 ROOT::Math::EulerAngles ea ;
00189 rtr.GetRotation( ea ) ;
00190 tVec.push_back( ea.Phi() ) ;
00191 tVec.push_back( ea.Theta() ) ;
00192 tVec.push_back( ea.Psi() ) ;
00193 }
00194
00195 const CCGFloat* par ( ptr->param() ) ;
00196
00197 unsigned int ishape ( 9999 ) ;
00198 for( unsigned int ivv ( 0 ) ; ivv != parVecVec().size() ; ++ivv )
00199 {
00200 bool ok ( true ) ;
00201 const CCGFloat* pv ( &(*parVecVec()[ivv].begin() ) ) ;
00202 for( unsigned int k ( 0 ) ; k != numberOfParametersPerShape() ; ++k )
00203 {
00204 ok = ok && ( fabs( par[k] - pv[k] ) < 1.e-6 ) ;
00205 }
00206 if( ok )
00207 {
00208 ishape = ivv ;
00209 break ;
00210 }
00211 }
00212 assert( 9999 != ishape ) ;
00213
00214 const unsigned int nn (( numberOfShapes()==1) ? (unsigned int)1 : m_validIds.size() ) ;
00215 if( iVec.size() < nn ) iVec.push_back( ishape ) ;
00216 }
00217 }
00218
00219 CCGFloat
00220 CaloSubdetectorGeometry::deltaPhi( const DetId& detId ) const
00221 {
00222 const CaloGenericDetId cgId ( detId ) ;
00223
00224 if( 0 == m_deltaPhi )
00225 {
00226 const uint32_t kSize ( sizeForDenseIndex(detId));
00227 m_deltaPhi = new std::vector<CCGFloat> ( kSize ) ;
00228 for( uint32_t i ( 0 ) ; i != kSize ; ++i )
00229 {
00230 const CaloCellGeometry* cellPtr ( cellGeomPtr( i ) ) ;
00231 if( 0 != cellPtr )
00232 {
00233 const CaloCellGeometry& cell ( *cellPtr ) ;
00234 CCGFloat dPhi1 ( fabs(
00235 GlobalPoint( ( cell.getCorners()[0].x() +
00236 cell.getCorners()[1].x() )/2. ,
00237 ( cell.getCorners()[0].y() +
00238 cell.getCorners()[1].y() )/2. ,
00239 ( cell.getCorners()[0].z() +
00240 cell.getCorners()[1].z() )/2. ).phi() -
00241 GlobalPoint( ( cell.getCorners()[2].x() +
00242 cell.getCorners()[3].x() )/2. ,
00243 ( cell.getCorners()[2].y() +
00244 cell.getCorners()[3].y() )/2. ,
00245 ( cell.getCorners()[2].z() +
00246 cell.getCorners()[3].z() )/2. ).phi() ) ) ;
00247 CCGFloat dPhi2 ( fabs(
00248 GlobalPoint( ( cell.getCorners()[0].x() +
00249 cell.getCorners()[3].x() )/2. ,
00250 ( cell.getCorners()[0].y() +
00251 cell.getCorners()[3].y() )/2. ,
00252 ( cell.getCorners()[0].z() +
00253 cell.getCorners()[3].z() )/2. ).phi() -
00254 GlobalPoint( ( cell.getCorners()[2].x() +
00255 cell.getCorners()[1].x() )/2. ,
00256 ( cell.getCorners()[2].y() +
00257 cell.getCorners()[1].y() )/2. ,
00258 ( cell.getCorners()[2].z() +
00259 cell.getCorners()[1].z() )/2. ).phi() ) ) ;
00260 if( M_PI < dPhi1 ) dPhi1 = fabs( dPhi1 - 2.*M_PI ) ;
00261 if( M_PI < dPhi2 ) dPhi2 = fabs( dPhi2 - 2.*M_PI ) ;
00262 (*m_deltaPhi)[i] = dPhi1>dPhi2 ? dPhi1 : dPhi2 ;
00263 }
00264 }
00265 }
00266 return (*m_deltaPhi)[ indexFor(detId) ] ;
00267 }
00268
00269 CCGFloat
00270 CaloSubdetectorGeometry::deltaEta( const DetId& detId ) const
00271 {
00272
00273 if( 0 == m_deltaEta )
00274 {
00275 const uint32_t kSize ( sizeForDenseIndex(detId));
00276 m_deltaEta = new std::vector<CCGFloat> ( kSize ) ;
00277 for( uint32_t i ( 0 ) ; i != kSize ; ++i )
00278 {
00279 const CaloCellGeometry* cellPtr ( cellGeomPtr( i ) ) ;
00280 if( 0 != cellPtr )
00281 {
00282 const CaloCellGeometry& cell ( *cellPtr ) ;
00283 const CCGFloat dEta1 ( fabs(
00284 GlobalPoint( ( cell.getCorners()[0].x() +
00285 cell.getCorners()[1].x() )/2. ,
00286 ( cell.getCorners()[0].y() +
00287 cell.getCorners()[1].y() )/2. ,
00288 ( cell.getCorners()[0].z() +
00289 cell.getCorners()[1].z() )/2. ).eta() -
00290 GlobalPoint( ( cell.getCorners()[2].x() +
00291 cell.getCorners()[3].x() )/2. ,
00292 ( cell.getCorners()[2].y() +
00293 cell.getCorners()[3].y() )/2. ,
00294 ( cell.getCorners()[2].z() +
00295 cell.getCorners()[3].z() )/2. ).eta() ) ) ;
00296 const CCGFloat dEta2 ( fabs(
00297 GlobalPoint( ( cell.getCorners()[0].x() +
00298 cell.getCorners()[3].x() )/2. ,
00299 ( cell.getCorners()[0].y() +
00300 cell.getCorners()[3].y() )/2. ,
00301 ( cell.getCorners()[0].z() +
00302 cell.getCorners()[3].z() )/2. ).eta() -
00303 GlobalPoint( ( cell.getCorners()[2].x() +
00304 cell.getCorners()[1].x() )/2. ,
00305 ( cell.getCorners()[2].y() +
00306 cell.getCorners()[1].y() )/2. ,
00307 ( cell.getCorners()[2].z() +
00308 cell.getCorners()[1].z() )/2. ).eta() ) ) ;
00309 (*m_deltaEta)[i] = dEta1>dEta2 ? dEta1 : dEta2 ;
00310 }
00311 }
00312 }
00313 return (*m_deltaEta)[ indexFor(detId)];
00314 }
00315
00316
00317 unsigned int CaloSubdetectorGeometry::indexFor(const DetId& id) const { return CaloGenericDetId(id).denseIndex(); }
00318
00319 unsigned int CaloSubdetectorGeometry::sizeForDenseIndex(const DetId& id) const { return CaloGenericDetId(id).sizeForDenseIndexing(); }
1.7.3