CaloSubdetectorGeometry.cc

Go to the documentation of this file.

#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloGenericDetId.h"

#include <Math/Transform3D.h>
#include <Math/EulerAngles.h>

#include <algorithm> 

typedef CaloCellGeometry::Pt3D     Pt3D     ;
typedef CaloCellGeometry::Pt3DVec  Pt3DVec  ;
typedef CaloCellGeometry::Tr3D     Tr3D     ;

typedef CaloSubdetectorGeometry::CCGFloat CCGFloat ;

CaloSubdetectorGeometry::CaloSubdetectorGeometry() : 
   m_parMgr ( nullptr ) ,
   m_cmgr   ( nullptr ) ,
   m_deltaPhi  (nullptr) ,
   m_deltaEta  (nullptr)
{}


CaloSubdetectorGeometry::~CaloSubdetectorGeometry() 
{ 
   delete m_cmgr ;
   delete m_parMgr ; 
   if (m_deltaPhi) delete m_deltaPhi.load() ;
   if (m_deltaEta) delete m_deltaEta.load() ;
}

void
CaloSubdetectorGeometry::addValidID(const DetId& id)
{
    auto pos = std::lower_bound(m_validIds.begin(), m_validIds.end(), id);
    m_validIds.insert(pos, id);
}

const std::vector<DetId>& 
CaloSubdetectorGeometry::getValidDetIds( DetId::Detector /*det*/    , 
                     int             /*subdet*/   ) const 
{
   return m_validIds ;
}

std::shared_ptr<const CaloCellGeometry> 
CaloSubdetectorGeometry::getGeometry( const DetId& id ) const
{
  return cellGeomPtr(CaloGenericDetId(id).denseIndex());
}

bool 
CaloSubdetectorGeometry::present( const DetId& id ) const {
  return std::find(m_validIds.begin(),m_validIds.end(),id)!=m_validIds.end();
}

DetId 
CaloSubdetectorGeometry::getClosestCell( const GlobalPoint& r ) const {
  const CCGFloat eta ( r.eta() ) ;
  const CCGFloat phi ( r.phi() ) ;
  uint32_t index ( ~0 ) ;
  CCGFloat closest ( 1e9 ) ;

  for( uint32_t i ( 0 ); i != m_validIds.size() ; ++i ) {
    std::shared_ptr<const CaloCellGeometry> cell ( getGeometry( m_validIds[ i ] ) ) ;
    if( nullptr != cell ) {
      const GlobalPoint& p ( cell->getPosition() ) ;
      const CCGFloat eta0 ( p.eta() ) ;
      const CCGFloat phi0 ( p.phi() ) ;
      const CCGFloat dR2 ( reco::deltaR2( eta0, phi0, eta, phi ) ) ;
      if( dR2 < closest ) {
    closest = dR2 ;
    index   = i   ;
      }
    }
  }
  return ( closest > 0.9e9 ||
       (uint32_t)(~0) == index       ? DetId(0) :
       m_validIds[index] ) ;
}

CaloSubdetectorGeometry::DetIdSet 
CaloSubdetectorGeometry::getCells(const GlobalPoint& r, double dR) const {
   const double dR2 ( dR*dR ) ;
   const double eta ( r.eta() ) ;
   const double phi ( r.phi() ) ;

   DetIdSet dss;
   
   if( 0.000001 < dR )
   {
      for( uint32_t i ( 0 ); i != m_validIds.size() ; ++i ) 
      {
     std::shared_ptr<const CaloCellGeometry> cell ( getGeometry( m_validIds[i] ) ) ;
     if( nullptr != cell )
     {
        const GlobalPoint& p ( cell->getPosition() ) ;
        const CCGFloat eta0 ( p.eta() ) ;
        if( fabs( eta - eta0 ) < dR )
        {
           const CCGFloat phi0 ( p.phi() ) ;
           CCGFloat delp ( fabs( phi - phi0 ) ) ;
           if( delp > M_PI ) delp = 2*M_PI - delp ;
           if( delp < dR )
           {
          const CCGFloat dist2 ( reco::deltaR2( eta0, phi0, eta, phi ) ) ;
          if( dist2 < dR2 ) dss.insert( m_validIds[i] ) ;
           }
        }
     }
      }   
   }
   return dss;
}

CaloSubdetectorGeometry::CellSet 
CaloSubdetectorGeometry::getCellSet( const GlobalPoint& r, double dR ) const {
  // stupid implementation not to be really used...
  DetIdSet ids = getCells(r, dR);
  CellSet cells; cells.reserve(ids.size());
  for ( auto id : ids) cells.emplace_back(getGeometry(id));
  return cells;
}

void 
CaloSubdetectorGeometry::allocateCorners( CaloCellGeometry::CornersVec::size_type n )
{
   assert( nullptr == m_cmgr ) ;
   m_cmgr = new CaloCellGeometry::CornersMgr( n*( CaloCellGeometry::k_cornerSize ),
                          CaloCellGeometry::k_cornerSize        ) ; 

   m_validIds.reserve( n ) ;
}

