/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/Geometry/HcalTowerAlgo/src/HcalGeometry.cc

Go to the documentation of this file.

#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
#include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
#include "Geometry/HcalTowerAlgo//src/HcalHardcodeGeometryData.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include <algorithm>

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

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

HcalGeometry::HcalGeometry( const HcalTopology& topology ) :
  theTopology( topology ) {
  init();
}
  

HcalGeometry::~HcalGeometry() {}


void HcalGeometry::init() {
  edm::LogInfo("HcalGeometry") << "HcalGeometry::init() "
                               << " HBSize " << theTopology.getHBSize() 
                               << " HESize " << theTopology.getHESize() 
                               << " HOSize " << theTopology.getHOSize() 
                               << " HFSize " << theTopology.getHFSize();
    
  m_hbCellVec = HBCellVec( theTopology.getHBSize() ) ;
  m_heCellVec = HECellVec( theTopology.getHESize() ) ;
  m_hoCellVec = HOCellVec( theTopology.getHOSize() ) ;
  m_hfCellVec = HFCellVec( theTopology.getHFSize() ) ;
}

void
HcalGeometry::fillDetIds() const
{
   const std::vector<DetId>& baseIds ( CaloSubdetectorGeometry::getValidDetIds() ) ;
   for( unsigned int i ( 0 ) ; i != baseIds.size() ; ++i ) 
   {
      const DetId id ( baseIds[i] );
      if( id.subdetId() == HcalBarrel )
      { 
         m_hbIds.push_back( id ) ;
      }
      else
      {
         if( id.subdetId() == HcalEndcap )
         { 
            m_heIds.push_back( id ) ;
         }
         else
         {
            if( id.subdetId() == HcalOuter )
            { 
               m_hoIds.push_back( id ) ;
            }
            else
            {
               if( id.subdetId() == HcalForward )
               { 
                  m_hfIds.push_back( id ) ;
               }
            }
         }
      }
   }
   std::sort( m_hbIds.begin(), m_hbIds.end() ) ;
   std::sort( m_heIds.begin(), m_heIds.end() ) ;
   std::sort( m_hoIds.begin(), m_hoIds.end() ) ;
   std::sort( m_hfIds.begin(), m_hfIds.end() ) ;
       
   m_emptyIds.resize( 0 ) ;
}

const std::vector<DetId>& 
HcalGeometry::getValidDetIds( DetId::Detector det,
                              int             subdet ) const 
{
   if( 0 != subdet &&
       0 == m_hbIds.size() ) fillDetIds() ;
   return ( 0 == subdet ? CaloSubdetectorGeometry::getValidDetIds() :
            ( HcalBarrel == subdet ? m_hbIds :
              ( HcalEndcap == subdet ? m_heIds :
                ( HcalOuter == subdet ? m_hoIds :
                  ( HcalForward == subdet ? m_hfIds : m_emptyIds ) ) ) ) ) ;
}

