displayEvent without .root file

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,
         )
 

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_<Sensor>(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);
 }
 
 

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;
 

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;

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;