added cherenkov example and fixed missing optical properties

maureeungaro · maureeungaro · commit e6e68c810db9 · 2026-05-11T16:05:53.000-04:00
diff --git a/api/gconfiguration.py b/api/gconfiguration.py
@@ -217,6 +217,7 @@ def init_variation(self, newVariation):
 			self.matFileName = self.system + "__materials_" + str(self.variation) + ".txt"
 			self.mirFileName = self.system + "__mirrors_" + str(self.variation) + ".txt"
 
+
 	# overwrites any existing geometry file.
 	def initialize_storage(self):
 		print()
@@ -256,7 +257,6 @@ def initialize_storage(self):
 			except OSError as e:
 				sys.exit(f"Error opening file {self.geoFileName}: {e}")
 			try:
-
 				with open(self.matFileName, "w") as file:
 					pass
 			except OSError as e:
diff --git a/examples/meson.build b/examples/meson.build
@@ -17,6 +17,9 @@ examples_map = {
         'b1' : { 'runs' : ['1'], 'variations' : ['default', 'test'] },
         'b2' : { 'runs' : ['1', '11'], 'variations' : ['default', 'alt'] },
         'b3' : { 'runs' : ['1', '11'], 'variations' : ['default', 'alt'] }
+    },
+    'optical' : {
+        'simple_flux' : { 'runs' : ['1'], 'variations' : ['default', 'CO2', 'C4F10'] },
     }
 }
 
@@ -44,7 +47,7 @@ foreach branch : examples_map.keys()
                      args : [example_dir + example + '.py', '-f', format, '-sql', sql_db, '-v', variation], # the sql flag is ignored for ascii format
                      workdir : example_dir,
                      is_parallel : false,
-                     env: test_env,
+                     env : test_env,
                      priority : 10)
                 test(run_test_name_asci,
                      gemc,
@@ -72,7 +75,7 @@ foreach branch : examples_map.keys()
                      args : [example_dir + example + '.py', '-f', format, '-sql', sql_db, '-r', run_value], # the sql flag is ignored for ascii format
                      workdir : example_dir,
                      is_parallel : false,
-                     env: test_env,
+                     env : test_env,
                      priority : 9)
                 test(run_test_name_asci,
                      gemc,
diff --git a/examples/optical/cherenkov/cherenkov.py b/examples/optical/cherenkov/cherenkov.py
@@ -0,0 +1,111 @@
+#!/usr/bin/env python3
+
+from gconfiguration import autogeometry
+from gvolume import GVolume
+from gmaterial import GMaterial
+
+cfg = autogeometry("examples", "cherenkov")
+
+# Using 3 gases with different refractive indices
+variation_material = {
+	"default": "CF4",
+	"CO2": "CO2",
+	"C4F10": "C4F10",
+}
+
+photon_energy = "2.0*eV 3.0*eV 4.0*eV 5.0*eV 6.0*eV"
+
+# Optical properties
+# Refractive index and absorbtion at 2, 3, 4, 5, 6 eV
+rindex = {
+	"CF4": "1.000481811 1.000487934 1.000496771 1.000508615 1.000523881",
+	"CO2": "1.000448007 1.000458006 1.000472449 1.000492454 1.000519671",
+	"C4F10": "1.001306001 1.001329786 1.001364578 1.001412079 1.001474826",
+}
+absorption = {
+	"CF4": "100*m       100*m       100*m       100*m       100*m",
+	"CO2": "50*m        50*m        40*m        30*m        20*m",
+	"C4F10": "6*m        6*m         6*m         5*m         4*m",
+}
+
+material_definitions = {
+	"CF4": {
+		"name": "CF4",
+		"description": "CF4 gas with optical properties",
+		# g/cm3, approximately 20 C and 1 atm
+		"density": 0.003658,
+		"atoms": [("C", 1), ("F", 4)],
+	},
+
+	"CO2": {
+		"name": "CO2",
+		"description": "CO2 gas at 1 atm with optical properties",
+		"density": 0.001842,
+		"atoms": [("C", 1), ("O", 2)],
+	},
+
+	"C4F10": {
+		"name": "C4F10",
+		"description": "C4F10 gas with optical properties",
+		"density": 0.00973,
+		"atoms": [("C", 4), ("F", 10)],
+	},
+}
+
+# common panel geometry
+panels = [
+	("BackplateLeft", -255, 0, 245, 500, 1),
+	("BackplateRight", 255, 0, 245, 500, 2),
+	("BackplateBottom", 0, -255, 500, 245, 3),
+	("BackplateTop", 0, 255, 500, 245, 4),
+]
+
+for variation, material in variation_material.items():
+	cfg.init_variation(variation)
+
+	mat = GMaterial(material_definitions[material]["name"])
+	mat.description = material_definitions[material]["description"]
+	mat.density = material_definitions[material]["density"]
+
+	for element, natoms in material_definitions[material]["atoms"]:
+		mat.addNAtoms(element, natoms)
+
+	mat.photonEnergy = photon_energy
+	mat.indexOfRefraction = rindex[material]
+	mat.absorptionLength = absorption[material]
+	mat.publish(cfg)
+
+	# World
+	world = GVolume("root")
+	world.description = "Air world for Cherenkov validation"
+	world.make_box(600, 600, 600)
+	world.material = "G4_AIR"
+	world.color = "ghostwhite"
+	world.style = 0
+	world.publish(cfg)
+
+	radiator = GVolume("radiator")
+	radiator.mother = "root"
+	radiator.description = f"1 m cube radiator using {material}"
+	radiator.make_box(500, 500, 500)
+	radiator.material = material
+	radiator.color = "66ccff"
+	radiator.opacity = 0.1
+	radiator.publish(cfg)
+
+	for name, x, y, half_x, half_y, panel_id in panels:
+		backplate = GVolume(name)
+		backplate.variation = variation
+		backplate.mother = "radiator"
+		backplate.description = "Downstream optical photon collection panel"
+		backplate.make_box(half_x, half_y, 0.5)
+		backplate.set_position(x, y, 499)
+
+		# Use same gas as the radiator so optical photons are not killed
+		# by crossing into a material without matching optical properties.
+		backplate.material = material
+
+		backplate.color = "lightgray"
+		backplate.digitization = "flux"
+		backplate.set_identifier("detector", 1)
+		backplate.publish(cfg)
diff --git a/examples/optical/cherenkov/cherenkov.yaml b/examples/optical/cherenkov/cherenkov.yaml
@@ -0,0 +1,23 @@
+experiment: examples
+runno: 1
+n: 100
+nthreads: 1
+
+# Required for optical-photon flux hits because optical photons deposit zero energy.
+recordZeroEdep: true
+
+gparticle:
+  - name: e-
+    p: 1000
+    vz: -50
+
+gsystem:
+  - name: cherenkov
+
+gstreamer:
+  - format: csv
+    filename: out
+  - format: root
+    filename: out
+
+phys_list: FTFP_BERT + G4OpticalPhysics
diff --git a/gdetector/gdetectorConstruction.cc b/gdetector/gdetectorConstruction.cc
@@ -20,7 +20,7 @@
 #include "G4RunManager.hh"
 #include "G4UserLimits.hh"
 
-G4ThreadLocal GMagneto* GDetectorConstruction::gmagneto = nullptr;
+G4ThreadLocal GMagneto *GDetectorConstruction::gmagneto = nullptr;
 
 GDetectorConstruction::GDetectorConstruction(std::shared_ptr<GOptions> gopts)
 	: GBase(gopts, GDETECTOR_LOGGER),
@@ -31,7 +31,7 @@ GDetectorConstruction::GDetectorConstruction(std::shared_ptr<GOptions> gopts)
 }
 
 // Builds (or rebuilds) the GEMC world and then the Geant4 world.
