97 lines
3.2 KiB
C++
97 lines
3.2 KiB
C++
#include "Geant/Solid.h"
|
|
#include "HEP/Geant/GeantEvent.h"
|
|
#include "HEP/Geant/Scene.h"
|
|
#include "HEP/Geant/EmitterPrimary.hh"
|
|
#include "Math/ContainerBox.h"
|
|
#include "Math/TriangleMesh.h"
|
|
#include "Math/Dense.h"
|
|
#include "Math/Units.h"
|
|
#include "testing-prototype.h"
|
|
#include <Geant4/G4Material.hh>
|
|
#include <Geant4/G4NistManager.hh>
|
|
#include <Geant4/G4LogicalVolume.hh>
|
|
#include <Geant4/G4TessellatedSolid.hh>
|
|
#include <string.h>
|
|
|
|
using namespace uLib;
|
|
|
|
int main() {
|
|
BEGIN_TESTING(Geant Event);
|
|
|
|
// Test: Scene with iron cube in air, launch muons, collect events //
|
|
{
|
|
// 1. Create world box (air, 30m x 30m x 30m)
|
|
Geant::Scene scene;
|
|
scene.ConstructWorldBox(Vector3f(30_m, 30_m, 30_m), "G4_AIR");
|
|
|
|
// 2. Create iron cube (1m x 1m x 1m) at center
|
|
ContainerBox *iron_box = new ContainerBox(Vector3f(1000, 1000, 1000)); // mm
|
|
Geant::BoxSolid *iron_cube = new Geant::BoxSolid("IronCube", iron_box);
|
|
|
|
Geant::Material *iron_mat = new Geant::Material("G4_Fe");
|
|
Geant::LogicalVolume *iron_lv = new Geant::LogicalVolume("IronCube_lv");
|
|
iron_lv->SetSolid(iron_cube);
|
|
iron_lv->SetMaterial(iron_mat);
|
|
iron_lv->Update();
|
|
|
|
Geant::PhysicalVolume *iron_pv = new Geant::PhysicalVolume("IronCube", iron_lv);
|
|
scene.AddVolume(iron_pv);
|
|
|
|
// 3. Set up emitter (default: mu- at 1 GeV, from z=+10m downward)
|
|
Geant::EmitterPrimary *emitter = new Geant::EmitterPrimary();
|
|
scene.SetEmitter(emitter);
|
|
|
|
// 4. Initialize Geant4
|
|
scene.Initialize();
|
|
|
|
// 5. Run simulation: 10 muons
|
|
int nEvents = 10;
|
|
Vector<Geant::GeantEvent> results;
|
|
scene.RunSimulation(nEvents, results);
|
|
|
|
// 6. Check results
|
|
printf(" Collected %zu events\n", results.size());
|
|
TEST1(results.size() > 0);
|
|
|
|
for (size_t i = 0; i < results.size(); ++i) {
|
|
const Geant::GeantEvent &ev = results[i];
|
|
bool hitIron = false;
|
|
for (const auto &d : ev.Path()) {
|
|
if (d.SolidName() == "IronCube") {
|
|
hitIron = true;
|
|
break;
|
|
}
|
|
}
|
|
|
|
printf(" Event %d: momentum=%.1f MeV, path steps=%zu, hitIron=%s\n",
|
|
ev.GetEventID(),
|
|
ev.GetMomentum(),
|
|
ev.Path().size(),
|
|
hitIron ? "YES" : "NO");
|
|
|
|
// Each event should have at least one step
|
|
TEST1(ev.Path().size() > 0);
|
|
|
|
// Print first few deltas
|
|
const auto &path = ev.Path();
|
|
for (size_t j = 0; j < path.size() && j < 5; ++j) {
|
|
const auto &d = path[j];
|
|
printf(" Delta[%zu]: solid=%s len=%.2f mm, p=%.1f MeV, "
|
|
"dir=(%.3f, %.3f, %.3f)\n",
|
|
j,
|
|
d.SolidName().c_str(),
|
|
d.GetLength(),
|
|
d.GetMomentum(),
|
|
d.Direction()(0),
|
|
d.Direction()(1),
|
|
d.Direction()(2));
|
|
}
|
|
if (path.size() > 5) {
|
|
printf(" ... (%zu more deltas)\n", path.size() - 5);
|
|
}
|
|
}
|
|
}
|
|
|
|
END_TESTING
|
|
}
|