Geo2G4Builder.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
*/
 
#include "Geo2G4Builder.h"
#include "Geo2G4OpticalSurfaceFactory.h"
 
#include "GeoModelKernel/GeoVDetectorManager.h"
#include "GeoModelKernel/GeoShapeUnion.h"
#include "GeoModelKernel/GeoShapeShift.h"
#include "GeoModelKernel/GeoMaterial.h"
#include "GeoModelKernel/GeoLogVol.h"
#include "GeoModelKernel/GeoPhysVol.h"
#include "GeoModelKernel/GeoTransform.h"
 
#include "GeoModelUtilities/GeoModelExperiment.h"
#include "GeoModelInterfaces/StoredMaterialManager.h"
 
#include "AthenaBaseComps/AthMsgStreamMacros.h"
 
#include "StoreGate/StoreGateSvc.h"
 
#include "G4ReflectionFactory.hh"
#include "G4LogicalBorderSurface.hh"
#include "G4OpticalSurface.hh"
#include "GeoModelKernel/throwExcept.h"
 
#include <map>
#include <iostream>
#include <utility>
 
Geo2G4Builder::Geo2G4Builder(const std::string& detectorName)
  : AthMessaging("Geo2G4Builder")
  , m_detectorName(detectorName)
  , m_motherTransform(GeoTrf::Transform3D::Identity()) {
  if (m_pDetStore.retrieve().isFailure()) {
    THROW_EXCEPTION("ERROR: Geo2G4Builder for detector "<< detectorName << " could not access the detector store.");
  }
 
  StatusCode sc = m_pDetStore->retrieve( m_theExpt, "ATLAS" );
  if(sc.isFailure()){
    THROW_EXCEPTION(detectorName<<" could not get GeoModelExperiment");
  }
 
  const GeoVDetectorManager *theManager = m_theExpt->getManager(detectorName);
  if (!theManager) {
      THROW_EXCEPTION("Failed to retrieve manager "<<detectorName);
  }
  ATH_MSG_DEBUG("Recieve detctor manager "<<detectorName);
  for(unsigned int i=0; i<theManager->getNumTreeTops(); ++i) {
    m_treeTops.push_back(theManager->getTreeTop(i));
  }
 
  ATH_MSG_DEBUG("Found detector: top volume(s)");
  for(unsigned int i=0; i<m_treeTops.size();++i) {
      ATH_MSG_DEBUG( "   Tree Top " << i << " " << m_treeTops[i]->getLogVol()->getName() );
  }
 
  if(m_treeTops.size()>1) {
      // -------- -------- MATERIAL MANAGER -------- ----------
      StoredMaterialManager* theMaterialManager = m_pDetStore->tryRetrieve<StoredMaterialManager>("MATERIALS");
      if(theMaterialManager) {
              m_matAir = theMaterialManager->getMaterial("std::Air");
      } else {
        m_matAir = m_treeTops[0]->getLogVol()->getMaterial();
      }
  }
 
  if (m_g2gSvc.retrieve().isFailure()) {
    THROW_EXCEPTION("Failed to retrieve manager Geo2G4Svc");
  }
 
  m_theBuilder = m_g2gSvc->GetDefaultBuilder();
  if(m_theBuilder) {
    ATH_MSG_INFO("Set volume builder ---> "<< m_theBuilder->GetKey());
  } else {
    ATH_MSG_WARNING("0 pointer to volume builder."
                    <<"\n Use 'DefaultBuilder' property of Geo2G4Svc or"
                    <<"\n 'GetVolumeBuilder' method of Geo2G4Builder");
  }
}
 