-G4VPhysicalVolume* GDetectorConstruction::Construct() {
+G4VPhysicalVolume *GDetectorConstruction::Construct() {
 	log->debug(NORMAL, FUNCTION_NAME);
 
 	// Clean any old geometry.
@@ -52,8 +52,7 @@ G4VPhysicalVolume* GDetectorConstruction::Construct() {
 	if (gsystems.empty()) {
 		log->debug(NORMAL, FUNCTION_NAME, "creating world from options");
 		gworld = std::make_shared<GWorld>(gopt);
-	}
-	else {
+	} else {
 		log->debug(NORMAL, FUNCTION_NAME, "creating world from a gsystem vector of size ", gsystems.size());
 		gworld = std::make_shared<GWorld>(gopt, gsystems);
 	}
@@ -65,8 +64,8 @@ G4VPhysicalVolume* GDetectorConstruction::Construct() {
 
 	// tally with number :
 	log->info(0, "Tally summary: \n - ", gworld->get_number_of_volumes() - 1, " volumes\n - ",
-			  g4world->number_of_volumes(), " geant4 built volumes\n - ",
-			  nsdetectors, " sensitive detectors\n");
+	          g4world->number_of_volumes(), " geant4 built volumes\n - ",
+	          nsdetectors, " sensitive detectors\n");
 
 
 	// Return the physical volume for the ROOT world volume.
@@ -81,28 +80,27 @@ void GDetectorConstruction::ConstructSDandField() {
 
 	// Local cache of sensitive detectors keyed by digitization name.
 	// Multiple volumes can share the same digitization name and therefore reuse one SD instance.
-	std::unordered_map<std::string, GSensitiveDetector*> sensitiveDetectorsMap;
+	std::unordered_map<std::string, GSensitiveDetector *> sensitiveDetectorsMap;
 
 	// Loop over all systems and their volumes.
-	for (const auto& [systemName, gsystemPtr] : *gworld->getSystemsMap()) {
-		for (const auto& [volumeName, gvolumePtr] : gsystemPtr->getGVolumesMap()) {
-			auto const& digitizationName = gvolumePtr->getDigitization();
-			auto const& g4name           = gvolumePtr->getG4Name();
-			auto*       g4volume         = g4world->getG4Volume(g4name)->getLogical();
+	for (const auto &[systemName, gsystemPtr]: *gworld->getSystemsMap()) {
+		for (const auto &[volumeName, gvolumePtr]: gsystemPtr->getGVolumesMap()) {
+			auto const &digitizationName = gvolumePtr->getDigitization();
+			auto const &g4name = gvolumePtr->getG4Name();
+			auto *g4volume = g4world->getG4Volume(g4name)->getLogical();
 
 			// Ensure the Geant4 logical volume exists.
 			// Some GEMC volumes can be "copy-of" another volume; in that case, reuse the
 			// referenced Geant4 logical volume rather than failing.
 			if (g4volume == nullptr) {
 				std::string copyOf = gvolumePtr->getCopyOf();
 				if (copyOf != "" && copyOf != UNINITIALIZEDSTRINGQUANTITY) {
-					auto gsystem      = gvolumePtr->getSystem();
-					auto volume_copy  = gsystem + "/" + copyOf;
+					auto gsystem = gvolumePtr->getSystem();
+					auto volume_copy = gsystem + "/" + copyOf;
 					auto copyG4Volume = g4world->getG4Volume(volume_copy)->getLogical();
-					if (copyG4Volume != nullptr) { g4volume = copyG4Volume; }
-					else {
+					if (copyG4Volume != nullptr) { g4volume = copyG4Volume; } else {
 						log->error(ERR_GVOLUMENOTFOUND, FUNCTION_NAME,
-								   " Logical volume copy <" + volume_copy + "> not found.");
+						           " Logical volume copy <" + volume_copy + "> not found.");
 					}
 				}
 			}
@@ -117,18 +115,17 @@ void GDetectorConstruction::ConstructSDandField() {
 					log->info(2, "Creating new sensitive detector <", digitizationName, "> for volume <", g4name, ">");
 
 					sensitiveDetectorsMap[digitizationName] = new GSensitiveDetector(digitizationName, gopt);
-				}
-				else {
+				} else {
 					log->info(2, "Sensitive detector <", digitizationName,
-							  "> is already created and available for volume <", g4name, ">");
+					          "> is already created and available for volume <", g4name, ">");
 				}
 
 				// Register the volume touchable with the sensitive detector.
 				// The touchable encodes identity and dimension metadata needed by digitization.
-				const auto& vdimensions     = gvolumePtr->getDetectorDimensions();
-				const auto& identity        = gvolumePtr->getGIdentity();
-				const auto& mass            = g4volume->GetMass();
-				auto        this_gtouchable = std::make_shared<
+				const auto &vdimensions = gvolumePtr->getDetectorDimensions();
+				const auto &identity = gvolumePtr->getGIdentity();
+				const auto &mass = g4volume->GetMass();
+				auto this_gtouchable = std::make_shared<
 					GTouchable>(gopt, digitizationName, identity, vdimensions, mass);
 				sensitiveDetectorsMap[digitizationName]->registerGVolumeTouchable(g4name, this_gtouchable);
 
@@ -141,17 +138,17 @@ void GDetectorConstruction::ConstructSDandField() {
 				//g4volume->SetUserLimits(new G4UserLimits(0.1*mm, 0.1*mm));
 
 				log->info(2, "Logical Volume  <" + g4name + "> has been successfully assigned to SD.",
-						  sensitiveDetectorsMap[digitizationName]);
+				          sensitiveDetectorsMap[digitizationName]);
 			}
 
 			// Process electromagnetic fields.
 			// If a volume declares an EM field, ensure the field container exists and install
 			// a per-volume field manager configured by the named field map.
-			const auto& field_name = gvolumePtr->getEMField();
+			const auto &field_name = gvolumePtr->getEMField();
 			if (field_name != "" && field_name != UNINITIALIZEDSTRINGQUANTITY) {
 				if (gmagneto == nullptr) { gmagneto = new GMagneto(gopt); }
 				log->info(2, "Volume <", volumeName, "> has field: <", field_name,
-						  ">. Looking into field map definitions.");
+				          ">. Looking into field map definitions.");
 				log->info(2, "Setting field manager for volume <", g4name, "> with field <", field_name, ">");
 				g4world->setFieldManagerForVolume(g4name, gmagneto->getFieldMgr(field_name).get(), true);
 			}
@@ -164,22 +161,22 @@ void GDetectorConstruction::ConstructSDandField() {
 
 	// Bind each digitization routine to its corresponding sensitive detector.
 	const auto sdetectors = gworld->getSensitiveDetectorsList();
-	for (auto& sdname : sdetectors) {
-		auto   digitization_routine = digitization_routines_map->at(sdname);
-		double maxStep              = digitization_routine->readoutSpecs->getMaxStep();
+	for (auto &sdname: sdetectors) {
+		auto digitization_routine = digitization_routines_map->at(sdname);
+		double maxStep = digitization_routine->readoutSpecs->getMaxStep();
 
 		sensitiveDetectorsMap[sdname]->assign_digi_routine(digitization_routine);
 		log->info(1, "Digitization routine <" + sdname + "> has been successfully assigned to SD.",
-				  sensitiveDetectorsMap[sdname]);
+		          sensitiveDetectorsMap[sdname]);
 
 		// Loop over all systems and their volumes.
 		// and assign max step to the corresponding logical volume
-		for (const auto& [systemName, gsystemPtr] : *gworld->getSystemsMap()) {
-			for (const auto& [volumeName, gvolumePtr] : gsystemPtr->getGVolumesMap()) {
-				auto const& digitizationName = gvolumePtr->getDigitization();
+		for (const auto &[systemName, gsystemPtr]: *gworld->getSystemsMap()) {
+			for (const auto &[volumeName, gvolumePtr]: gsystemPtr->getGVolumesMap()) {
+				auto const &digitizationName = gvolumePtr->getDigitization();
 				if (digitizationName == sdname) {
-					auto const& g4name   = gvolumePtr->getG4Name();
-					auto*       g4volume = g4world->getG4Volume(g4name)->getLogical();
+					auto const &g4name = gvolumePtr->getG4Name();
+					auto *g4volume = g4world->getG4Volume(g4name)->getLogical();
 
 					// g4volume->SetUserLimits(new G4UserLimits(maxStep, maxStep)); // this will also kill track cause
 					// the second argument is max track length
@@ -189,28 +186,24 @@ void GDetectorConstruction::ConstructSDandField() {
 				}
 			}
 		}
-
 	}
 }
 
 // Loads (or dynamically resolves) the digitization routine for each sensitive detector.
 void GDetectorConstruction::loadDigitizationPlugins() {
 	const auto sdetectors = gworld->getSensitiveDetectorsList();
 
-	for (auto& sdname : sdetectors) {
+	for (auto &sdname: sdetectors) {
 		if (sdname == FLUXNAME) {
 			log->info(1, "Loading flux digitization plugin for routine <" + sdname + ">");
 			digitization_routines_map->emplace(sdname, std::make_shared<GFluxDigitization>(gopt));
-		}
-		else if (sdname == COUNTERNAME) {
+		} else if (sdname == COUNTERNAME) {
 			log->info(1, "Loading particle counter digitization plugin for routine <" + sdname + ">");
 			digitization_routines_map->emplace(sdname, std::make_shared<GParticleCounterDigitization>(gopt));
-		}
-		else if (sdname == DOSIMETERNAME) {
+		} else if (sdname == DOSIMETERNAME) {
 			log->info(1, "Loading dosimeter digitization plugin for routine <" + sdname + ">");
 			digitization_routines_map->emplace(sdname, std::make_shared<GDosimeterDigitization>(gopt));
-		}
-		else {
+		} else {
 			// if it's not in the map already, add it
 			log->info(0, "Loading new digitization plugin for routine <" + sdname + ">");
 			digitization_routines_map->emplace(sdname, gdynamicdigitization::load_dynamicRoutine(sdname, gopt));
@@ -221,8 +214,7 @@ void GDetectorConstruction::loadDigitizationPlugins() {
 
 		if (digitization_routines_map->at(sdname)->defineReadoutSpecs()) {
 			log->info(1, "Digitization routine <" + sdname + "> has been successfully defined.");
-		}
-		else { log->error(ERR_DEFINESPECFAIL, "defineReadoutSpecs failure for <" + sdname + ">"); }
+		} else { log->error(ERR_DEFINESPECFAIL, "defineReadoutSpecs failure for <" + sdname + ">"); }
 	}
 }
 
@@ -242,6 +234,5 @@ void GDetectorConstruction::reload_geometry(SystemList sl) {
 	if (rm) {
 		rm->DefineWorldVolume(Construct());
 		ConstructSDandField();
-	}
-	else { log->error(1, "GDetectorConstruction::reload_geometry", "Geant4 Run manager not found."); }
+	} else { log->error(1, "GDetectorConstruction::reload_geometry", "Geant4 Run manager not found."); }
 }
diff --git a/releases/0.2.md b/releases/0.2.md
diff --git a/releases/dev.md b/releases/dev.md