void 
CaloSubdetectorGeometry::allocatePar( ParVec::size_type n,
                      unsigned int      m     )
{
   assert( nullptr == m_parMgr ) ;
   m_parMgr = new ParMgr( n*m, m ) ;
}

void
CaloSubdetectorGeometry::getSummary( CaloSubdetectorGeometry::TrVec&  tVec ,
                     CaloSubdetectorGeometry::IVec&   iVec ,   
                     CaloSubdetectorGeometry::DimVec& dVec ,
                     CaloSubdetectorGeometry::IVec& /*dins*/) const {
  tVec.reserve( m_validIds.size()*numberOfTransformParms() ) ;
  iVec.reserve( numberOfShapes()==1 ? 1 : m_validIds.size() ) ;
  dVec.reserve( numberOfShapes()*numberOfParametersPerShape() ) ;

  for(const auto & pv : parVecVec())  {
    for(float iv : pv) {
      dVec.emplace_back( iv ) ;
    }
  }

  for( uint32_t i ( 0 ) ; i != m_validIds.size() ; ++i ) {
    Tr3D tr ;
    std::shared_ptr<const CaloCellGeometry> ptr ( cellGeomPtr( i ) ) ;
    assert( nullptr != ptr ) ;
    ptr->getTransform( tr, ( Pt3DVec* ) nullptr ) ;
    
    if( Tr3D() == tr ) {// for preshower there is no rotation
      const GlobalPoint& gp ( ptr->getPosition() ) ; 
      tr = HepGeom::Translate3D( gp.x(), gp.y(), gp.z() ) ;
    }

    const CLHEP::Hep3Vector  tt ( tr.getTranslation() ) ;
    tVec.emplace_back( tt.x() ) ;
    tVec.emplace_back( tt.y() ) ;
    tVec.emplace_back( tt.z() ) ;
    if( 6 == numberOfTransformParms() ) {
      const CLHEP::HepRotation rr ( tr.getRotation() ) ;
      const ROOT::Math::Transform3D rtr ( rr.xx(), rr.xy(), rr.xz(), tt.x(),
                      rr.yx(), rr.yy(), rr.yz(), tt.y(),
                      rr.zx(), rr.zy(), rr.zz(), tt.z()  ) ;
      ROOT::Math::EulerAngles ea ;
      rtr.GetRotation( ea ) ;
      tVec.emplace_back( ea.Phi() ) ;
      tVec.emplace_back( ea.Theta() ) ;
      tVec.emplace_back( ea.Psi() ) ;
    }

    const CCGFloat* par ( ptr->param() ) ;
    
    unsigned int ishape ( 9999 ) ;
    for( unsigned int ivv ( 0 ) ; ivv != parVecVec().size() ; ++ivv ) {
      bool ok ( true ) ;
      const CCGFloat* pv ( &(*parVecVec()[ivv].begin() ) ) ;
      for( unsigned int k ( 0 ) ; k != numberOfParametersPerShape() ; ++k ) {
    ok = ok && ( fabs( par[k] - pv[k] ) < 1.e-6 ) ;
      }
      if( ok ) {
    ishape = ivv ;
    break ;
      }
    }
    assert( 9999 != ishape ) ;

    const unsigned int nn (( numberOfShapes()==1) ? (unsigned int)1 : m_validIds.size() ) ; 
    if( iVec.size() < nn ) iVec.emplace_back( ishape ) ;
  }
}

