GGSDetectorConstruction.cpp

Go to the documentation of this file.

/*
 * GGSDetectorConstruction.cpp
 *
 *  Created on: 2010-09-29
 *     Authors: Emiliano Mocchiutti and Cecilia Pizzolotto
 */

// Geant4 headers
#include "G4GeometryManager.hh"
#include "G4LogicalVolumeStore.hh"
#include "G4PhysicalVolumeStore.hh"
#include "G4SDManager.hh"
#include "G4SolidStore.hh"
#include "G4UIcmdWithAString.hh"
#include "G4UImanager.hh"
#ifdef USE_GDML
#include "G4GDMLParser.hh"
#endif

// DD4HEP headers
#ifdef GGS_USE_DD4HEP
#include "DD4hep/Detector.h"
#include "DDG4/Geant4Converter.h"
#include "DDG4/Geant4DetectorConstruction.h"
#endif

// GGS headers
#include "geometry/GGSVGeometryConstruction.h"
#include "geometry/pluginmanagers/GGSGeoPluginManager.h"
#include "montecarlo/mccore/GGSDetectorConstruction.h"
#include "montecarlo/scoring/GGSIntHitSD.h"
#include "montecarlo/scoring/GGSMultiSensitiveDetector.h"
#include "utils/GGSNameDecoder.h"
#include "utils/GGSSmartLog.h"

// C++ headers
#include <stdexcept>

namespace {
// Messenger for build commands for sensitive detectors
class SDBuildMessenger : public G4UImessenger {
public:
  SDBuildMessenger(const G4String &path, const G4String &dsc);
  void AddSDBuildCommand(const std::string &newCommand);
  void SetNewValue(G4UIcommand *command, G4String newValue);
  void BuildSD(const std::string &sdClassName, const std::string &logVolNameAndParams);

private:
  std::vector<G4UIcmdWithAString> _sdBuildCommands;
};

G4ThreadLocal std::unique_ptr<SDBuildMessenger> sdBuildMessenger = nullptr;

SDBuildMessenger::SDBuildMessenger(const G4String &path, const G4String &dsc) : G4UImessenger(path, dsc) {}

void SDBuildMessenger::AddSDBuildCommand(const std::string &newCommand) {
  _sdBuildCommands.emplace_back(newCommand.c_str(), this);
}

void SDBuildMessenger::SetNewValue(G4UIcommand *command, G4String newValue) {
  std::string path = command->GetCommandPath();
  if (path.size() > 16 && path.substr(0, 16) == "/GGS/scoring/add") {
    // Build user actions
    auto sdClassName = path.substr(16);
    BuildSD(sdClassName, newValue);
  }
}

void SDBuildMessenger::BuildSD(const std::string &sdClassName, const std::string &logVolNameAndParams) {
  static const std::string routineName("GGSScoringManager::AddSDToLogVol");

  // Attach the sensitive detector to all the logical volumes
  G4String logVolName = logVolNameAndParams.substr(0, logVolNameAndParams.find_first_of(' '));
  auto separatorPos = logVolNameAndParams.find_first_of(' ');
  G4String params;
  if (separatorPos != std::string::npos) {
    params = logVolNameAndParams.substr(logVolNameAndParams.find_first_of(' ') + 1);
  }

  auto AddSDToLogVol = [](G4VSensitiveDetector *sd, G4LogicalVolume *logVol) {
    // Add sd to logical volume
    if (logVol->GetSensitiveDetector() == nullptr) {
      // Case 1: logical has not a sd. Add the newly built sd to it
      logVol->SetSensitiveDetector(sd);
    } else {
      GGSMultiSensitiveDetector *multiSD = dynamic_cast<GGSMultiSensitiveDetector *>(logVol->GetSensitiveDetector());
      if (multiSD) {
        // Case 2: the current sd is a multisd. Add the newly built sd to it
        multiSD->AddSensitiveDetector(sd);
      } else {
        // Case 3: the current sd is a standard sd. Replace it with a multisd and add it and the newly built sd to
        // the multisd
        multiSD = new GGSMultiSensitiveDetector(G4String(logVol->GetName()).append(".GGSMultiSD"));
        multiSD->AddSensitiveDetector(logVol->GetSensitiveDetector());
        multiSD->AddSensitiveDetector(sd);
        G4SDManager::GetSDMpointer()->AddNewDetector(multiSD);
        logVol->SetSensitiveDetector(multiSD);
      }
    }

    // Add the sd to the sdmanager if it has not already been added
    if (G4SDManager::GetSDMpointer()->GetHCtable()->GetCollectionID(sd->GetName()) == -1) {
      G4SDManager::GetSDMpointer()->AddNewDetector(sd);
    }
  };

  G4LogicalVolumeStore::const_iterator logVol = G4LogicalVolumeStore::GetInstance()->begin();
  G4VSensitiveDetector *sd = nullptr;
  for (; logVol != G4LogicalVolumeStore::GetInstance()->end(); ++logVol) {
    if ((*logVol)->GetName() == logVolName) {
      if (!sd) {
        // Build the sd
        G4String sdName = sdClassName;
        if (params != "") {
          sdName.append(".").append(params);
        }

        sd = GGSFactory<G4VSensitiveDetector, std::string>::GetInstance()
                 .CreateObject(sdClassName, logVolName + '.' + sdName)
                 .release();
        if (!sd) {
          G4ExceptionDescription ed;
          ed << "No registered builder for sensitive detector " << sdClassName << std::endl;
          G4Exception("GGSUserActionInitialization::UAIMessenger::SetNewValue", "Invalid generator",
                      G4ExceptionSeverity::FatalException, ed);
        }
      }
      // Add sd to logvol
      if (G4Threading::IsMultithreadedApplication() && G4Threading::IsMasterThread()) {
        // Add the sd to the sdmanager if it has not already been added
        if (G4SDManager::GetSDMpointer()->GetHCtable()->GetCollectionID(sd->GetName()) == -1) {
          G4SDManager::GetSDMpointer()->AddNewDetector(sd);
        }
      } else {
        AddSDToLogVol(sd, *logVol);
      }
    }
  }

  if (!sd) {
    G4ExceptionDescription ed;
    ed << routineName << ": logical volume " << logVolName << " not present. Sensitive detector " << sdClassName
       << " not added." << std::endl;
    G4Exception("GGSUserActionInitialization::UAIMessenger::SetNewValue", "Invalid generator",
                G4ExceptionSeverity::FatalException, ed);
  }
}

} // namespace

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

