dd4hepplugins/examples/test_raindrop_dd4hep_cpu.py

0001 #!/usr/bin/env python3
0002 """
0003 CPU-only raindrop optical photon test via DD4hep.
0004
0005 Uses standard G4OpticalPhysics (no eic-opticks GPU plugins).
0006 Optical photons are tracked on CPU by Geant4 and collected as
0007 Geant4Tracker::Hit via the sensitive detector.
0008
0009 Compare hit counts with test_raindrop_dd4hep_gpu.py to validate.
0010
0011 Prerequisites:
0012   - Spack environment activated (ROOT, DD4hep on PYTHONPATH/LD_LIBRARY_PATH)
0013   - DD4hepINSTALL set (for elements.xml lookup)
0014   - libRaindropGeo.so on DD4HEP_LIBRARY_PATH
0015 """
0016 import os
0017 import sys
0018
0019 import cppyy
0020 import DDG4
0021 from g4units import MeV
0022
0023 cppyy.include("G4OpticalParameters.hh")
0024
0025 _SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
0026
0027
0028 def run():
0029     kernel = DDG4.Kernel()
0030
0031     compact_file = os.path.join(_SCRIPT_DIR, "geometry", "raindrop_dd4hep.xml")
0032     if not os.path.exists(compact_file):
0033         print(f"ERROR: Compact file not found: {compact_file}")
0034         sys.exit(1)
0035
0036     if "DD4hepINSTALL" not in os.environ:
0037         print("ERROR: DD4hepINSTALL not set. Activate the spack environment first.")
0038         sys.exit(1)
0039
0040     print(f"Loading geometry: {compact_file}")
0041     kernel.loadGeometry(str("file:" + compact_file))
0042
0043     geant4 = DDG4.Geant4(kernel)
0044     geant4.printDetectors()
0045     geant4.setupUI(typ='tcsh', vis=False, ui=False)
0046
0047     # Same particle gun as the GPU test: e- at 10 MeV inside water box
0048     geant4.setupGun(
0049         "Gun",
0050         particle='e-',
0051         energy=10 * MeV,
0052         position=(0, 0, 0),
0053         isotrop=False,
0054         direction=(1.0, 0.0, 0.0),
0055         multiplicity=1,
0056     )
0057
0058     # Register optical photon sensitive detector on lead container.
0059     # Use Geant4OpticalTrackerAction (not default Geant4SimpleTrackerAction
0060     # which has a 1 keV min energy deposit cut -- optical photons are ~eV).
0061     seq, act = geant4.setupDetector('Raindrop', 'Geant4OpticalTrackerAction')
0062     filt = DDG4.Filter(kernel, 'ParticleSelectFilter/OpticalPhotonSelector')
0063     filt.particle = 'opticalphoton'
0064     seq.adopt(filt)
0065
0066     # Physics: QGSP_BERT + standard Geant4 optical physics (CPU tracking)
0067     geant4.setupPhysics('QGSP_BERT')
0068     ph = DDG4.PhysicsList(kernel, 'Geant4PhysicsList/OpticalPhys')
0069     ph.addPhysicsConstructor(str('G4OpticalPhysics'))
0070     kernel.physicsList().adopt(ph)
0071
0072     # Save hits to ROOT file
0073     output_file = "/tmp/raindrop_hits_cpu.root"
0074     geant4.setupROOTOutput('RootOutput', output_file)
0075
0076     # Run 2 events (same as GPU test)
0077     kernel.NumEvents = 2
0078     kernel.configure()
0079     kernel.initialize()
0080
0081     # Enable boundary SD invocation so photons detected at the air/lead
0082     # surface (EFFICIENCY=1.0) trigger hits in the sensitive lead volume.
0083     # Geant4 11 defaults BoundaryInvokeSD to false.
0084     from cppyy.gbl import G4OpticalParameters
0085     G4OpticalParameters.Instance().SetBoundaryInvokeSD(True)
0086
0087     kernel.run()
0088     kernel.terminate()
0089
0090     # Verify hits
0091     verify_hits(output_file)
0092
0093
0094 def verify_hits(root_file):
0095     """Read back detected photon hits and print summary."""
0096     import ROOT
0097     f = ROOT.TFile.Open(root_file)
0098     if not f or f.IsZombie():
0099         print(f"ERROR: Cannot open {root_file}")
0100         return
0101
0102     tree = f.Get("EVENT")
0103     if not tree:
0104         print("ERROR: No EVENT tree in ROOT file")
0105         f.Close()
0106         return
0107
0108     print(f"\n{'='*60}")
0109     print(f"  CPU hit verification from {root_file}")
0110     print(f"{'='*60}")
0111     print(f"  Events in file: {tree.GetEntries()}")
0112
0113     for evt_idx in range(tree.GetEntries()):
0114         tree.GetEntry(evt_idx)
0115         hits = tree.RaindropHits
0116         nhits = hits.size()
0117
0118         print(f"\n  Event {evt_idx}:")
0119         print(f"    Detected photons: {nhits}")
0120         for i in range(min(3, nhits)):
0121             h = hits[i]
0122             pos = h.position
0123             mom = h.momentum
0124             print(f"      [{i}] pos=({pos.x():.1f}, {pos.y():.1f}, {pos.z():.1f}) "
0125                   f"mom=({mom.x():.3f}, {mom.y():.3f}, {mom.z():.3f}) "
0126                   f"time={h.truth.time:.2f}ns")
0127
0128     f.Close()
0129     print(f"\n  CPU test complete")
0130     print(f"{'='*60}\n")
0131
0132
0133 if __name__ == "__main__":
0134     run()