uLib/src/HEP/Geant/Scene.cpp

/*//////////////////////////////////////////////////////////////////////////////
// CMT Cosmic Muon Tomography project //////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////

  Copyright (c) 2014, Universita' degli Studi di Padova, INFN sez. di Padova
  All rights reserved

  Authors: Andrea Rigoni Garola < andrea.rigoni@pd.infn.it >

  ------------------------------------------------------------------
  This library is free software;  you  can  redistribute  it  and/or
  modify it  under the  terms  of  the  GNU  Lesser  General  Public
  License as published  by  the  Free  Software  Foundation;  either
  version 3.0 of the License, or (at your option) any later version.

  This library is  distributed in  the hope that it will  be useful,
  but  WITHOUT ANY WARRANTY;  without  even  the implied warranty of
  MERCHANTABILITY  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.

  You should have received a copy of  the GNU Lesser General  Public
  License along with this library.

//////////////////////////////////////////////////////////////////////////////*/

#include <Geant4/G4Box.hh>
#include <Geant4/G4LogicalVolume.hh>
#include <Geant4/G4Material.hh>
#include <Geant4/G4NistManager.hh>
#include <Geant4/G4PVPlacement.hh>
#include <Geant4/G4RunManager.hh>
#include <Geant4/G4RunManagerFactory.hh>
#include <Geant4/G4SystemOfUnits.hh>
#include <Geant4/G4VPhysicalVolume.hh>

#include "Core/Vector.h"
#include "HEP/Geant/DetectorConstruction.hh"
#include "Math/ContainerBox.h"
#include "Math/Dense.h"
#include "Solid.h"

#include "Scene.h"
#include "PhysicsList.hh"
#include "ActionInitialization.hh"

namespace uLib {
namespace Geant {

class SceneDetectorConstruction : public DetectorConstruction {
public:
  SceneDetectorConstruction(class SceneImpl *owner);
  G4VPhysicalVolume *Construct() override;
private:
  class SceneImpl *m_Owner;
};

static void CheckGeant4Environment() {
  static bool checked = false;
  if (checked) return;
  checked = true;
  if (!std::getenv("G4ENSDFSTATEDATA")) {
    std::cerr << "********************************************************" << std::endl;
    std::cerr << " WARNING: Geant4 environment variables are not set!" << std::endl;
    std::cerr << " Please activate the environment before running:" << std::endl;
    std::cerr << "   micromamba activate mutom" << std::endl;
    std::cerr << "********************************************************" << std::endl;
  }
}

class SceneImpl {
public:
  // constructor //
  SceneImpl() : m_RunManager(G4RunManagerFactory::CreateRunManager(G4RunManagerType::Serial)),
                m_Emitter(nullptr),
                m_Output(nullptr) {
    m_RunManager->SetUserInitialization(new PhysicsList);
  }

  // destructor //
  ~SceneImpl() { 
    if (m_RunManager) delete m_RunManager; 
    if (m_World) delete m_World; 
  }

  void Initialize() {
    // Set mandatory initialization classes for Geant4
    m_RunManager->SetUserInitialization(new SceneDetectorConstruction(this));
    m_RunManager->SetUserInitialization(
        new ActionInitialization(m_Emitter, m_Output));

    // Initialize Geant4
    m_RunManager->Initialize();
  }

  // members //
  Vector<Solid *> m_Solids;
  Solid *m_World = nullptr;
  ContainerBox m_WorldBox;
  G4RunManager *m_RunManager;
  EmitterPrimary *m_Emitter;
  Vector<GeantEvent> *m_Output;
};

SceneDetectorConstruction::SceneDetectorConstruction(SceneImpl *owner) 
  : DetectorConstruction("Scene"), m_Owner(owner) {}

G4VPhysicalVolume *SceneDetectorConstruction::Construct() {
  return m_Owner->m_World->GetPhysical();
}




Scene::Scene() {
  CheckGeant4Environment();
  d = new SceneImpl();
}
Scene::~Scene() { delete d; }

void Scene::AddSolid(Solid *solid, Solid *parent) { 
  d->m_Solids.push_back(solid);
  if (!d->m_World) {
    d->m_World = solid;
  } else {
    solid->SetParent(parent ? parent : d->m_World);
  }
}

const Solid* Scene::GetWorld() const { return d->m_World; }

ContainerBox* Scene::GetWorldBox() const { return &d->m_WorldBox; }

void Scene::ConstructWorldBox(const Vector3f &size, const char *material) {
  // Get nist material manager

  d->m_WorldBox.Scale(size);
  d->m_WorldBox.SetPosition(-size/2.0f);

  if (!d->m_World) {
    d->m_World = new Solid("World");
    d->m_World->SetNistMaterial(material);
    AddSolid(d->m_World);
  }

  G4Box *solidWorld = new G4Box("World", 
                                0.5 * size(0), 
                                0.5 * size(1),
                                0.5 * size(2));

  G4LogicalVolume *logicWorld = new G4LogicalVolume(solidWorld,
                                                    d->m_World->GetMaterial(),
                                                    d->m_World->GetName());

  d->m_World->SetLogical(logicWorld);
  
  G4PVPlacement *physWorld = new G4PVPlacement(
    nullptr,
    G4ThreeVector(0, 0, 0),
    logicWorld,
    d->m_World->GetName(),
    0,
    false,
    0,
    true);

  d->m_World->SetPhysical(physWorld);

  // no transforms are allowed for the world box
  // Matrix4f transform = box->GetMatrix(); 
  // d->m_World->SetTransform(transform);
}

void Scene::SetEmitter(EmitterPrimary *emitter) {
  d->m_Emitter = emitter;
}

void Scene::Initialize() {
  d->Initialize();
}

void Scene::RunSimulation(int nEvents, Vector<GeantEvent> &results) {
  d->m_Output = &results;

  // Re-initialize ActionInitialization with the output buffer
  // (ActionInitialization was already set during Initialize, but we need
  //  to ensure the output pointer is current)
  d->m_RunManager->SetUserInitialization(
      new ActionInitialization(d->m_Emitter, &results));

  d->m_RunManager->BeamOn(nEvents);
}

} // namespace Geant
} // namespace uLib