DetId HcalGeometry::getClosestCell(const GlobalPoint& r) const {

  // Now find the closest eta_bin, eta value of a bin i is average
  // of eta[i] and eta[i-1]
  double abseta = fabs(r.eta());
  
  // figure out subdetector, giving preference to HE in HE/HF overlap region
  HcalSubdetector bc= HcalEmpty;
  if (abseta <= theHBHEEtaBounds[theTopology.lastHBRing()] ) {
    bc = HcalBarrel;
  } else if (abseta <= theHBHEEtaBounds[theTopology.lastHERing()] ) {
    bc = HcalEndcap;
  } else {
    bc = HcalForward;
  }

  if (bc == HcalForward) {
    static const double z_short=1137.0;
    int etaring = etaRing(bc, abseta);  // This is safer
    /*
      static const double z_long=1115.0;
      // determine front-face eta
      double radius=sqrt(pow(r.x(),2)+pow(r.y(),2));
      double trueAeta=asinh(z_long/radius);
      // find eta bin
      int etaring = etaRing(bc, trueAeta);
    */
    if (etaring>theTopology.lastHFRing()) etaring=theTopology.lastHFRing(); 
  
    int phibin = phiBin(r.phi(), etaring);

    // add a sign to the etaring
    int etabin = (r.z() > 0) ? etaring : -etaring;
    // Next line is premature depth 1 and 2 can coexist for large z-extent
    //    HcalDetId bestId(bc,etabin,phibin,((fabs(r.z())>=z_short)?(2):(1)));
    // above line is no good with finite precision
    HcalDetId bestId(bc,etabin,phibin,((fabs(r.z()) - z_short >-0.1)?(2):(1)));
    return bestId;
  } else {

    // find eta bin
    int etaring = etaRing(bc, abseta);
    
    int phibin = phiBin(r.phi(), etaring);
    
    // add a sign to the etaring
    int etabin = (r.z() > 0) ? etaring : -etaring;
    
    //Now do depth if required
    int dbin = 1;
    double pointrz=0, drz=99999.;
    HcalDetId currentId(bc, etabin, phibin, dbin);
    if (bc == HcalBarrel) pointrz = r.mag();
    else                  pointrz = std::abs(r.z());
    HcalDetId bestId;
//    std::cout << "Current ID " << currentId << std::endl;
    for ( ; currentId != HcalDetId(); theTopology.incrementDepth(currentId)) {
//      std::cout << "Incremented Current ID " << currentId << std::endl;
      const CaloCellGeometry * cell = getGeometry(currentId);
      if (cell == 0) {
//      std::cout << "Cell 0 for " << currentId << " Best " << bestId << " dummy " << HcalDetId() << std::endl;
        assert (bestId != HcalDetId());
        break;
      } else {
        double rz;
        if (bc == HcalEndcap) rz = std::abs(cell->getPosition().z());
        else                  rz = cell->getPosition().mag();
        if (std::abs(pointrz-rz)<drz) {
          bestId = currentId;
          drz    = std::abs(pointrz-rz);
        }
      }
    }
    
    return bestId;
  }
}


int HcalGeometry::etaRing(HcalSubdetector bc, double abseta) const
{
  int etaring;
  if( bc == HcalForward ) {
    for(etaring = theTopology.firstHFRing();
        etaring <= theTopology.lastHFRing(); ++etaring)
    {
      if(theHFEtaBounds[etaring-theTopology.firstHFRing()+1] > abseta) break;
    }
  }
  else
  {
    for(etaring = 1;
        etaring <= theTopology.lastHERing(); ++etaring)
    {
      if(theHBHEEtaBounds[etaring] >= abseta) break;
    }
  }

  return etaring;
}


int HcalGeometry::phiBin(double phi, int etaring) const
{
   static const double twopi = M_PI+M_PI;
  //put phi in correct range (0->2pi)
  if(phi<0.0) phi += twopi;
  if(phi>twopi) phi -= twopi;
  int nphibins = theTopology.nPhiBins(etaring);
  int phibin= static_cast<int>(phi/twopi*nphibins)+1;
  int iphi;

  // rings 40 and 41 are offset wrt the other phi numbering
  //  1        1         1         2
  //  ------------------------------
  //  72       36        36        1
  if(etaring >= theTopology.firstHFQuadPhiRing())
  {
    phi+=(twopi/36); //shift by half tower.    
    phibin=static_cast<int>(phi/twopi*nphibins);
    if (phibin==0) phibin=18;
    iphi=phibin*4-1; // 71,3,5,
  } else {
    // convert to the convention of numbering 1,3,5, in 36 phi bins
    iphi=(phibin-1)*(72/nphibins) + 1;
  }

  return iphi;
}

