lars/minicalosim
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'main' of gitlab.etp.kit.edu:jkiesele/minicalosim

Browse Source
This commit is contained in:
lars
2023-10-07 15:31:02 +02:00
parent a70bb85247 6e591c6255
commit e7b26da63f
15 changed files with 110 additions and 1059 deletions
Download Patch File Download Diff File Expand all files Collapse all files
bind/G4Calo.py
+8 -13
View File
@@ -1,4 +1,4 @@
from minicalo import ConstructionWrapper
from minicalo import GeometryDescriptor
from minicalo import G4System as _G4System
import plotly.graph_objects as go
@@ -10,15 +10,8 @@ import uproot
from IPython.display import Image, display
class G4System(_G4System):
def run_visualize(
self, particleSpec: str, minEnergy_GeV: float, maxEnergy_GeV: float = -1.0
):
if maxEnergy_GeV < 0:
maxEnergy_GeV = minEnergy_GeV
# anpassen: run 1 event and get stuff from construction wrapper
_G4System.run_visualize(self, particleSpec, minEnergy_GeV, maxEnergy_GeV)
self.displayEvent()
class __G4System(_G4System):
def run_batch(
self,
@@ -36,7 +29,7 @@ class G4System(_G4System):
df = ttree["Hits;1"].arrays(library="pd")
return df
def displayEvent(self, gd):
def displayEvent(self):
# for loop over all layers
to_plot = []
material_dict = {}
@@ -44,7 +37,7 @@ class G4System(_G4System):
# sum up total deposited energy
total_dep_energy = 0
for layer in gd.getLayers():
for layer in self.getGeometryDescriptor().getLayers():
for sensor in layer.sensors:
total_dep_energy += sensor.getEnergy()
@@ -53,7 +46,7 @@ class G4System(_G4System):
total_dep_energy = 10**-8 # to avoid division by zero
# loop over materials
for layer in gd.getLayers():
for layer in self.getGeometryDescriptor().getLayers():
layer_width = layer.nx * layer.sens_xwidth
@@ -246,3 +239,5 @@ col_dict = {
"G4_Cu": 'dimgray',
"G4_BRASS": 'slategrey',
}
G4System = __G4System()#singleton instance
bind/bindings.cpp
+13 -11
View File
@@ -3,22 +3,23 @@
#include <pybind11/stl.h>
#include <pybind11/operators.h>
#include "ConstructionWrapper.hh"
#include "GeometryDescriptor.hh"
#include "G4System.hh"
namespace py = pybind11;
template<class M>
void makeConstructionWrapper(M & m, std::string name){
void makeGeometryDescriptor(M & m, std::string name){
py::class_<ConstructionWrapper>(m, name.data()).def(py::init())
.def("addLayer", &ConstructionWrapper::addLayer, py::arg("thickness"), py::arg("material"), py::arg("isActive")=true, py::arg("nx")=1, py::arg("ny")=-1)
.def("getXYWidth", &ConstructionWrapper::getXYWidth)
py::class_<GeometryDescriptor>(m, name.data()).def(py::init())
.def("addLayer", &GeometryDescriptor::addLayer, py::arg("thickness"), py::arg("material"), py::arg("isActive")=true, py::arg("nx")=1, py::arg("ny")=-1)
.def("getXYWidth", &GeometryDescriptor::getXYWidth)
// bind overloaded getLayers function
.def("getLayers", (std::vector<Layer> & (ConstructionWrapper::*)()) &ConstructionWrapper::getLayers)
.def("getLayers", (const std::vector<Layer> & (ConstructionWrapper::*)() const) &ConstructionWrapper::getLayers)
.def("getNSensors", &ConstructionWrapper::getNSensors);
.def("getLayers", (std::vector<Layer> & (GeometryDescriptor::*)()) &GeometryDescriptor::getLayers)
.def("getLayers", (const std::vector<Layer> & (GeometryDescriptor::*)() const) &GeometryDescriptor::getLayers)
.def("isAssigned", &GeometryDescriptor::isAssigned)
.def("getNSensors", &GeometryDescriptor::getNSensors);
}
@@ -31,8 +32,9 @@ void makeG4System(M &m, std::string name)
.def("run_gui", &G4System::run_gui).def("run_batch", &G4System::run_batch, py::arg("nEvents"), py::arg("partSpecies"), py::arg("minEnergy_GeV"), py::arg("maxEnergy_GeV"))
.def("applyUICommand", &G4System::applyUICommand, py::arg("command"))
.def("displayEvent", &G4System::displayEvent)
.def("printMaterial", &G4System::printMaterial, py::arg("name"))
.def("check", &G4System::check);
// .def("printMaterial", &G4System::printMaterial, py::arg("name"))
//.def("check", &G4System::check)
.def("getGeometryDescriptor", &G4System::getGeometryDescriptor);
}
// create bindings for Layer class
@@ -69,7 +71,7 @@ void makeSensor(M &m, std::string name){
PYBIND11_MODULE(minicalo, m)
{
m.doc() = "pybind11 plugin"; // optional module docstring
makeConstructionWrapper(m, "ConstructionWrapper");
makeGeometryDescriptor(m, "GeometryDescriptor");
makeG4System(m, "G4System");
makeSensor(m, "Sensor");
makeLayer(m, "Layer");
bind/example.ipynb
+7 -978
View File
File diff suppressed because one or more lines are too long
bind/example.py
+6 -6
View File
@@ -1,14 +1,14 @@
from G4Calo import ConstructionWrapper, G4System
from G4Calo import GeometryDescriptor, G4System
cw = ConstructionWrapper()#the width of the calorimeter is 50 cm times 50 cm (also steerable, but I'd leave it)
cw = GeometryDescriptor()#the width of the calorimeter is 50 cm times 50 cm (also steerable, but I'd leave it)
cw.addLayer(10,"G4_Si",True) #10 cm of Silicon, active
cw.addLayer(1,"G4_Pb",False) #1 cm of Lead, passive
cw.addLayer(10,"G4_Si",True, 3) #10 cm of Silicon, active, 3x3 segmentation
G4s = G4System()
G4s.init(cw)
G4s.run_batch(100, "gamma", 1, 5) #100 events, photons, 1000 MeV to 5000 MeV
G4System = G4System
G4System.init(cw)
G4System.run_batch(100, "gamma", 1, 5) #100 events, photons, 1000 MeV to 5000 MeV
#G4s.run_visualize("gamma", 1, 5)
G4s.printMaterial("G4_Pb")
G4System.printMaterial("G4_Pb")
exampleB4a.cc
+6 -6
View File
@@ -39,7 +39,7 @@
#include "FTFP_BERT.hh"
#include "Randomize.hh"
#include "ConstructionWrapper.hh"
#include "GeometryDescriptor.hh"
#include "G4System.hh"
@@ -121,7 +121,7 @@ int main2(int argc,char** argv)
}
#endif
ConstructionWrapper * cw = new ConstructionWrapper();
GeometryDescriptor * cw = new GeometryDescriptor();
cw->addLayer(1, "G4_Pb", false);
cw->addLayer(10, "G4_Si", true, 9);
cw->addLayer(1, "G4_Pb", false);
@@ -185,11 +185,11 @@ int main2(int argc,char** argv)
//later
//class Runner {
//public:
// void initialize(ConstructionWrapper CW, bool gui=false);
// void initialize(GeometryDescriptor CW, bool gui=false);
// void run(int nEvents, std::string partSpecies, double minEnergy, double maxEnergy);
//}
void run(ConstructionWrapper CW, int nEvents, std::string partSpecies, double minEnergy, double maxEnergy, bool gui=false){
void run(GeometryDescriptor CW, int nEvents, std::string partSpecies, double minEnergy, double maxEnergy, bool gui=false){
char* argv[]={(char*)"dummy"};
@@ -199,7 +199,7 @@ void run(ConstructionWrapper CW, int nEvents, std::string partSpecies, double mi
G4RunManagerFactory::CreateRunManager(G4RunManagerType::Default);
runManager->SetNumberOfThreads(1);
ConstructionWrapper * cw = &CW;
GeometryDescriptor * cw = &CW;
G4String session;
G4UIExecutive* ui = nullptr;
if ( gui ) {
@@ -261,7 +261,7 @@ void run(ConstructionWrapper CW, int nEvents, std::string partSpecies, double mi
int main(){
ConstructionWrapper cw;
GeometryDescriptor cw;
//to be moved
cw.addLayer(1, "G4_Pb", false);
cw.addLayer(10, "G4_Si", true, 9);
include/DetectorConstruction.hh
+6 -6
View File
@@ -35,7 +35,7 @@
class G4VPhysicalVolume;
class G4GlobalMagFieldMessenger;
class ConstructionWrapper;
class GeometryDescriptor;
namespace B4
{
@@ -67,7 +67,7 @@ class DetectorConstruction : public G4VUserDetectorConstruction
return _global_detector_construction;
}
DetectorConstruction( ConstructionWrapper* Cw): G4VUserDetectorConstruction(), cw(Cw){
DetectorConstruction( GeometryDescriptor* Cw): G4VUserDetectorConstruction(), cw(Cw){
_global_detector_construction=this;
}
~DetectorConstruction() override = default;
@@ -76,14 +76,14 @@ class DetectorConstruction : public G4VUserDetectorConstruction
G4VPhysicalVolume* Construct() override;
void ConstructSDandField() override;
void setConstructionWrapper(ConstructionWrapper* Cw){
void setGeometryDescriptor(GeometryDescriptor* Cw){
this->cw = Cw;
}
ConstructionWrapper* getConstructionWrapper(){
GeometryDescriptor* getGeometryDescriptor(){
return cw;
}
const ConstructionWrapper* getConstructionWrapper()const{
const GeometryDescriptor* getGeometryDescriptor()const{
return cw;
}
@@ -100,7 +100,7 @@ class DetectorConstruction : public G4VUserDetectorConstruction
G4bool fCheckOverlaps = true; // option to activate checking of volumes overlaps
ConstructionWrapper* cw=nullptr;
GeometryDescriptor* cw=nullptr;
};
include/EventAction.hh
+1 -1
View File
@@ -34,7 +34,7 @@
#include "RunAction.hh"
#include "PrimaryGeneratorAction.hh"
#include "globals.hh"
class ConstructionWrapper;
class GeometryDescriptor;
namespace B4a
{
include/G4System.hh
+14 -4
View File
@@ -2,7 +2,7 @@
#ifndef G4SYSTEM_HH
#define G4SYSTEM_HH
#include "ConstructionWrapper.hh"
#include "GeometryDescriptor.hh"
#include <stdexcept>
#include "G4RunManagerFactory.hh"
#include "DetectorConstruction.hh"
@@ -28,10 +28,11 @@ G4System(bool Gui=false):gui(Gui){};
delete runManager;
delete ui;
}
};
void init(ConstructionWrapper &cw);
void init(GeometryDescriptor &cw);
//will use dawn for visualization, also wrap more in python
void run_visualize(const std::string& partSpecies, double minEnergy_GeV, double maxEnergy_GeV);
//runs the whole gui if available
@@ -47,9 +48,18 @@ void displayEvent()const{}; //just a placeholder, this will be implemented in py
void printMaterial(const std::string& name)const;
GeometryDescriptor& getGeometryDescriptor(){
check();
if ( assigned_cw == nullptr ){
throw std::runtime_error("GeometryDescriptor not assigned");
}
return *assigned_cw;
}
private:
void check()const;
GeometryDescriptor * assigned_cw=nullptr;
bool gui;
@@ -67,7 +77,7 @@ G4UIExecutive * ui=nullptr;
/*
void run(ConstructionWrapper CW, int nEvents, std::string partSpecies, double minEnergy, double maxEnergy, bool gui=false){
void run(GeometryDescriptor CW, int nEvents, std::string partSpecies, double minEnergy, double maxEnergy, bool gui=false){
char* argv[]={(char*)"dummy"};
@@ -77,7 +87,7 @@ void run(ConstructionWrapper CW, int nEvents, std::string partSpecies, double mi
G4RunManagerFactory::CreateRunManager(G4RunManagerType::Default);
runManager->SetNumberOfThreads(1);
ConstructionWrapper * cw = &CW;
GeometryDescriptor * cw = &CW;
G4String session;
G4UIExecutive* ui = nullptr;
if ( gui ) {
include/ConstructionWrapper.hh → include/GeometryDescriptor.hh
+16 -5
View File
@@ -1,6 +1,6 @@
#ifndef CONSTRUCTIONWRAPPER_HH
#define CONSTRUCTIONWRAPPER_HH
#ifndef GeometryDescriptor_HH
#define GeometryDescriptor_HH
#include "G4ThreeVector.hh"
#include <vector>
@@ -66,6 +66,11 @@ public:
void setIsActive(bool isActive);
void assignPhysicalVolume(G4VPhysicalVolume* physicalVolume);
void unAssign(){
physicalVolume = nullptr;
sensors.clear();
}
double thickness;
double sens_xwidth;
double sens_ywidth;
@@ -81,11 +86,11 @@ public:
};
class ConstructionWrapper{
class GeometryDescriptor{
public:
ConstructionWrapper(double xy_width=50);
~ConstructionWrapper() {};
GeometryDescriptor(double xy_width=50);
~GeometryDescriptor() {};
void addLayer(double thickness_cm, std::string material, bool isActive=true, int nx=1, int ny=-1);
@@ -123,6 +128,12 @@ public:
}
}
void unAssign(){
for(auto& layer: layers){
layer.unAssign();
}
}
private:
double xywidth;
include/SteppingAction.hh
+1 -1
View File
@@ -31,7 +31,7 @@
#define B4aSteppingAction_h 1
#include "G4UserSteppingAction.hh"
#include "ConstructionWrapper.hh"
#include "GeometryDescriptor.hh"
#include "DetectorConstruction.hh"
src/DetectorConstruction.cc
+1 -1
View File
@@ -28,7 +28,7 @@
/// \brief Implementation of the B4::DetectorConstruction class
#include "DetectorConstruction.hh"
#include "ConstructionWrapper.hh"
#include "GeometryDescriptor.hh"
#include "G4Material.hh"
#include "G4NistManager.hh"
src/EventAction.cc
+5 -5
View File
@@ -27,7 +27,7 @@
/// \file B4/B4a/src/EventAction.cc
/// \brief Implementation of the B4a::EventAction class
#include "ConstructionWrapper.hh"
#include "GeometryDescriptor.hh"
#include "EventAction.hh"
#include "RunAction.hh"
@@ -49,11 +49,11 @@ namespace B4a
void EventAction::BeginOfEventAction(const G4Event* /*event*/)
{
auto cw = B4::DetectorConstruction::getDetectorConstruction()->getConstructionWrapper();
auto cw = B4::DetectorConstruction::getDetectorConstruction()->getGeometryDescriptor();
// initialisation per event
if(cw == nullptr){
G4cout << "ConstructionWrapper not found" << G4endl;
throw std::runtime_error("ConstructionWrapper not found");
G4cout << "GeometryDescriptor not found" << G4endl;
throw std::runtime_error("GeometryDescriptor not found");
}
if(gen == nullptr){
G4cout << "PrimaryGeneratorAction not found" << G4endl;
@@ -84,7 +84,7 @@ void EventAction::EndOfEventAction(const G4Event* event)
// get analysis manager
auto analysisManager = G4AnalysisManager::Instance();
auto cw = B4::DetectorConstruction::getDetectorConstruction()->getConstructionWrapper();
auto cw = B4::DetectorConstruction::getDetectorConstruction()->getGeometryDescriptor();
//cw->printSensorEnergies();//DEBUG
src/G4System.cc
+8 -4
View File
@@ -4,11 +4,16 @@
#include <G4NistManager.hh>
#include <G4String.hh>
void G4System::init(ConstructionWrapper &CW){
void G4System::init(GeometryDescriptor &CW){
if(CW.isAssigned()){
throw std::runtime_error("ConstructionWrapper already assigned");
throw std::runtime_error("GeometryDescriptor already assigned");
}
if(assigned_cw != nullptr){
//reassign
assigned_cw->unAssign();
}
assigned_cw = &CW;
char* argv[]={(char*)"dummy"};
@@ -33,7 +38,6 @@ void G4System::init(ConstructionWrapper &CW){
//delete actionInitialization;
}
G4String session;
ui = nullptr;
if ( gui ) {
@@ -54,7 +58,7 @@ void G4System::init(ConstructionWrapper &CW){
runManager->SetUserInitialization(actionInitialization);
}
else{
detConstruction->setConstructionWrapper(&CW);
detConstruction->setGeometryDescriptor(&CW);
runManager->ReinitializeGeometry(true);
}
src/ConstructionWrapper.cc → src/GeometryDescriptor.cc
+12 -12
View File
@@ -1,6 +1,6 @@
#include "ConstructionWrapper.hh"
#include "GeometryDescriptor.hh"
#include "G4VPhysicalVolume.hh"
Layer::Layer(){
@@ -36,11 +36,11 @@ void Layer::assignPhysicalVolume(G4VPhysicalVolume* physicalVolume){
this->physicalVolume = physicalVolume;
}
ConstructionWrapper::ConstructionWrapper(double xy_width){
GeometryDescriptor::GeometryDescriptor(double xy_width){
xywidth = xy_width;
}
void ConstructionWrapper::addLayer(double thickness, std::string material, bool isActive, int nx, int ny){
void GeometryDescriptor::addLayer(double thickness, std::string material, bool isActive, int nx, int ny){
if(ny<0){
ny = nx;
}
@@ -60,20 +60,20 @@ void ConstructionWrapper::addLayer(double thickness, std::string material, bool
layers.push_back(layer);
}
const std::vector<Layer> & ConstructionWrapper::getLayers() const{
const std::vector<Layer> & GeometryDescriptor::getLayers() const{
return layers;
}
std::vector<Layer> & ConstructionWrapper::getLayers(){
std::vector<Layer> & GeometryDescriptor::getLayers(){
return layers;
}
double ConstructionWrapper::getXYWidth() const{
double GeometryDescriptor::getXYWidth() const{
return xywidth;
}
void ConstructionWrapper::resetSensorEnergies()const {
void GeometryDescriptor::resetSensorEnergies()const {
for(const auto & layer : layers){
for(const auto & sensor : layer.sensors){
sensor.energy = 0;
@@ -81,7 +81,7 @@ void ConstructionWrapper::resetSensorEnergies()const {
}
}
Layer* ConstructionWrapper::getLayerByVolume(G4VPhysicalVolume* volume){
Layer* GeometryDescriptor::getLayerByVolume(G4VPhysicalVolume* volume){
for(auto & layer : layers){
if(layer.physicalVolume == volume){
return &layer;
@@ -90,7 +90,7 @@ Layer* ConstructionWrapper::getLayerByVolume(G4VPhysicalVolume* volume){
return nullptr;
}
const Layer* ConstructionWrapper::getLayerByVolume(G4VPhysicalVolume* volume)const{
const Layer* GeometryDescriptor::getLayerByVolume(G4VPhysicalVolume* volume)const{
for(auto & layer : layers){
if(layer.physicalVolume == volume){
return &layer;
@@ -99,7 +99,7 @@ const Layer* ConstructionWrapper::getLayerByVolume(G4VPhysicalVolume* volume)con
return nullptr;
}
Sensor* ConstructionWrapper::getSensorByVolume(G4VPhysicalVolume* volume){
Sensor* GeometryDescriptor::getSensorByVolume(G4VPhysicalVolume* volume){
int copyNo = volume->GetCopyNo();
auto layer = getLayerByVolume(volume);
if(layer == nullptr){
@@ -111,7 +111,7 @@ Sensor* ConstructionWrapper::getSensorByVolume(G4VPhysicalVolume* volume){
return &layer->sensors[copyNo];
}
const Sensor* ConstructionWrapper::getSensorByVolume(G4VPhysicalVolume* volume)const{
const Sensor* GeometryDescriptor::getSensorByVolume(G4VPhysicalVolume* volume)const{
int copyNo = volume->GetCopyNo();
auto layer = getLayerByVolume(volume);
if(layer == nullptr){
@@ -124,7 +124,7 @@ const Sensor* ConstructionWrapper::getSensorByVolume(G4VPhysicalVolume* volume)c
}
void ConstructionWrapper::printSensorEnergies()const{
void GeometryDescriptor::printSensorEnergies()const{
for(const auto& layer: layers){
for(const auto& sensor: layer.sensors){
G4cout << sensor.energy << " ";
src/SteppingAction.cc
+5 -5
View File
@@ -59,15 +59,15 @@ void SteppingAction::UserSteppingAction(const G4Step* step)
//return;
//find the layer through the constructionwrapper that corresponds to the volume
auto cw = DetectorConstruction::getDetectorConstruction()->getConstructionWrapper();
//find the layer through the GeometryDescriptor that corresponds to the volume
auto cw = DetectorConstruction::getDetectorConstruction()->getGeometryDescriptor();
if(cw == nullptr){
G4cout << "ConstructionWrapper not found" << G4endl;
throw std::runtime_error("ConstructionWrapper not found");
G4cout << "GeometryDescriptor not found" << G4endl;
throw std::runtime_error("GeometryDescriptor not found");
}
auto sensor = DetectorConstruction::getDetectorConstruction()->getConstructionWrapper()->getSensorByVolume(volume);
auto sensor = DetectorConstruction::getDetectorConstruction()->getGeometryDescriptor()->getSensorByVolume(volume);
if(sensor == nullptr){
//can simply be not an active volume, so no need to throw an error
//G4cout << "Sensor not found in" << volume->GetName() << G4endl; //DEBUG