CCGFloat CaloSubdetectorGeometry::deltaPhi( const DetId& detId ) const {
  const CaloGenericDetId cgId ( detId ) ;
  
  if(!m_deltaPhi.load(std::memory_order_acquire))  {
    const uint32_t kSize ( sizeForDenseIndex(detId));
    auto ptr = new std::vector<CCGFloat>(kSize);
    for( uint32_t i ( 0 ) ; i != kSize ; ++i ) {
      std::shared_ptr<const CaloCellGeometry> cellPtr ( cellGeomPtr( i ) ) ;
      if( nullptr != cellPtr ) {
    CCGFloat dPhi1 (fabs(
                 GlobalPoint((cellPtr->getCorners()[0].x() + 
                      cellPtr->getCorners()[1].x() )/2.,
                     (cellPtr->getCorners()[0].y() + 
                      cellPtr->getCorners()[1].y() )/2.,
                     (cellPtr->getCorners()[0].z() + 
                      cellPtr->getCorners()[1].z() )/2.).phi() -
                 GlobalPoint((cellPtr->getCorners()[2].x() + 
                      cellPtr->getCorners()[3].x() )/2.,
                     (cellPtr->getCorners()[2].y() + 
                      cellPtr->getCorners()[3].y() )/2.,
                     (cellPtr->getCorners()[2].z() + 
                      cellPtr->getCorners()[3].z() )/2.).phi() ) ) ;
    CCGFloat dPhi2 (fabs(
                 GlobalPoint((cellPtr->getCorners()[0].x() + 
                      cellPtr->getCorners()[3].x() )/2.,
                     (cellPtr->getCorners()[0].y() + 
                      cellPtr->getCorners()[3].y() )/2.,
                     (cellPtr->getCorners()[0].z() + 
                      cellPtr->getCorners()[3].z() )/2. ).phi() -
                 GlobalPoint((cellPtr->getCorners()[2].x() + 
                      cellPtr->getCorners()[1].x() )/2.,
                     (cellPtr->getCorners()[2].y() + 
                      cellPtr->getCorners()[1].y() )/2.,
                     (cellPtr->getCorners()[2].z() + 
                      cellPtr->getCorners()[1].z() )/2.).phi() ) ) ;
        if( M_PI < dPhi1 ) dPhi1 = fabs( dPhi1 - 2.*M_PI ) ;
        if( M_PI < dPhi2 ) dPhi2 = fabs( dPhi2 - 2.*M_PI ) ;
        (*ptr)[i] = dPhi1>dPhi2 ? dPhi1 : dPhi2 ;
     }
      }
    std::vector<CCGFloat>* expect = nullptr;
    bool exchanged = m_deltaPhi.compare_exchange_strong(expect, ptr, std::memory_order_acq_rel);
    if (!exchanged) delete ptr;
   }
   return (*m_deltaPhi.load(std::memory_order_acquire))[ indexFor(detId) ] ;
}

CCGFloat 
CaloSubdetectorGeometry::deltaEta( const DetId& detId ) const {

  if(!m_deltaEta.load(std::memory_order_acquire)) {
    const uint32_t kSize ( sizeForDenseIndex(detId));
    auto ptr = new std::vector<CCGFloat> ( kSize ) ;
    for( uint32_t i ( 0 ) ; i != kSize ; ++i ) {
      std::shared_ptr<const CaloCellGeometry> cellPtr ( cellGeomPtr( i ) ) ;
      if( nullptr != cellPtr ) {
    const CCGFloat dEta1(fabs(
                  GlobalPoint((cellPtr->getCorners()[0].x() + 
                           cellPtr->getCorners()[1].x())/2.,
                          (cellPtr->getCorners()[0].y() + 
                           cellPtr->getCorners()[1].y())/2.,
                          (cellPtr->getCorners()[0].z() + 
                           cellPtr->getCorners()[1].z())/2.).eta() -
                  GlobalPoint((cellPtr->getCorners()[2].x() + 
                           cellPtr->getCorners()[3].x())/2.,
                          (cellPtr->getCorners()[2].y() + 
                           cellPtr->getCorners()[3].y())/2. ,
                          (cellPtr->getCorners()[2].z() + 
                           cellPtr->getCorners()[3].z())/2.).eta() ) ) ;
    const CCGFloat dEta2(fabs(
                  GlobalPoint((cellPtr->getCorners()[0].x() + 
                           cellPtr->getCorners()[3].x())/2.,
                          (cellPtr->getCorners()[0].y() + 
                           cellPtr->getCorners()[3].y())/2.,
                          (cellPtr->getCorners()[0].z() + 
                           cellPtr->getCorners()[3].z())/2.).eta() -
                  GlobalPoint((cellPtr->getCorners()[2].x() + 
                           cellPtr->getCorners()[1].x())/2.,
                          (cellPtr->getCorners()[2].y() + 
                           cellPtr->getCorners()[1].y())/2.,
                          (cellPtr->getCorners()[2].z() + 
                           cellPtr->getCorners()[1].z())/2.).eta() ) ) ;
    (*ptr)[i] = dEta1>dEta2 ? dEta1 : dEta2 ;
      }
    }
    std::vector<CCGFloat>* expect = nullptr;
    bool exchanged = m_deltaEta.compare_exchange_strong(expect, ptr, std::memory_order_acq_rel);
    if (!exchanged) delete ptr;
   }
   return (*m_deltaEta.load(std::memory_order_acquire))[ indexFor(detId)];
}


unsigned int CaloSubdetectorGeometry::indexFor(const DetId& id) const { return CaloGenericDetId(id).denseIndex(); }

unsigned int CaloSubdetectorGeometry::sizeForDenseIndex(const DetId& id) const { return CaloGenericDetId(id).sizeForDenseIndexing(); }

std::shared_ptr<const CaloCellGeometry> 
CaloSubdetectorGeometry::cellGeomPtr(uint32_t index) const {
  // Default version
  auto ptr = getGeometryRawPtr(index);
  static const auto do_not_delete = [](const void*){};
  return ptr == nullptr ? nullptr : std::shared_ptr<const CaloCellGeometry>(ptr, do_not_delete);
}