CaloSubdetectorGeometry::DetIdSet 
HcalGeometry::getCells( const GlobalPoint& r, 
                        double             dR ) const 
{
   CaloSubdetectorGeometry::DetIdSet dis;  // this is the return object

   if( 0.000001 < dR )
   {
      if( dR > M_PI/2. ) // this version needs "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 lowEta  ( reta - dR ) ;
         const double highEta ( reta + dR ) ;
         const double lowPhi  ( rphi - dR ) ;
         const double highPhi ( rphi + dR ) ;
         
         const double hfEtaHi ( theHFEtaBounds[ theTopology.lastHFRing() -
                                                theTopology.firstHFRing() + 1 ] ) ;
         
         if( highEta > -hfEtaHi &&
             lowEta  <  hfEtaHi    ) // in hcal
         {
            const HcalSubdetector hs[] = { HcalBarrel, HcalOuter, HcalEndcap, HcalForward } ;

            for( unsigned int is ( 0 ) ; is != 4 ; ++is )
            {
               const int sign        (  reta>0 ? 1 : -1 ) ;
               const int ieta_center ( sign*etaRing( hs[is], fabs( reta ) ) ) ;
               const int ieta_lo     ( ( 0 < lowEta*sign ? sign : -sign )*etaRing( hs[is], fabs( lowEta ) ) ) ;
               const int ieta_hi     ( ( 0 < highEta*sign ? sign : -sign )*etaRing( hs[is], fabs( highEta ) ) ) ;
               const int iphi_lo     ( phiBin( lowPhi , ieta_center ) ) ;
               const int iphi_hi     ( phiBin( highPhi, ieta_center ) ) ;
               const int jphi_lo     ( iphi_lo>iphi_hi ? iphi_lo - 72 : iphi_lo ) ;
               const int jphi_hi     ( iphi_hi ) ;

               const int idep_lo     ( 1 == is ? 4 : 1 ) ;
               const int idep_hi     ( 1 == is ? 4 :
                                       ( 2 == is ? 3 : 2 ) ) ;
               for( int ieta ( ieta_lo ) ; ieta <= ieta_hi ; ++ieta ) // over eta limits
               {
                  if( ieta != 0 )
                  {
                     for( int jphi ( jphi_lo ) ; jphi <= jphi_hi ; ++jphi )  // over phi limits
                     {
                        const int iphi ( 1 > jphi ? jphi+72 : jphi ) ;

                        for( int idep ( idep_lo ) ; idep <= idep_hi ; ++idep )
                          {
                            const HcalDetId did ( hs[is], ieta, iphi, idep ) ;
                            if( theTopology.valid(did) ) 

                              {
                              const CaloCellGeometry* cell ( getGeometry( did ) );
                              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( did ) ;
                              }
                           }
                        }
                     }
                  }
               }
            }
         }
      }
   }
   return dis;
}



DetId 
HcalGeometry::detIdFromBarrelAlignmentIndex( unsigned int i )
{
   assert( i < numberOfBarrelAlignments() ) ;
   const int ieta  ( i < numberOfBarrelAlignments()/2 ? -1 : 1 ) ;
   const int iphi ( 1 + (4*i)%72 ) ;
   return HcalDetId( HcalBarrel, ieta, iphi, 1 ) ;
}

DetId 
HcalGeometry::detIdFromEndcapAlignmentIndex( unsigned int i )
{
   assert( i < numberOfEndcapAlignments() ) ;
   const int ieta  ( i < numberOfEndcapAlignments()/2 ? -16 : 16 ) ;
   const int iphi ( 1 + (4*i)%72 ) ;
   return HcalDetId( HcalEndcap, ieta, iphi, 1 ) ;
}

DetId 
HcalGeometry::detIdFromForwardAlignmentIndex(    unsigned int i )
{
   assert( i < numberOfForwardAlignments() ) ;
   const int ieta ( i < numberOfForwardAlignments()/2 ? -29 : 29 ) ;
   const int iphi ( 1 + (4*i)%72 ) ;
   return HcalDetId( HcalForward, ieta, iphi, 1 ) ;
}

DetId 
HcalGeometry::detIdFromOuterAlignmentIndex( unsigned int i )
{
   assert( i < numberOfOuterAlignments() ) ;
   const int ring ( i/12 ) ;
   const int ieta ( 0 == ring ? -11 :
                    1 == ring ? -5  :
                    2 == ring ?  1  :
                    3 == ring ?  5  : 11 ) ;
   const int iphi ( 1 + ( i - ring*12 )*6 ) ;
   return HcalDetId( HcalOuter, ieta, iphi, 4 ) ;
}

DetId 
HcalGeometry::detIdFromLocalAlignmentIndex( unsigned int i )
{
   assert( i < numberOfAlignments() ) ;

   const unsigned int nB ( numberOfBarrelAlignments()  ) ;
   const unsigned int nE ( numberOfEndcapAlignments()  ) ;
   const unsigned int nF ( numberOfForwardAlignments() ) ;
//   const unsigned int nO ( numberOfOuterAlignments()   ) ;

   return (  i < nB       ? detIdFromBarrelAlignmentIndex( i ) :
             i < nB+nE    ? detIdFromEndcapAlignmentIndex( i - nB ) :
             i < nB+nE+nF ? detIdFromForwardAlignmentIndex( i - nB - nE ) :
             detIdFromOuterAlignmentIndex( i - nB - nE - nF ) ) ;
}

