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

better plotting

Browse Source
This commit is contained in:
Jan Kieseler
2024-08-23 11:22:50 +02:00
parent 905d5408ed
commit 9038c9521a
5 changed files with 51 additions and 191 deletions
Download Patch File Download Diff File Expand all files Collapse all files
bind/.gitignore
+1
View File
@@ -0,0 +1 @@
minicalo.cpython-38-x86_64-linux-gnu.so
bind/G4Calo.py
+27 -187
View File
@@ -72,177 +72,7 @@ class __G4System(_G4System):
def displayEvent(self, logE = False, renderer=None):
# for loop over all layers
to_plot = []
material_dict = {}
z0=0
# sum up total deposited energy
total_dep_energy = 0
for layer in self.getGeometryDescriptor().getLayers():
for sensor in layer.sensors:
total_dep_energy += sensor.getEnergy()
if total_dep_energy == 0:
print("No energy deposited in calorimeter!")
total_dep_energy = 10**-8 # to avoid division by zero
# loop over materials
for layer in self.getGeometryDescriptor().getLayers():
layer_width = layer.nx * layer.sens_xwidth * 10. # in mm
#
# plot layers
#
layer_hx = layer_width / 2.
layer_hy = layer_width / 2.
layer_z = layer.thickness * 10. # in mm
layer_material = layer.material
# add material to materials if it is not already in there
if layer_material not in material_dict.keys():
material_dict[layer_material] = {'name': layer_material,
'color': col_dict[layer_material],
'showlegend': False,
'flatshading': True,
'opacity': 0.2}
# add legend entry
to_plot.append(go.Mesh3d(x=[None], y=[None], z=[None], i=[0], j=[0], k=[0],
color=material_dict[layer_material]['color'],
showlegend=True, name=layer_material))
to_plot.append(go.Mesh3d(
# 8 vertices of a cube
x = np.array([-1, -1, 1, 1, -1, -1, 1, 1]) * layer_hx,
z = np.array([-1, 1, 1, -1, -1, 1, 1, -1]) * layer_hy,
y = np.array([0, 0, 0, 0, layer_z, layer_z, layer_z, layer_z]) + z0,
**ijk_cube,
**material_dict[layer_material]
))
#
# sensors
#
# if there are sensors in the current layer, add them
if layer.sensors != []:
z = layer.sensors[0].getdz()
corr = 0. #layer.sensors[0].getX() - layer.sensors[0].getdx()/2. + layer_width/2.
# loop over all sensors in current layer and add them to plot
for sensor in layer.sensors:
x_center = sensor.getX() - corr
y_center = sensor.getY() - corr
hwidth = sensor.getdx() /2.
energy = sensor.getEnergy()
use_energy = float(energy / total_dep_energy)
if logE:
raise NotImplementedError
use_energy = np.log(use_energy+1.) # - np.log(total_dep_energy)
to_plot.append(go.Mesh3d(
# 8 vertices of a cube
x = np.array([-1, -1, 1, 1, -1, -1, 1, 1]) * hwidth + x_center,
z = np.array([-1, 1, 1, -1, -1, 1, 1, -1]) * hwidth + y_center,
y = np.array([0, 0, 0, 0, z, z, z, z]) + z0,
**ijk_cube,
flatshading=True,
color='black',
name='Sensor',
opacity= max(0.03, use_energy),
showlegend=False,
))
z0 += layer_z
# add legend entry for sensors
to_plot.append(go.Mesh3d(x=[None], y=[None], z=[None], i=[0], j=[0], k=[0],
color='black', showlegend=True, name='Sensors'))
# add black-white colorbar for sensor hits
to_plot.append(go.Surface(
z=[[0, 0], [0, 0]],
x=[[0, 0], [0, 0]],
y=[[0, 0], [0, 0]],
colorscale=[[0, 'white'], [1, 'black']],
showscale=True,
cmin=0,
cmax=1,
colorbar=dict(
title='Fraction of total deposited Energy',
tickvals=[0, 1],
ticktext=['0', '1'],
ticks='outside',
ticklen=10,
),
))
#
# add red arrow for incoming particle
#
# Define the start and end points of the line
start_point = [0, 0, - z0*0.1]
end_point = [0, 0, - z0*0.25]
# Create the line trace
line_trace = go.Scatter3d(
x=[start_point[0], end_point[0]],
z=[start_point[1], end_point[1]],
y=[start_point[2], end_point[2]],
mode='lines',
line=dict(color='red', width=5),
name='Incoming particle',
showlegend=True,
)
# Calculate the direction vector for the arrow
direction_vector = [(end_point[0] - start_point[0]), (end_point[1] - start_point[1]), (end_point[2] - start_point[2])]
# Create the arrowhead at the start point with the opposite direction
arrowhead_trace = go.Cone(
x=[start_point[0]],
z=[start_point[1]],
y=[start_point[2]],
u=[-direction_vector[0]],
w=[-direction_vector[1]],
v=[-direction_vector[2]],
sizemode='scaled',
sizeref=0.8,
showscale=False,
colorscale='Reds',
opacity=1.0,
anchor='tail',
)
# Create the 3D scatter plot with both traces
to_plot.append(line_trace)
to_plot.append(arrowhead_trace)
#
# finally show plot
#
fig = go.Figure(data=[
*to_plot
])
fig.update_layout(legend=dict(x=0))
# add legend
if renderer is not None:
fig.show(renderer=renderer)
else:
fig.show()
raise NotImplementedError("This method is not implemented anymore. Use the function display_event instead.")
@@ -408,8 +238,8 @@ def display_event(gd : GeometryDescriptor,
if total_dep_energy == 0:
print("No energy deposited in calorimeter!")
total_dep_energy = 10**-8 # to avoid division by zero
# loop over materials
for layer in gd.getLayers():
# loop over layers, invert order
for layer in gd.getLayers()[::-1]:
layer_width = layer.nx * layer.sens_xwidth * 10. # in mm
@@ -418,8 +248,9 @@ def display_event(gd : GeometryDescriptor,
#
layer_hx = layer_width / 2.
layer_hy = layer_width / 2.
layer_z = layer.thickness * 10. # in mm
layer_dz = layer.thickness * 10. # in mm
layer_material = layer.material
z0 = layer.getZ()
# add material to materials if it is not already in there
if layer_material not in material_dict.keys():
@@ -436,7 +267,8 @@ def display_event(gd : GeometryDescriptor,
# 8 vertices of a cube
x = np.array([-1, -1, 1, 1, -1, -1, 1, 1]) * layer_hx,
z = np.array([-1, 1, 1, -1, -1, 1, 1, -1]) * layer_hy,
y = np.array([0, 0, 0, 0, layer_z, layer_z, layer_z, layer_z]) + z0,
#y = np.array([0, 0, 0, 0, layer_dz, layer_dz, layer_dz, layer_dz]) + z0,
y = (layer_dz / 2) * np.array([-1, -1, -1, -1, 1, 1, 1, 1]) + z0,
**ijk_cube,
**material_dict[layer_material]
))
@@ -462,7 +294,8 @@ def display_event(gd : GeometryDescriptor,
# 8 vertices of a cube
x = np.array([-1, -1, 1, 1, -1, -1, 1, 1]) * hwidth + x_center,
z = np.array([-1, 1, 1, -1, -1, 1, 1, -1]) * hwidth + y_center,
y = np.array([0, 0, 0, 0, z, z, z, z]) + z0,
y = (layer_dz / 2) * np.array([-1, -1, -1, -1, 1, 1, 1, 1]) + z0,
#y = np.array([0, 0, 0, 0, z, z, z, z]) + z0,
**ijk_cube,
flatshading=True,
color='black',
@@ -470,9 +303,8 @@ def display_event(gd : GeometryDescriptor,
opacity= max(0.03, use_energy),
showlegend=False,
))
z0 += layer_z
# add legend entry for sensors
to_plot.append(go.Mesh3d(x=[None], y=[None], z=[None], i=[0], j=[0], k=[0],
color='black', showlegend=True, name='Sensors'))
@@ -498,8 +330,8 @@ def display_event(gd : GeometryDescriptor,
# add red arrow for incoming particle
#
# Define the start and end points of the line
start_point = [0, 0, - z0*0.1]
end_point = [0, 0, - z0*0.25]
start_point = [0, 0, gd.getLayers()[0].getZ()-20.]
end_point = [0, 0, gd.getLayers()[0].getZ()-10.]
# Create the line trace
line_trace = go.Scatter3d(
x=[start_point[0], end_point[0]],
@@ -514,12 +346,12 @@ def display_event(gd : GeometryDescriptor,
direction_vector = [(end_point[0] - start_point[0]), (end_point[1] - start_point[1]), (end_point[2] - start_point[2])]
# Create the arrowhead at the start point with the opposite direction
arrowhead_trace = go.Cone(
x=[start_point[0]],
z=[start_point[1]],
y=[start_point[2]],
u=[-direction_vector[0]],
w=[-direction_vector[1]],
v=[-direction_vector[2]],
x=[end_point[0]],
z=[end_point[1]],
y=[end_point[2]],
u=[direction_vector[0]],
w=[direction_vector[1]],
v=[direction_vector[2]],
sizemode='scaled',
sizeref=0.8,
showscale=False,
@@ -533,12 +365,20 @@ def display_event(gd : GeometryDescriptor,
#
# finally show plot
#
fig = go.Figure(data=[
*to_plot
])
fig.update_layout(legend=dict(x=0))
# add legend
#rotate standard view point
# use layer_width as it is constant for all layers and the largest dimension
fig.update_layout(scene_camera=dict(eye=dict(x=1 * 1.5, y=-1.5, z=1 * 1.5)))
#name the axes in the HEP way, so y and z switch names
fig.update_layout(scene=dict(xaxis_title='x [mm]', yaxis_title='z [mm]', zaxis_title='y [mm]'))
if renderer is not None:
fig.show(renderer=renderer)
else:
bind/bindings.cpp
+3
View File
@@ -45,6 +45,9 @@ PYBIND11_MODULE(minicalo, m) {
.def("assignPhysicalVolume", &Layer::assignPhysicalVolume)
.def("unAssign", &Layer::unAssign)
.def_readwrite("sensors", &Layer::sensors)
.def("getX", &Layer::getX)
.def("getY", &Layer::getY)
.def("getZ", &Layer::getZ)
.def(py::pickle(
[](const Layer &l) { // __getstate__
return l.__getstate__();
include/GeometryDescriptor.hh
+19 -4
View File
@@ -69,7 +69,7 @@ public:
class Layer {
public:
Layer() : thickness(0),sens_xwidth(0), sens_ywidth(0), material(""), nx(1), ny(1), isActive(false), physicalVolume(nullptr) {};
Layer() : thickness(0),sens_xwidth(0), sens_ywidth(0), material(""), nx(1), ny(1), isActive(false), physicalVolume(nullptr), position(0,0,0) {};
~Layer() {};
void setThickness(double thickness_cm);
@@ -84,6 +84,18 @@ public:
sensors.clear();
}
double getX() const {
return position.x();
}
double getY() const {
return position.y();
}
double getZ() const {
return position.z();
}
double thickness;
double sens_xwidth;
double sens_ywidth;
@@ -92,17 +104,19 @@ public:
int ny;
bool isActive;
G4VPhysicalVolume* physicalVolume;
std::string name;
G4ThreeVector position;
std::vector<Sensor> sensors;
pybind11::tuple __getstate__() const {
return pybind11::make_tuple(thickness, sens_xwidth, sens_ywidth, material, nx, ny, isActive, sensors);
return pybind11::make_tuple(thickness, sens_xwidth, sens_ywidth, material, nx, ny, isActive, position.x(), position.y(), position.z(), sensors);
}
static Layer __setstate__(pybind11::tuple t) {
if (t.size() != 8) throw std::runtime_error("Invalid state!");
if (t.size() != 11) throw std::runtime_error("Invalid state!");
Layer layer;
layer.thickness = t[0].cast<double>();
layer.sens_xwidth = t[1].cast<double>();
@@ -111,7 +125,8 @@ public:
layer.nx = t[4].cast<int>();
layer.ny = t[5].cast<int>();
layer.isActive = t[6].cast<bool>();
layer.sensors = t[7].cast<std::vector<Sensor>>();
layer.position = G4ThreeVector(t[7].cast<double>(), t[8].cast<double>(), t[9].cast<double>());
layer.sensors = t[10].cast<std::vector<Sensor>>();
layer.physicalVolume = nullptr; // Reset pointer
return layer;
}
src/DetectorConstruction.cc
+1
View File
@@ -290,6 +290,7 @@ G4VPhysicalVolume* DetectorConstruction::DefineVolumes()
layer.sensors.push_back(sensor);//this should now be aligned with copy number
}
}
layer.position = G4ThreeVector(0,0,position);
position += layer.thickness / 2 *cm; //assign the physical volume to the layer