add back EMCal

ref/geant4_simulation/EMCal.py

+__author__ = 'Dowling Wong'
+
+from Geant4 import *
+import geant4_simulation.materials as materials
+
+class EMCal(G4VUserDetectorConstruction):
+    def __init__(self, tower_x=5.535, tower_y=5.535, tower_z=36.96, n_towers_x=72, n_towers_y=36, n_layers=66):
+        """Initialize the EMCal with given or default parameters."""
+        G4VUserDetectorConstruction.__init__(self)
+        self.tower_size_x = tower_x * cm
+        self.tower_size_y = tower_y * cm
+        self.tower_size_z = tower_z * cm
+        self.n_towers_x = n_towers_x
+        self.n_towers_y = n_towers_y
+        self.n_layers = n_layers
+
+        # Calculate world dimensions based on EMCal dimensions
+        self.worldX = self.tower_size_x * self.n_towers_x / 2
+        self.worldY = self.tower_size_y * self.n_towers_y / 2
+        self.worldZ = self.tower_size_z / 2
+
+        # Define step limits for energy deposition readout
+        self.stepLimit = G4UserLimits(4. * mm)
+
+    def Construct(self):
+        global solidWorld, logicalWorld, physicalWorld
+        global solidEMCal, physicalEMCal
+
+        # World volume setup
+        world_size_x = self.tower_size_x * self.n_towers_x
+        world_size_y = self.tower_size_y * self.n_towers_y
+        world_size_z = self.tower_size_z
+
+        solidWorld = G4Box("World", world_size_x * 5, world_size_y * 5, world_size_z * 5)
+        logicalWorld = G4LogicalVolume(solidWorld, G4Material.GetMaterial("G4_AIR"), "World")
+        physicalWorld = G4PVPlacement(None,
+                                      G4ThreeVector(),
+                                      logicalWorld,
+                                      "World",
+                                      None,
+                                      False,
+                                      0)
+
+        # EMCal envelope setup
+        solidEMCal = G4Box("EMCal", world_size_x / 2, world_size_y / 2, world_size_z / 2)
+        self.logicalEMCal = G4LogicalVolume(solidEMCal, G4Material.GetMaterial("G4_AIR"), "EMCal")
+        physicalEMCal = G4PVPlacement(None,
+                                      G4ThreeVector(),
+                                      self.logicalEMCal,
+                                      "EMCal",
+                                      logicalWorld,
+                                      False,
+                                      0)
+
+        # Return the top-level physical world
+        return physicalWorld
+
+
+    def construct_single_tower(self):
+        """Construct the logical volume for a single EMCal tower with layers."""
+        # Define the solid and logical volume for a single tower
+        solidTower = G4Box("Tower", self.tower_size_x / 2, self.tower_size_y / 2, self.tower_size_z / 2)
+        logicalTower = G4LogicalVolume(solidTower, materials.vac, "Tower")
+
+        # Layer thicknesses and materials
+        thickness_layer = self.tower_size_z / float(self.n_layers)
+        thickness_absorber = 1.5 * mm
+        thickness_scint = 4.0 * mm
+        thickness_airgap = thickness_layer - thickness_absorber - thickness_scint
+
+        # Create absorber and scintillator plate solids
+        absorber_solid = G4Box("AbsorberPlate", self.tower_size_x / 2, self.tower_size_y / 2, thickness_absorber / 2)
+        scintillator_solid = G4Box("ScintillatorPlate", self.tower_size_x / 2, self.tower_size_y / 2, thickness_scint / 2)
+        
+        # Logical volumes with appropriate materials
+        absorber_logic = G4LogicalVolume(absorber_solid, materials.Pb, "AbsorberPlate")
+        scintillator_logic = G4LogicalVolume(scintillator_solid, materials.pet, "ScintillatorPlate")
+
+        # Apply step limits
+        absorber_logic.SetUserLimits(self.stepLimit)
+        scintillator_logic.SetUserLimits(self.stepLimit)
+
+        # Position layers within the tower
+        zpos = -self.tower_size_z / 2 + thickness_absorber / 2
+        for i in range(self.n_layers):
+            # Place absorber layer
+            G4PVPlacement(None, G4ThreeVector(0, 0, zpos), absorber_logic, f"Absorber_{i}", logicalTower, False, i)
+            zpos += (thickness_absorber / 2 + thickness_airgap / 2 + thickness_scint / 2)
+
+            # Place scintillator layer
+            G4PVPlacement(None, G4ThreeVector(0, 0, zpos), scintillator_logic, f"Scintillator_{i}", logicalTower, False, i)
+            zpos += (thickness_scint / 2 + thickness_airgap / 2 + thickness_absorber / 2)
+
+        return logicalTower
+
+    def place_towers(self, emcal_logic, tower_logic):
+        """Place towers in a grid within the EMCal volume."""
+        copyNo = 0
+        for ix in range(self.n_towers_x):
+            x_pos = -0.5 * self.n_towers_x * self.tower_size_x + (ix + 0.5) * self.tower_size_x
+            for iy in range(self.n_towers_y):
+                y_pos = -0.5 * self.n_towers_y * self.tower_size_y + (iy + 0.5) * self.tower_size_y
+                G4PVPlacement(None,
+                              G4ThreeVector(x_pos, y_pos, 0),
+                              tower_logic,
+                              f"Tower_{ix}_{iy}",
+                              emcal_logic,
+                              False,
+                              copyNo)
+                copyNo += 1