unsigned int 
HcalGeometry::alignmentBarEndForIndexLocal(    const DetId& id ,
                                               unsigned int nD   )
{
   const HcalDetId hid ( id ) ;
   const unsigned int iphi ( hid.iphi() ) ;
   const int ieta ( hid.ieta() ) ;
   const unsigned int index ( ( 0 < ieta ? nD/2 : 0 ) + ( iphi + 1 )%72/4 ) ;
   assert( index < nD ) ;
   return index ;
}

unsigned int 
HcalGeometry::alignmentBarrelIndexLocal(    const DetId& id )
{
  return alignmentBarEndForIndexLocal( id, numberOfBarrelAlignments() ) ;
}
unsigned int 
HcalGeometry::alignmentEndcapIndexLocal(    const DetId& id )
{
   return alignmentBarEndForIndexLocal( id, numberOfEndcapAlignments() ) ;
}

unsigned int 
HcalGeometry::alignmentForwardIndexLocal(   const DetId& id )
{
   return alignmentBarEndForIndexLocal( id, numberOfForwardAlignments() ) ;
}

unsigned int 
HcalGeometry::alignmentOuterIndexLocal(     const DetId& id )
{
   const HcalDetId hid ( id ) ;
   const int ieta ( hid.ieta() ) ;
   const int iphi ( hid.iphi() ) ;
   const int ring ( ieta < -10 ? 0 :
                    ( ieta < -4 ? 1 :
                      ( ieta < 5 ? 2 :
                        ( ieta < 11 ? 3 : 4 ) ) ) ) ;

   const unsigned int index ( 12*ring + ( iphi - 1 )/6 ) ;
   assert( index < numberOfOuterAlignments() ) ;
   return index ;
}

unsigned int
HcalGeometry::alignmentTransformIndexLocal( const DetId& id )
{
   assert(id.det() == DetId::Hcal) ;

   const HcalDetId hid ( id ) ;
   bool isHB = (hid.subdet() == HcalBarrel);
   bool isHE = (hid.subdet() == HcalEndcap);
   bool isHF = (hid.subdet() == HcalForward);
   // bool isHO = (hid.subdet() == HcalOuter);

   const unsigned int nB ( numberOfBarrelAlignments()  ) ;
   const unsigned int nE ( numberOfEndcapAlignments()  ) ;
   const unsigned int nF ( numberOfForwardAlignments() ) ;
   // const unsigned int nO ( numberOfOuterAlignments()   ) ;

   const unsigned int index ( 
      isHB ? alignmentBarrelIndexLocal(id) :
      isHE ? alignmentEndcapIndexLocal(id) + nB :
      isHF ? alignmentForwardIndexLocal( id ) + nB + nE :
      alignmentOuterIndexLocal(id) + nB + nE + nF
                              );

   assert( index < numberOfAlignments() ) ;
   return index ;
}

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

void
HcalGeometry::localCorners( Pt3DVec&        lc  ,
                            const CCGFloat* pv  ,
                            unsigned int    i   ,
                            Pt3D&           ref  )
{
  HcalDetId hid=HcalDetId(theTopology.denseId2detId(i));

   if( hid.subdet() == HcalForward )
   {
      IdealZPrism::localCorners( lc, pv, ref ) ;
   }
   else
   {
      IdealObliquePrism::localCorners( lc, pv, ref ) ;
   }
}