G4LogicalVolume* Geo2G4Builder::BuildTree()
{
  ATH_MSG_DEBUG("Entering Geo2G4Builder::BuildTree()...");
  G4LogicalVolume* result = nullptr;
  std::unique_ptr<OpticalVolumesMap> optical_volumes{};
  const GeoBorderSurfaceContainer* surface_container = nullptr;
 
  // Check whether we have to deal with optical surfaces
  if(m_pDetStore->contains<GeoBorderSurfaceContainer>(m_detectorName)) {
    StatusCode sc = m_pDetStore->retrieve(surface_container,m_detectorName);
    if(sc.isSuccess() && surface_container->size()>0) {
      optical_volumes = std::make_unique<OpticalVolumesMap>();
    }
  }
 
  if(m_theBuilder) {
    if(m_treeTops.size()==1) {
      m_motherTransform = m_treeTops[0]->getX();
      result = m_theBuilder->Build(m_treeTops[0], optical_volumes.get());
    } else if (!m_treeTops.empty()) {
      // Create temporary GeoModel physical volume
      // The shape is composed by TreeTop shapes + their transforms
      const GeoShape& shFirst = (*(m_treeTops[0]->getLogVol()->getShape()))<<(m_treeTops[0]->getX());
      GeoIntrusivePtr<const GeoShape> shResult{&shFirst};
 
      for(unsigned int i=1; i<m_treeTops.size(); i++) {
        GeoIntrusivePtr<const GeoShape> booleanOp{&shResult->add((*(m_treeTops[i]->getLogVol()->getShape()))<<(m_treeTops[i]->getX()))};
        shResult = booleanOp;
      }
 
      auto lvEnvelope = make_intrusive<GeoLogVol>(m_detectorName, shResult, m_matAir);
      auto pvEnvelope = make_intrusive<GeoPhysVol>(lvEnvelope);
      m_theExpt->addTmpVolume(pvEnvelope);
      result = m_theBuilder->Build(pvEnvelope);
 
      // Get pointer to the World
      PVConstLink world = m_treeTops[0]->getParent();
 
      // Add all tree tops to the result
      for(unsigned int i=0; i<m_treeTops.size(); i++) {
        // Current Tree Top and its index
        PVConstLink pv = m_treeTops[i];
        std::optional<unsigned int> childIndx = world->indexOf(pv);
 
        // Tree Top transformation
        G4Transform3D theG4Position(Amg::EigenTransformToCLHEP(world->getXToChildVol(*childIndx)));
 
        // Copy number
        int id = 16969;
        std::optional<int> Qint = world->getIdOfChildVol(*childIndx);
        if(Qint) id = *Qint;
 
        // PV Tree Top name
        std::string nameTT =  world->getNameOfChildVol(*childIndx);
        if (nameTT == "ANON") nameTT = pv->getLogVol()->getName();
 
 
        G4LogicalVolume* g4LV = m_theBuilder->Build(pv, optical_volumes.get());
        G4ReflectionFactory::Instance()->Place(theG4Position,
                                               nameTT,
                                               g4LV,
                                               result,
                                               false,
                                               id);
      }
    }
  }
 
  // build optical surfaces if necessary
  if(optical_volumes) {
    if(optical_volumes->size()>0) {
      BuildOpticalSurfaces(surface_container,optical_volumes.get());
    } else {
      ATH_MSG_WARNING("Optical volumes apparently requested, but none found!  Deleting temps");
    }
  }
 
  return result;
}
 
 
VolumeBuilder*  Geo2G4Builder::GetVolumeBuilder(std::string bname)
{
  m_theBuilder = m_g2gSvc->GetVolumeBuilder(std::move(bname));
  return m_theBuilder;
}
 
void Geo2G4Builder::BuildOpticalSurfaces(const GeoBorderSurfaceContainer* surface_container,
                                         const OpticalVolumesMap* optical_volumes)
{
  Geo2G4OpticalSurfaceFactory surfaceFactory;
 
  for (const GeoBorderSurface& border_surface : *surface_container)
    {
      // Build Optical Surface
      G4OpticalSurface* g4OptSurface = surfaceFactory.Build(border_surface.getOptSurface());
 
      G4VPhysicalVolume* g4PV1 = 0;
      G4VPhysicalVolume* g4PV2 = 0;
      OpticalVolumesMap::const_iterator volIt;
 
      // First physical volume
      volIt = optical_volumes->find(border_surface.getPV1());
      if(volIt == optical_volumes->end())
        {
          ATH_MSG_WARNING("Unable to find " << border_surface.getPV1()->getLogVol()->getName() << " in Optical Volumes map");
          continue;
        }
      g4PV1 = volIt.operator->()->second;
 
      // Second physical volume
      volIt = optical_volumes->find(border_surface.getPV2());
      if(volIt == optical_volumes->end())
        {
          ATH_MSG_WARNING("Unable to find " << border_surface.getPV2()->getLogVol()->getName() << " in Optical Volumes map");
          continue;
        }
      g4PV2 = volIt.operator->()->second;
 
      // G4LogicalBorderSurface
      G4LogicalBorderSurface* g4BorderSurface __attribute__((unused)) = new G4LogicalBorderSurface(border_surface.getName(),g4PV1,g4PV2,g4OptSurface);
    }
}