Geant4DetectorGeometryConstruction.cpp

//==========================================================================
//  AIDA Detector description implementation 
//--------------------------------------------------------------------------
// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
// All rights reserved.
//
// For the licensing terms see $DD4hepINSTALL/LICENSE.
// For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
//
//  \author Markus Frank
//  \date   2015-11-09
//
//==========================================================================

// Framework include files
#include <DDG4/Geant4DetectorConstruction.h>

/// Namespace for the AIDA detector description toolkit
namespace dd4hep {

  /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
  namespace sim {

    /// Class to create Geant4 detector geometry from TGeo representation in memory
    /**
     *  On demand (ie. when calling "Construct") the dd4hep geometry is converted
     *  to Geant4 with all volumes, assemblies, shapes, materials etc.
     *  The actuak work is performed by the Geant4Converter class called by this method.
     *
     *  \author  M.Frank
     *  \version 1.0
     *  \ingroup DD4HEP_SIMULATION
     */
    class Geant4DetectorGeometryConstruction : public Geant4DetectorConstruction   {
      /// Property: Dump geometry hierarchy if not NULL. Flags can steer actions. 
      unsigned long m_dumpHierarchy = 0;
      /// Property: Flag to debug materials during conversion mechanism
      bool m_debugMaterials         = false;
      /// Property: Flag to debug elements during conversion mechanism
      bool m_debugElements          = false;
      /// Property: Flag to debug shapes during conversion mechanism
      bool m_debugShapes            = false;
      /// Property: Flag to debug volumes during conversion mechanism
      bool m_debugVolumes           = false;
      /// Property: Flag to debug placements during conversion mechanism
      bool m_debugPlacements        = false;
      /// Property: Flag to debug reflections during conversion mechanism
      bool m_debugReflections       = false;
      /// Property: Flag to debug regions during conversion mechanism
      bool m_debugRegions           = false;
      /// Property: Flag to debug regions during conversion mechanism
      bool m_debugSurfaces          = false;

      /// Property: Flag to dump all placements after the conversion procedure
      bool m_printPlacements        = false;
      /// Property: Flag to dump all sensitives after the conversion procedure
      bool m_printSensitives        = false;

      /// Property: Printout level of info object
      int  m_geoInfoPrintLevel;
      /// Property: G4 GDML dump file name (default: empty. If non empty, dump)
      std::string m_dumpGDML;

      /// Write GDML file
      int writeGDML(const char* gdml_output);
      /// Print geant4 volume 
      int printVolumeObj(const char* vol_path, PlacedVolume pv, int flg);
      /// Print volume tree with attributes
      int printVolumeTree(const char* vol_path);
      /// Print volume tree WITHOUT attributes
      int printVolTree(const char* vol_path);
      /// Print geant4 volume tree
      int printG4Tree(const char* vol_path);
      /// Print geant4 volume
      int printVolume(const char* vol_path);
      /// Check geant4 volume
      int checkVolume(const char* vol_path);
      /// Print geant4 material
      int printMaterial(const char* mat_name);

      std::pair<std::string, PlacedVolume> resolve_path(const char* vol_path)   const;
      void printG4(const std::string& prefix, const G4VPhysicalVolume* g4pv)  const;

    public:
      /// Initializing constructor for DDG4
      Geant4DetectorGeometryConstruction(Geant4Context* ctxt, const std::string& nam);
      /// Default destructor
      virtual ~Geant4DetectorGeometryConstruction();
      /// Geometry construction callback. Called at "Construct()"
      void constructGeo(Geant4DetectorConstructionContext* ctxt)  override;
      /// Install command control messenger to write GDML file from command prompt.
      virtual void installCommandMessenger()   override;
    };
  }    // End namespace sim
}      // End namespace dd4hep


// Framework include files
#include <DD4hep/InstanceCount.h>
#include <DD4hep/DetectorTools.h>
#include <DD4hep/DD4hepUnits.h>
#include <DD4hep/Printout.h>
#include <DD4hep/Detector.h>

#include <DDG4/Geant4HierarchyDump.h>
#include <DDG4/Geant4UIMessenger.h>
#include <DDG4/Geant4Converter.h>
#include <DDG4/Geant4Kernel.h>
#include <DDG4/Factories.h>

