diff --git a/bind/.gitignore b/bind/.gitignore new file mode 100644 index 0000000..e443899 --- /dev/null +++ b/bind/.gitignore @@ -0,0 +1 @@ +minicalo.cpython-38-x86_64-linux-gnu.so \ No newline at end of file diff --git a/bind/G4Calo.py b/bind/G4Calo.py index 45af3a1..b4663fd 100644 --- a/bind/G4Calo.py +++ b/bind/G4Calo.py @@ -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: diff --git a/bind/bindings.cpp b/bind/bindings.cpp index 4079dc9..3c767f6 100644 --- a/bind/bindings.cpp +++ b/bind/bindings.cpp @@ -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__(); diff --git a/include/GeometryDescriptor.hh b/include/GeometryDescriptor.hh index 9bf9bd2..eca710b 100644 --- a/include/GeometryDescriptor.hh +++ b/include/GeometryDescriptor.hh @@ -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 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(); layer.sens_xwidth = t[1].cast(); @@ -111,7 +125,8 @@ public: layer.nx = t[4].cast(); layer.ny = t[5].cast(); layer.isActive = t[6].cast(); - layer.sensors = t[7].cast>(); + layer.position = G4ThreeVector(t[7].cast(), t[8].cast(), t[9].cast()); + layer.sensors = t[10].cast>(); layer.physicalVolume = nullptr; // Reset pointer return layer; } diff --git a/src/DetectorConstruction.cc b/src/DetectorConstruction.cc index 1a079fb..71f1c17 100644 --- a/src/DetectorConstruction.cc +++ b/src/DetectorConstruction.cc @@ -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