diff --git a/bind/G4Calo.py b/bind/G4Calo.py index 3907630..f2f2ed6 100644 --- a/bind/G4Calo.py +++ b/bind/G4Calo.py @@ -36,29 +36,34 @@ class G4System(_G4System): df = ttree["Hits;1"].arrays(library="pd") return df - def displayEvent(self, particleSpec, minEnergy_GeV, maxEnergy_GeV=-1, sensor_width=50): - if not isinstance(sensor_width, np.ndarray): - sensor_width = np.array([sensor_width]) - + def displayEvent(self, gd): # for loop over all layers - event = self.run_batch(1, particleSpec, minEnergy_GeV, maxEnergy_GeV) to_plot = [] material_dict = {} - z0=0 - # loop over materials - for layer_i, layer in enumerate(self.cw.getLayers()): - #layer_i = len(self.cw.getLayers()) - layer_i -1 - print(layer_i) + # sum up total deposited energy + total_dep_energy = 0 + for layer in gd.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 gd.getLayers(): + + layer_width = layer.nx * layer.sens_xwidth + # # plot layers # - layer_hx = sensor_width / 2. - layer_hy = sensor_width / 2. + layer_hx = layer_width / 2. + layer_hy = layer_width / 2. layer_z = layer.thickness layer_material = layer.material - print(layer_material) # add material to materials if it is not already in there if layer_material not in material_dict.keys(): @@ -88,17 +93,16 @@ class G4System(_G4System): # # if there are sensors in the current layer, add them - if layer_i in event['sensor_layer'].to_numpy(): - is_in_layer = event['sensor_layer'].to_numpy() == layer_i - print(event['sensor_energy'].to_numpy()[is_in_layer].sum()) - z = event['sensor_dz'].to_numpy()[is_in_layer] - n_sensors = len(z) - z=z[0] - - xy_centers, hwidth = calculate_sensor_centers(n_sensors, sensor_width) - + if layer.sensors != []: + z = layer.sensors[0].getdz() + corr = 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 (x_center, y_center), energy in zip(xy_centers, event['sensor_energy'].to_numpy()[is_in_layer]): + for sensor in layer.sensors: + x_center = sensor.getX() - corr + y_center = sensor.getY() - corr + hwidth = sensor.getdx() /2. + energy = sensor.getEnergy() + to_plot.append(go.Mesh3d( # 8 vertices of a cube x = np.array([-1, -1, 1, 1, -1, -1, 1, 1]) * hwidth + x_center, @@ -108,7 +112,7 @@ class G4System(_G4System): flatshading=True, color='black', name='Sensor', - opacity= max(0.03, float(energy / event['total_dep_energy'].to_numpy())), + opacity= max(0.03, float(energy / total_dep_energy)), showlegend=False, )) @@ -155,7 +159,7 @@ class G4System(_G4System): y=[start_point[2], end_point[2]], mode='lines', line=dict(color='red', width=5), - name='Incoming ' + particleSpec + ' @{} GeV'.format(minEnergy_GeV), + name='Incoming particle', showlegend=True, ) diff --git a/bind/bindings.cpp b/bind/bindings.cpp index 304eaa3..587449b 100644 --- a/bind/bindings.cpp +++ b/bind/bindings.cpp @@ -32,9 +32,7 @@ void makeG4System(M &m, std::string name) .def("applyUICommand", &G4System::applyUICommand, py::arg("command")) .def("displayEvent", &G4System::displayEvent) .def("printMaterial", &G4System::printMaterial, py::arg("name")) - .def("check", &G4System::check) - .def_readwrite("cw", &G4System::cw); - + .def("check", &G4System::check); } // create bindings for Layer class @@ -57,7 +55,13 @@ void makeSensor(M &m, std::string name){ py::class_(m, name.data()).def(py::init()) .def("getEnergy", &Sensor::getEnergy) .def("getPos", &Sensor::getPos) - .def("getSize", &Sensor::getSize); + .def("getSize", &Sensor::getSize) + .def("getX", &Sensor::getX) + .def("getY", &Sensor::getY) + .def("getZ", &Sensor::getZ) + .def("getdx", &Sensor::getdx) + .def("getdy", &Sensor::getdy) + .def("getdz", &Sensor::getdz); } diff --git a/include/ConstructionWrapper.hh b/include/ConstructionWrapper.hh index 0ad2378..7a4ed49 100644 --- a/include/ConstructionWrapper.hh +++ b/include/ConstructionWrapper.hh @@ -24,7 +24,30 @@ public: return size; } + double getX()const{ + return position.x(); + } + + double getY()const{ + return position.y(); + } + + double getZ()const{ + return position.z(); + } + + double getdx()const{ + return size.x(); + } + double getdy()const{ + return size.y(); + } + + double getdz()const{ + return size.z(); + } + G4ThreeVector position; G4ThreeVector size; mutable G4double energy; @@ -50,6 +73,7 @@ public: int nx; int ny; bool isActive; + G4VPhysicalVolume* physicalVolume; std::string name; diff --git a/include/G4System.hh b/include/G4System.hh index 9f59b67..3e85711 100644 --- a/include/G4System.hh +++ b/include/G4System.hh @@ -43,8 +43,6 @@ void applyUICommand(const std::string& command){ UImanager->ApplyCommand(command); } -ConstructionWrapper cw; - void displayEvent()const{}; //just a placeholder, this will be implemented in python void printMaterial(const std::string& name)const; diff --git a/src/G4System.cc b/src/G4System.cc index 2e82f0b..bb43651 100644 --- a/src/G4System.cc +++ b/src/G4System.cc @@ -33,9 +33,6 @@ void G4System::init(ConstructionWrapper &CW){ //delete actionInitialization; } - // save layers of CW to G4System class - G4System::cw = CW; - G4String session; ui = nullptr;