#include <TGeoScaledShape.h>

// Geant4 include files
#include <G4LogicalVolume.hh>
#include <G4PVPlacement.hh>
#include <G4Material.hh>
#include <G4Version.hh>
#include <G4VSolid.hh>
#include "CLHEP/Units/SystemOfUnits.h"

//#ifdef GEANT4_HAS_GDML
#include <G4GDMLParser.hh>
//#endif

#include <cmath>

using namespace std;
using namespace dd4hep;
using namespace dd4hep::sim;
DECLARE_GEANT4ACTION(Geant4DetectorGeometryConstruction)

/// Initializing constructor for other clients
Geant4DetectorGeometryConstruction::Geant4DetectorGeometryConstruction(Geant4Context* ctxt, const string& nam)
: Geant4DetectorConstruction(ctxt,nam)
{
  declareProperty("DebugMaterials",    m_debugMaterials);
  declareProperty("DebugElements",     m_debugElements);
  declareProperty("DebugShapes",       m_debugShapes);
  declareProperty("DebugVolumes",      m_debugVolumes);
  declareProperty("DebugPlacements",   m_debugPlacements);
  declareProperty("DebugReflections",  m_debugReflections);
  declareProperty("DebugRegions",      m_debugRegions);
  declareProperty("DebugSurfaces",     m_debugSurfaces);

  declareProperty("PrintPlacements",   m_printPlacements);
  declareProperty("PrintSensitives",   m_printSensitives);
  declareProperty("GeoInfoPrintLevel", m_geoInfoPrintLevel = DEBUG);

  declareProperty("DumpHierarchy",     m_dumpHierarchy);
  declareProperty("DumpGDML",          m_dumpGDML="");
  InstanceCount::increment(this);
}

/// Default destructor
Geant4DetectorGeometryConstruction::~Geant4DetectorGeometryConstruction() {
  InstanceCount::decrement(this);
}

/// Geometry construction callback. Called at "Construct()"
void Geant4DetectorGeometryConstruction::constructGeo(Geant4DetectorConstructionContext* ctxt)   {
  Geant4Mapping&  g4map = Geant4Mapping::instance();
  DetElement      world = ctxt->description.world();
  Geant4Converter conv(ctxt->description, outputLevel());
  conv.debugMaterials   = m_debugMaterials;
  conv.debugElements    = m_debugElements;
  conv.debugShapes      = m_debugShapes;
  conv.debugVolumes     = m_debugVolumes;
  conv.debugRegions     = m_debugRegions;
  conv.debugSurfaces    = m_debugSurfaces;
  conv.debugPlacements  = m_debugPlacements;
  conv.debugReflections = m_debugReflections;

  ctxt->geometry = conv.create(world).detach();
  ctxt->geometry->printLevel = outputLevel();
  g4map.attach(ctxt->geometry);
  G4VPhysicalVolume* w = ctxt->geometry->world();
  // Save away the reference to the world volume
  context()->kernel().setWorld(w);
  // Create Geant4 volume manager only if not yet available
  g4map.volumeManager();
  if ( m_dumpHierarchy != 0 )   {
    Geant4HierarchyDump dmp(ctxt->description, m_dumpHierarchy);
    dmp.dump("",w);
  }
  ctxt->world = w;
  if ( !m_dumpGDML.empty() ) writeGDML(m_dumpGDML.c_str());
  else if ( ::getenv("DUMP_GDML") ) writeGDML(::getenv("DUMP_GDML"));
  enableUI();
}

pair<string, PlacedVolume>
Geant4DetectorGeometryConstruction::resolve_path(const char* vol_path)  const {
  string       p   = vol_path;
  Detector&    det = context()->kernel().detectorDescription();
  PlacedVolume top = det.world().placement();
  PlacedVolume pv  = detail::tools::findNode(top, p);
  if ( !pv.isValid() )    {
    DetElement de = detail::tools::findElement(det, p);
    if ( de.isValid() )  {
      pv = de.placement();
      p = detail::tools::placementPath(de);
    }
  }
  return make_pair(p,pv);
}