GGSDetectorConstruction::GGSDetectorConstruction(const G4String &library, const G4String &configDataCard)
    : _physicalWorld(nullptr), _geometry(nullptr), _library(library), _configDataCard(configDataCard), _gdxml(""),
      _messenger(this, "/GGS/"), _nToPrintForVectGeoParams{10} {
  _messenger.DeclareProperty("nToPrintForVectGeoParams", _nToPrintForVectGeoParams,
                             "Maximum number of elemts to print for vector parameters (default: 10; 0 -> all)");
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

GGSDetectorConstruction::GGSDetectorConstruction(const G4String &library, const G4String &,
                                                 const G4String &configDataCard, bool)
    : _physicalWorld(nullptr), _geometry(nullptr), _library(library), _configDataCard(configDataCard), _gdxml(""),
      _messenger(this, "/GGS/"), _nToPrintForVectGeoParams{10} {
  _messenger.DeclareProperty("nToPrintForVectGeoParams", _nToPrintForVectGeoParams,
                             "Maximum number of elemts to print for vector parameters (default: 10; 0 -> all)");
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

GGSDetectorConstruction::GGSDetectorConstruction(const G4String &gdxml)
    : _physicalWorld(nullptr), _geometry(nullptr), _library(""), _configDataCard(""), _gdxml(gdxml),
      _messenger(this, "/GGS/"), _nToPrintForVectGeoParams{10} {
  _messenger.DeclareProperty("nToPrintForVectGeoParams", _nToPrintForVectGeoParams,
                             "Maximum number of elemts to print for vector parameters (default: 10; 0 -> all)");
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

GGSDetectorConstruction::~GGSDetectorConstruction() { delete _geometry; }

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

namespace {
#ifdef USE_GDML
// Geometry construction class for top-level GDML file.
// This class provides a thin and feature-less geometry construction to be used when
// a GDML file is used instead of a geometry library.
class GGSGDMLGeometryConstruction : public GGSVGeometryConstruction {
public:
  // The world parameter must be a pointer to the already constructed world volume.
  GGSGDMLGeometryConstruction(G4VPhysicalVolume *world) : _world(world) {}
  G4VPhysicalVolume *Construct() { return _world; }
  G4VPhysicalVolume *GetVolume() { return _world; };

private:
  G4VPhysicalVolume *_world;
};
#endif

#ifdef GGS_USE_DD4HEP
// Geometry construction class for DD4hep geometries.
// This class provides a thin and feature-less geometry construction to be used when
// a DD4hep geometry is used instead of a geometry library.
class GGSDD4hepGeometryConstruction : public GGSVGeometryConstruction {
public:
  // The world parameter must be a pointer to the already constructed world volume.
  GGSDD4hepGeometryConstruction(G4VPhysicalVolume *world) : _world(world) {}
  G4VPhysicalVolume *Construct() { return _world; }
  G4VPhysicalVolume *GetVolume() { return _world; };

private:
  G4VPhysicalVolume *_world;
};
#endif

template <typename T> void PrintVectParams(const T &paramsVect, size_t nToPrint) {
  for (auto &par : paramsVect) {
    GGSCCOUT(INFO) << par.first << ": {";
    for (size_t iVal = 0; iVal < std::min(par.second.size() - 1, nToPrint - 1); ++iVal) {
      GGSCOUT_BACKEND << par.second[iVal] << ", ";
    }
    if (nToPrint != 0 && par.second.size() > nToPrint) {
      GGSCOUT_BACKEND << "..., ";
    }
    GGSCOUT_BACKEND << par.second.back() << "}\n";
  }
}

} // namespace

G4VPhysicalVolume *GGSDetectorConstruction::Construct() {

  static const std::string routineName("GGSDetectorConstruction::Construct");

  // 1. Construct the geometry
  if (_library != "") {
    // 1.1 Load the concrete GGSVGeometryConstruction from a shared library
    if (!(GGSGeoPluginManager::GetInstance().LoadGeoPlugin(_library))) {
      GGSCOUT(ERROR) << "Impossible to load geometry library " << _library << GGSENDL;
      throw std::runtime_error("Impossible to load geometry library");
    }
    _geometry = GGSGeoPluginManager::GetInstance().GetGeoConstruction();
    if (!_geometry) {
      GGSCOUT(ERROR) << "Can't create the geometry." << GGSENDL;
    }
    GGSCOUT(INFO) << "Construct the detector." << GGSENDL;
    _geometry->SetGeoDataCard(_configDataCard);
    _physicalWorld = _geometry->Construct();
    // Print geometry version
    const std::string geoVersion = _geometry->GetVersion();
    if (geoVersion != "") {
      GGSCOUT(INFO) << "Geometry version: " << geoVersion << GGSENDL;
    } else {
      GGSCOUT(INFO) << "No geometry version available" << GGSENDL;
    }
    // Print geometry parameters
    _geometry->ExportParameters(); // Will throw if something goes wrong, so don't check the return value
    auto intParams = _geometry->GetIntParameters();
    auto boolParams = _geometry->GetBoolParameters();
    auto realParams = _geometry->GetRealParameters();
    auto stringParams = _geometry->GetStringParameters();
    auto vectIntParams = _geometry->GetVectIntParameters();
    auto vectBoolParams = _geometry->GetVectBoolParameters();
    auto vectRealParams = _geometry->GetVectRealParameters();
    auto vectStringParams = _geometry->GetVectStringParameters();
    if (intParams.size() != 0 || boolParams.size() != 0 || realParams.size() != 0 || stringParams.size() != 0 ||
        vectIntParams.size() != 0 || vectBoolParams.size() != 0 || vectRealParams.size() != 0 ||
        vectStringParams.size() != 0) {
      GGSCOUT(INFO) << "Geometry parameters:\n";
      for (auto &par : intParams) {
        GGSCCOUT(INFO) << par.first << ": " << par.second << "\n";
      }
      for (auto &par : boolParams) {
        GGSCCOUT(INFO) << par.first << ": " << (par.second ? "true" : "false") << "\n";
      }
      for (auto &par : realParams) {
        GGSCCOUT(INFO) << par.first << ": " << par.second << "\n";
      }
      for (auto &par : stringParams) {
        GGSCCOUT(INFO) << par.first << ": " << par.second << "\n";
      }
      PrintVectParams(vectIntParams, _nToPrintForVectGeoParams);
      PrintVectParams(vectRealParams, _nToPrintForVectGeoParams);
      PrintVectParams(vectBoolParams, _nToPrintForVectGeoParams);
      PrintVectParams(vectStringParams, _nToPrintForVectGeoParams);
    }
  }
#ifdef USE_GDML
  else if (_gdxml.length() > 4 && _gdxml.substr(_gdxml.length() - 5) == ".gdml") {
    // 1.2 Build the geometry from a GDML file
    G4GDMLParser gdmlParser;
    gdmlParser.SetOverlapCheck(false);
    try {
      gdmlParser.Read(_gdxml, false);
    } catch (int &ge) {
      GGSCOUT(ERROR) << "G4GDML: " << ge << ", missing network connection? wrong schema URL?" << GGSENDL;
      // G4GGSCOUT << "Try to read GDML file _WITHOUT_ schema validation" << G4GGSENDL;
      // error, try without schema validation
      // gdmlParser.ValidateSchema(false);
      // gdmlParser.Read(_gdxml, false);
    }
    gdmlParser.StripNamePointers();
    _physicalWorld = gdmlParser.GetWorldVolume();
    if (!_physicalWorld) {
      GGSCOUT(ERROR) << "Cannot build the GMDL geometry." << GGSENDL;
      throw std::runtime_error("Cannot build the GMDL geometry.");
    }
    _geometry = new GGSGDMLGeometryConstruction(_physicalWorld);
  }
#endif
#ifdef GGS_USE_DD4HEP
  else if (_gdxml.length() > 3 && _gdxml.substr(_gdxml.length() - 4) == ".xml") {
    // 1.3 Build the DD4hep geometry
    dd4hep::Detector &detDesc = dd4hep::Detector::getInstance();
    detDesc.fromXML(_gdxml);
    dd4hep::DetElement world = detDesc.world();
    dd4hep::sim::Geant4Converter conv(detDesc, dd4hep::NOLOG);
    dd4hep::sim::Geant4GeometryInfo *geoInfo = conv.create(world).detach();
    _physicalWorld = geoInfo->world();
    _geometry = new GGSDD4hepGeometryConstruction(_physicalWorld);
  }
#endif

  // 2. Construct the standard SDs
  // This feature is removed from GGS, so print an error message and throw an exception.

  const G4LogicalVolumeStore *logicalVolumeStore = G4LogicalVolumeStore::GetInstance();
  std::vector<G4LogicalVolume *>::const_iterator i;
  GGSNameDecoder &nameDecoder = GGSNameDecoder::GetInstance();
  for (i = logicalVolumeStore->begin(); i != logicalVolumeStore->end(); i++) {
    std::string logVolName = (*i)->GetName();
    auto detVolName = logVolName.substr(0, logVolName.find_first_of(' '));
    auto spacePos = logVolName.find_last_of(' ');
    std::string scorerName("");
    if (spacePos != std::string::npos)
      scorerName = logVolName.substr(spacePos + 1);
    // check if it is a sensitive volume
    if (nameDecoder.IsSensitive(detVolName)) {
      GGSCOUT(WARNING) << "This version of GGS does not support the definition of sensitive volumes by name ("
                       << detVolName << ").\n";
      GGSCCOUT(WARNING) << "Please change the name of the volume and make it active using a datacard command."
                        << GGSENDL;
      throw std::runtime_error("Use of removed feature: definition of sensitive volume by name.");
    }
  }

  return _physicalWorld;
}

void GGSDetectorConstruction::ConstructSDandField() {
  // Construct the build commands for the sensitive detectors
  sdBuildMessenger = std::make_unique<SDBuildMessenger>("/GGS/scoring/", "Messenger for SD build commands");
  auto &availableSDBuilders =
      GGSFactory<G4VSensitiveDetector, std::string>::GetInstance().GetListOfRegisteredBuilders();
  for (auto &builder : availableSDBuilders) {
    // Create the command
    sdBuildMessenger->AddSDBuildCommand("/GGS/scoring/add" + builder);
  }
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......