void
HcalGeometry::newCell( const GlobalPoint& f1 ,
                       const GlobalPoint& f2 ,
                       const GlobalPoint& f3 ,
                       const CCGFloat*    parm ,
                       const DetId&       detId   ) {

  assert (detId.det()==DetId::Hcal);
    
  const HcalDetId hid ( detId ) ;
  unsigned int din=theTopology.detId2denseId(detId);

  edm::LogInfo("HcalGeometry") << " newCell subdet "
            << detId.subdetId() << ", raw ID " 
            << detId.rawId() << ", hid " << hid << ", din " 
            << din << ", index ";
  
  if( hid.subdet()==HcalBarrel) {
    m_hbCellVec[ din ] = IdealObliquePrism( f1, cornersMgr(), parm ) ;
  } else if( hid.subdet()==HcalEndcap ) {
    const unsigned int index ( din - m_hbCellVec.size() ) ;
    m_heCellVec[ index ] = IdealObliquePrism( f1, cornersMgr(), parm ) ;
  } else if( hid.subdet()==HcalOuter ) {
    const unsigned int index ( din 
                               - m_hbCellVec.size()
                               - m_heCellVec.size() ) ;
    m_hoCellVec[ index ] = IdealObliquePrism( f1, cornersMgr(), parm ) ;
  } else {
    const unsigned int index ( din 
                               - m_hbCellVec.size()
                               - m_heCellVec.size()
                               - m_hoCellVec.size() ) ;
    m_hfCellVec[ index ] = IdealZPrism( f1, cornersMgr(), parm ) ;
  }

  m_validIds.push_back( detId ) ; 
  m_dins.push_back( din );
}

const CaloCellGeometry* 
HcalGeometry::cellGeomPtr( unsigned int din ) const
{
   const CaloCellGeometry* cell ( 0 ) ;
   if( m_hbCellVec.size() > din )
   {
      cell = &m_hbCellVec[ din ] ;
   }
   else
   {
      if( m_hbCellVec.size() +
          m_heCellVec.size() > din )
      {
         const unsigned int index ( din - m_hbCellVec.size() ) ;
         cell = &m_heCellVec[ index ] ;
      }
      else
      {
         if( m_hbCellVec.size() +
             m_heCellVec.size() +
             m_hoCellVec.size() > din )
         {
            const unsigned int index ( din 
                                       - m_hbCellVec.size() 
                                       - m_heCellVec.size() ) ;
            cell = &m_hoCellVec[ index ] ;
         }
         else
         {
            if( m_hbCellVec.size() +
                m_heCellVec.size() +
                m_hoCellVec.size() +
                m_hfCellVec.size() > din )
            {
               const unsigned int index ( din 
                                          - m_hbCellVec.size() 
                                          - m_heCellVec.size() 
                                          - m_hoCellVec.size() ) ;
               cell = &m_hfCellVec[ index ] ;
            }
         }
      }
   }
   
   return (( 0 == cell || 0 == cell->param()) ? 0 : cell ) ;
}

void
HcalGeometry::getSummary( CaloSubdetectorGeometry::TrVec&  tVec,
                          CaloSubdetectorGeometry::IVec&   iVec,
                          CaloSubdetectorGeometry::DimVec& dVec,
                          CaloSubdetectorGeometry::IVec& dinsVec ) const
{
   tVec.reserve( theTopology.ncells()*numberOfTransformParms() ) ;
   iVec.reserve( numberOfShapes()==1 ? 1 : theTopology.ncells() ) ;
   dVec.reserve( numberOfShapes()*numberOfParametersPerShape() ) ;
   dinsVec.reserve( m_dins.size() );
   
   for( ParVecVec::const_iterator ivv ( parVecVec().begin() ) ; ivv != parVecVec().end() ; ++ivv )
   {
      const ParVec& pv ( *ivv ) ;
      for( ParVec::const_iterator iv ( pv.begin() ) ; iv != pv.end() ; ++iv )
      {
         dVec.push_back( *iv ) ;
      }
   }

   for( unsigned int i ( 0 ) ; i < theTopology.ncells() ; ++i )
   {
      Tr3D tr ;
      const CaloCellGeometry* ptr ( cellGeomPtr( i ) ) ;
      dinsVec.push_back( i );
      
      assert( 0 != ptr ) ;
      ptr->getTransform( tr, ( Pt3DVec* ) 0 ) ;

      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.push_back( tt.x() ) ;
      tVec.push_back( tt.y() ) ;
      tVec.push_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.push_back( ea.Phi() ) ;
         tVec.push_back( ea.Theta() ) ;
         tVec.push_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_dins.size() ) ; 
      if( iVec.size() < nn ) iVec.push_back( ishape ) ;
   }
}