/// Print geant4 material
int Geant4DetectorGeometryConstruction::printMaterial(const char* mat_name)  {
  if ( mat_name )   {
    auto& g4map = Geant4Mapping::instance().data().g4Materials;
    for ( auto it = g4map.begin(); it != g4map.end(); ++it )  {
      const auto* mat = (*it).second;
      if ( mat->GetName() == mat_name )   {
	stringstream output;
        const auto* ion = mat->GetIonisation();
        printP2("+++  Dump of GEANT4 material: %s", mat_name);
        output << mat;
        if ( ion )   {
          output << "          MEE:  ";
          output << setprecision(12);
          output << ion->GetMeanExcitationEnergy()/CLHEP::eV;
          output << " [eV]";
        }
        else
          output << "          MEE: UNKNOWN";
	always("+++ printMaterial: \n%s\n", output.str().c_str());
        return 1;
      }
    }
    warning("+++ printMaterial: FAILED to find the material %s", mat_name);
  }
  warning("+++ printMaterial: Property materialName not set!");
  return 0;
}

/// Print geant4 volume
int Geant4DetectorGeometryConstruction::printVolumeObj(const char* vol_path, PlacedVolume pv, int flg)   {
  if ( pv.isValid() )   {
    const G4LogicalVolume* vol = 0;
    auto& g4map = Geant4Mapping::instance().data();
    auto pit = g4map.g4Placements.find(pv.ptr());
    auto vit = g4map.g4Volumes.find(pv.volume());
    warning("+++ printVolume: %s", vol_path);
    if ( vit != g4map.g4Volumes.end() )   {
      vol = (*vit).second;
      auto* sol = vol->GetSolid();
      const auto* mat = vol->GetMaterial();
      const auto* ion = mat->GetIonisation();
      Solid sh  = pv.volume().solid();
      if ( flg )  {
	printP2("++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++");
	printP2("+++  Dump of GEANT4 solid: %s", vol_path);
      }
      stringstream output;
      if ( flg )   {
	output << mat;
	if ( ion )   {
	  output << "          MEE:  ";
	  output << setprecision(12);
	  output << ion->GetMeanExcitationEnergy()/CLHEP::eV;
	  output << " [eV]";
	}
	else
	  output << "          MEE: UNKNOWN";
      }
      if ( flg )   {
	output << endl << *sol;
	printP2("%s", output.str().c_str());
	printP2("+++  Dump of ROOT   solid: %s", vol_path);
	sh->InspectShape();
	if ( sh->IsA() == TGeoScaledShape::Class() )    {
	  TGeoScaledShape* scaled = (TGeoScaledShape*)sh.ptr();
	  const Double_t* scale = scaled->GetScale()->GetScale();
	  double dot = scale[0]*scale[1]*scale[2];
	  printP2("+++ TGeoScaledShape: %8.3g  %8.3g  %8.3g  [%s]", scale[0], scale[1], scale[2],
		  dot > 0e0 ? "RIGHT handed" : "LEFT handed");
	}
	else if ( pit != g4map.g4Placements.end() )   {
	  const G4VPhysicalVolume* pl = (*pit).second;
	  const G4RotationMatrix* rot = pl->GetRotation();
	  const G4ThreeVector& tr = pl->GetTranslation();
	  G4Transform3D transform(rot ? *rot : G4RotationMatrix(), tr);
	  HepGeom::Scale3D  sc;
	  HepGeom::Rotate3D rr;
	  G4Translate3D     tt;
	  transform.getDecomposition(sc,rr,tt);
	  double dot = sc(0,0)*sc(1,1)*sc(2,2);
	  printP2("+++ TGeoShape:       %8.3g  %8.3g  %8.3g  [%s]", sc(0,0), sc(1,1), sc(2,2),
		  dot > 0e0 ? "RIGHT handed" : "LEFT handed");
	}
	const TGeoMatrix* matrix = pv->GetMatrix();
	printP2("+++ TGeoMatrix:      %s",
		matrix->TestBit(TGeoMatrix::kGeoReflection) ? "LEFT handed" : "RIGHT handed");        
	printP2("+++ Shape: %s  cubic volume: %8.3g mm^3  area: %8.3g mm^2",
		sol->GetName().c_str(), sol->GetCubicVolume(), sol->GetSurfaceArea());
      }
      return 1;
    }
    else   {
      auto ai = g4map.g4AssemblyVolumes.find(pv.ptr());
      if ( ai != g4map.g4AssemblyVolumes.end() )   {
        Volume v = pv.volume();
        warning("+++ printVolume: volume %s is an assembly...need to resolve imprint",vol_path);
        for(Int_t i=0; i < v->GetNdaughters(); ++i)   {
          TGeoNode*   dau_nod = v->GetNode(i);
          string p = vol_path + string("/") + dau_nod->GetName();
          printVolumeObj(p.c_str(), dau_nod, flg);
        }
        return 0;
      }
    }
    warning("+++ printVolume: FAILED to find the volume %s in geant4 mapping...",vol_path);
    return 0;
  }
  warning("+++ printVolume: FAILED to dump invalid volume",vol_path);
  return 0;
}

