SCT_OverlapDescriptor.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
*/
 
// SCT_OverlapDescriptor.cxx, (c) ATLAS Detector software
 
// Amg
#include "GeoPrimitives/GeoPrimitives.h"
// InDet
#include "InDetReadoutGeometry/SiDetectorElement.h"
#include "InDetTrackingGeometry/SCT_OverlapDescriptor.h"
#include "InDetIdentifier/SCT_ID.h"
// Trk
#include "TrkSurfaces/Surface.h"
// 
#include "StoreGate/StoreGateSvc.h"
 
InDet::SCT_OverlapDescriptor::SCT_OverlapDescriptor(bool addMoreSurfaces, int eta_slices)
  : m_robustMode(true),
  m_addMoreSurfaces(addMoreSurfaces), 
  m_etaSlices(eta_slices)
{}
 
bool
InDet::SCT_OverlapDescriptor::reachableSurfaces(
  std::vector<Trk::SurfaceIntersection>& cSurfaces,
  const Trk::Surface& tsf,
  const Amg::Vector3D& pos,
  const Amg::Vector3D&) const
 
{
  // first add the target surface and the backside surface (in the if statement)
  cSurfaces.emplace_back(Trk::Intersection(pos, 0., true), &tsf);
 
  // make sure we have the correct associated element
  const auto* tmp = tsf.associatedDetectorElement();
  const InDetDD::SiDetectorElement* sElement =
    tmp->detectorType() == Trk::DetectorElemType::Silicon
      ? static_cast<const InDetDD::SiDetectorElement*>(tmp)
      : nullptr;
      
  if (sElement) {
    
    size_t newCapacity = cSurfaces.size() + 2;
    if (m_robustMode and m_addMoreSurfaces and sElement->isBarrel()) {
      // sum up the defined slices to evaluate the new capacity
      // only uses one additional slice in phi
      newCapacity += (2 + 6*m_etaSlices)*2;
    } else if (m_robustMode) {
      newCapacity += 16;
    } else {
      newCapacity += 6;
    }
    cSurfaces.reserve(newCapacity);
    
    double surfacePhi = tsf.center().phi() + M_PI;
    double positionPhi = pos.phi() + M_PI;
 
    // now get the overlap options
    addOtherSide(sElement, cSurfaces);
 
    // 8-cell-connectivity depending on track/surface geometry
    // nPhi - can be jump + or -
    const InDetDD::SiDetectorElement* nElement = nullptr;
 
    // robust mode --> return 9 (*2) surfaces
    if (m_robustMode) {
      addNextInPhiOS(sElement, cSurfaces);
      addNextInEtaOS(sElement, cSurfaces);
 
      addPrevInPhiOS(sElement, cSurfaces);
      addPrevInEtaOS(sElement, cSurfaces);
 
      nElement = sElement->nextInPhi();
      addNextInEtaOS(nElement, cSurfaces);
      addPrevInEtaOS(nElement, cSurfaces);
 
      nElement = sElement->prevInPhi();
      addNextInEtaOS(nElement, cSurfaces);
      addPrevInEtaOS(nElement, cSurfaces);
      
      if (m_addMoreSurfaces and sElement->isBarrel()) {  
        unsigned int next = 1;
        const InDetDD::SiDetectorElement* currentElement = sElement->nextInEta();
        while (currentElement and next<(unsigned int)m_etaSlices) {
          addNextInEtaOS(currentElement,cSurfaces);
          currentElement = currentElement->nextInEta();
          if (currentElement) {
            addNextInPhiOS(currentElement,cSurfaces);
            addPrevInPhiOS(currentElement,cSurfaces);
          }
          next++;
        }
        
        unsigned int prev = 1;
        currentElement = sElement->prevInEta();
        while (currentElement and prev<(unsigned int)m_etaSlices) {
          addPrevInEtaOS(currentElement,cSurfaces);
          currentElement = currentElement->prevInEta();
          if (currentElement) {
            addNextInPhiOS(currentElement,cSurfaces);
            addPrevInPhiOS(currentElement,cSurfaces);
          }
          prev++;
        }
      }      
      
    } else {
      // get the phi information
      if (surfacePhi < positionPhi) {
        addNextInPhiOS(sElement, cSurfaces);
        nElement = sElement->nextInPhi();
      } else {
        addPrevInPhiOS(sElement, cSurfaces);
        nElement = sElement->prevInPhi();
      }
      // get the eta information - also possible
      double positionEta = pos.eta();
      double surfaceEta = tsf.center().eta();
      int side = (sElement->isBarrel() || surfaceEta > 0.) ? 1 : -1;
      // check the surface / position eta values
      if (side * surfaceEta < side * positionEta) {
        // we go next in eta for both, the original and the phi jumped one
        addNextInEtaOS(sElement, cSurfaces);
        addNextInEtaOS(nElement, cSurfaces);
      } else {
        // opposite direction
        addPrevInEtaOS(sElement, cSurfaces);
        addPrevInEtaOS(nElement, cSurfaces);
      }
    }
  }
 
  return false;
}
 
bool InDet::SCT_OverlapDescriptor::dumpSurfaces(std::vector<Trk::SurfaceIntersection>& surfaces) const {
  
  if (m_sctIdHelper==nullptr) {
    ISvcLocator* svcLocator = Gaudi::svcLocator();
    
    // get DetectorStore service
    StoreGateSvc* detStore = nullptr;
    if (svcLocator->service("DetectorStore", detStore).isFailure()) {
      return false;
    }
 
    const SCT_ID* sctIdHelper = nullptr;
    if (detStore->retrieve(sctIdHelper, "SCT_ID").isFailure()) {
      return false;
    }
    m_sctIdHelper = sctIdHelper;
  }
 
  std::cout << "Dumping Surfaces for SCT with size = " << surfaces.size() << std::endl;
  for (auto & surface : surfaces) {
    Identifier hitId = (surface.object)->associatedDetectorElementIdentifier(); 
    std::cout <<  "barrel_ec " << m_sctIdHelper.load()->barrel_ec(hitId) 
              << ", layer_disk " << m_sctIdHelper.load()->layer_disk(hitId) << ", phi_module " << m_sctIdHelper.load()->phi_module(hitId) 
              << ", eta_module " << m_sctIdHelper.load()->eta_module(hitId) << ", side " << m_sctIdHelper.load()->side(hitId) <<std::endl;
  }  
  return true;
}