/// Print geant4 volume
int Geant4DetectorGeometryConstruction::printVolume(const char* vol_path)  {
  if ( vol_path )   {
    auto [p, pv] = resolve_path(vol_path);
    if ( pv.isValid() )    {
      return printVolumeObj(vol_path, pv, ~0x0);
    }
  }
  warning("+++ printVolume: Property VolumePath not set. [Ignored]");
  return 0;
}

/// Print geant4 volume
int Geant4DetectorGeometryConstruction::printVolumeTree(const char* vol_path)  {
  if ( vol_path )   {
    auto [p, pv] = resolve_path(vol_path);
    if ( pv.isValid() )    {
      if ( printVolumeObj(p.c_str(), pv, ~0x0) )     {
	TGeoVolume* vol = pv->GetVolume();
	for(Int_t i=0; i < vol->GetNdaughters(); ++i)   {
	  PlacedVolume dau_pv(vol->GetNode(i));
	  string path = (p + "/") + dau_pv.name();
	  if ( printVolumeTree(path.c_str()) )     {
	  }
	}
      }
      return 1;
    }
  }
  warning("+++ printVolume: Could not access Volume/DetElement '%s'", vol_path ? vol_path : "UNKNOWN");
  return 0;
}

int Geant4DetectorGeometryConstruction::printVolTree(const char* vol_path)  {
  if ( vol_path )   {
    auto [p, pv] = resolve_path(vol_path);
    if ( pv.isValid() )    {
      if ( printVolumeObj(p.c_str(), pv, 0) )     {
	TGeoVolume* vol = pv->GetVolume();
	for(Int_t i=0; i < vol->GetNdaughters(); ++i)   {
	  PlacedVolume dau_pv(vol->GetNode(i));
	  string path = (p + "/") + dau_pv.name();
	  if ( printVolTree(path.c_str()) )     {
	  }
	}
      }
      return 1;
    }
  }
  warning("+++ printVolume: Could not access Volume/DetElement '%s'", vol_path ? vol_path : "UNKNOWN");
  return 0;
}

/// Check geant4 volume
int Geant4DetectorGeometryConstruction::checkVolume(const char* vol_path)  {
  if ( vol_path )   {
    auto [p, pv] = resolve_path(vol_path);
    if ( pv.isValid() )   {
      auto& g4map = Geant4Mapping::instance().data().g4Volumes;
      auto it = g4map.find(pv.volume());
      if ( it != g4map.end() )   {
        const G4LogicalVolume* vol = (*it).second;
        auto* g4_sol = vol->GetSolid();
        Box   rt_sol = pv.volume().solid();
        printP2("Geant4 Shape: %s  cubic volume: %8.3g mm^3  area: %8.3g mm^2",
                g4_sol->GetName().c_str(), g4_sol->GetCubicVolume(), g4_sol->GetSurfaceArea());
#if G4VERSION_NUMBER>=1040
        G4ThreeVector pMin, pMax;
        double conv = (dd4hep::centimeter/CLHEP::centimeter)/2.0;
        g4_sol->BoundingLimits(pMin,pMax);
        printP2("Geant4 Bounding box extends:    %8.3g  %8.3g %8.3g",
                (pMax.x()-pMin.x())*conv, (pMax.y()-pMin.y())*conv, (pMax.z()-pMin.z())*conv);
#endif
        printP2("ROOT   Bounding box dimensions: %8.3g  %8.3g %8.3g",
                rt_sol->GetDX(), rt_sol->GetDY(), rt_sol->GetDZ());
        
        return 1;
      }
    }
    warning("+++ checkVolume: FAILED to find the volume %s from the top volume",vol_path);
  }
  warning("+++ checkVolume: Property VolumePath not set. [Ignored]");
  return 0;
}

/// Write GDML file
int Geant4DetectorGeometryConstruction::writeGDML(const char* output)  {
  G4VPhysicalVolume* w  = context()->world();
  if ( output && ::strlen(output) > 0 && output != m_dumpGDML.c_str() )
    m_dumpGDML = output;

  //#ifdef GEANT4_HAS_GDML
  if ( !m_dumpGDML.empty() ) {
    G4GDMLParser parser;
    parser.Write(m_dumpGDML.c_str(), w);
    info("+++ writeGDML: Wrote GDML file: %s", m_dumpGDML.c_str());
    return 1;
  }
  else {
    const char* gdml_dmp = ::getenv("DUMP_GDML");
    if ( gdml_dmp )    {
      G4GDMLParser parser;
      parser.Write(gdml_dmp, w);
      info("+++ writeGDML: Wrote GDML file: %s", gdml_dmp);
      return 1;
    }
  }
  warning("+++ writeGDML: Neither property DumpGDML nor environment DUMP_GDML set. No file written!");
  //#endif
  return 0;
}

void Geant4DetectorGeometryConstruction::printG4(const string& prefix, const G4VPhysicalVolume* g4pv)    const   {
  string path = prefix + "/";
  printP2(  "+++  GEANT4 volume: %s", prefix.c_str());
  auto* g4v = g4pv->GetLogicalVolume();
  for(size_t i=0, n=g4v->GetNoDaughters(); i<n; ++i)    {
    auto* dau = g4v->GetDaughter(i);
    printG4(path + dau->GetName(), dau);
  }
}

int Geant4DetectorGeometryConstruction::printG4Tree(const char* vol_path)  {
  if ( vol_path )   {
    auto [p, pv] = resolve_path(vol_path);
    if ( pv.isValid() )    {
      auto& g4map = Geant4Mapping::instance().data().g4Placements;
      auto it = g4map.find(pv);
      if ( it != g4map.end() )   {
	printG4(p, (*it).second);
      }
      return 1;
    }
  }
  warning("+++ printVolume: Could not access Volume/DetElement '%s'", vol_path ? vol_path : "UNKNOWN");
  return 0;
}

/// Install command control messenger to write GDML file from command prompt.
void Geant4DetectorGeometryConstruction::installCommandMessenger()   {
  this->Geant4DetectorConstruction::installCommandMessenger();
  m_control->addCall("writeGDML", "GDML: write geometry to file: '"+m_dumpGDML+
                     "' [uses argument - or - property DumpGDML]",
                     Callback(this).make(&Geant4DetectorGeometryConstruction::writeGDML),1);
  m_control->addCall("printVolume", "Print Geant4 volume properties [uses argument]",
                     Callback(this).make(&Geant4DetectorGeometryConstruction::printVolume),1);
  m_control->addCall("printTree",   "Print volume tree WITHOUT properties [uses argument]",
                     Callback(this).make(&Geant4DetectorGeometryConstruction::printVolTree),1);
  m_control->addCall("printG4Tree",   "Print Geant4 volume tree [uses argument]",
                     Callback(this).make(&Geant4DetectorGeometryConstruction::printG4Tree),1);
  m_control->addCall("printVolumeTree", "Print volume tree with properties [uses argument]",
                     Callback(this).make(&Geant4DetectorGeometryConstruction::printVolumeTree),1);
  m_control->addCall("checkVolume", "Check Geant4 volume properties [uses argument]",
                     Callback(this).make(&Geant4DetectorGeometryConstruction::checkVolume),1);
  m_control->addCall("printMaterial", "Print Geant4 material properties [uses argument]",
                     Callback(this).make(&Geant4DetectorGeometryConstruction::printMaterial),1);
}