1 files changed, 125 insertions, 0 deletions
diff --git a/g4gen.py b/g4gen.py
new file mode 100644
index 0000000..409570f
--- /dev/null
+++ b/g4gen.py
@@ -0,0 +1,125 @@
+from Geant4 import *
+import g4py.ezgeom
+import g4py.NISTmaterials
+import g4py.ParticleGun
+import pyublas
+import numpy
+
+try:
+    import G4chroma
+except:
+    # Try building the module
+    import subprocess
+    import sys, os
+    module_dir = os.path.split(os.path.realpath(__file__))[0]
+    print >>sys.stderr, 'Compiling G4chroma.so...'
+    retcode = subprocess.call('g++ -o \'%s/G4chroma.so\' -shared \'%s/G4chroma.cc\' -fPIC `geant4-config --cflags --libs` `python-config --cflags --libs --ldflags` -lboost_python' % (module_dir, module_dir), shell=True)
+    assert retcode == 0
+    import G4chroma
+
+class G4Generator(object):
+    def __init__(self, material):
+        '''Create generator to produce photons inside the specified material.
+
+           material: chroma.geometry.Material object with density, 
+                     composition dict and refractive_index.
+
+                     composition dictionary should be 
+                        { element_symbol : fraction_by_weight, ... }.
+        '''
+        g4py.NISTmaterials.Construct()
+        g4py.ezgeom.Construct()
+        self.physics_list = G4chroma.ChromaPhysicsList()
+        gRunManager.SetUserInitialization(self.physics_list)
+        self.particle_gun = g4py.ParticleGun.Construct()
+        
+        self.world_material = self.create_g4material(material)
+        g4py.ezgeom.SetWorldMaterial(self.world_material)
+
+        self.world = g4py.ezgeom.G4EzVolume('world')
+        self.world.CreateBoxVolume(self.world_material, 100*m, 100*m, 100*m)
+        self.world.PlaceIt(G4ThreeVector(0,0,0))
+
+        self.tracking_action = G4chroma.PhotonTrackingAction()
+        gRunManager.SetUserAction(self.tracking_action)
+        gRunManager.Initialize()
+
+    def create_g4material(self, material):
+        g4material = G4Material('world_material', material.density * g / cm3,
+                                len(material.composition))
+
+        # Add elements
+        for element_name, element_frac_by_weight in material.composition.items():
+            g4material.AddElement(G4Element.GetElement(element_name, True),
+                                  element_frac_by_weight)
+
+        # Set index of refraction
+        prop_table = G4MaterialPropertiesTable()
+        # Reverse entries so they are in ascending energy order rather
+        # than wavelength
+        energy = list((2*pi*hbarc / (material.refractive_index[::-1,0] * nanometer)).astype(float))
+        values = list(material.refractive_index[::-1, 1].astype(float))
+        prop_table.AddProperty('RINDEX', energy, values)
+
+        # Load properties
+        g4material.SetMaterialPropertiesTable(prop_table)
+        return g4material
+
+
+    def generate_photons(self, particle_name, total_energy, position, direction):
+        '''Use GEANT4 to generate photons produced by the given particle.
+           
+               particle_name: GEANT4 name of particle.  'e-', 'mu-', etc
+               total_energy: Total energy of particle (incl rest mass) in MeV
+               position: 3-tuple of position of particle in global coordinates
+               direction: 3-tuple direction vector.  
+                          Does not have to be normalized.
+        '''
+        self.particle_gun.SetParticleByName(particle_name)
+        self.particle_gun.SetParticleEnergy(total_energy * MeV)
+        self.particle_gun.SetParticlePosition(G4ThreeVector(*position))
+        self.particle_gun.SetParticleMomentumDirection(G4ThreeVector(*direction).unit())
+
+        self.tracking_action.Clear()
+
+        gRunManager.BeamOn(1)
+        n = self.tracking_action.GetNumPhotons()        
+        pos = numpy.zeros(shape=(n,3), dtype=numpy.float32)
+        pos[:,0] = self.tracking_action.GetX()
+        pos[:,1] = self.tracking_action.GetY()
+        pos[:,2] = self.tracking_action.GetZ()
+
+        dir = numpy.zeros(shape=(n,3), dtype=numpy.float32)
+        dir[:,0] = self.tracking_action.GetDirX()
+        dir[:,1] = self.tracking_action.GetDirY()
+        dir[:,2] = self.tracking_action.GetDirZ()
+
+        pol = numpy.zeros(shape=(n,3), dtype=numpy.float32)
+        pol[:,0] = self.tracking_action.GetPolX()
+        pol[:,1] = self.tracking_action.GetPolY()
+        pol[:,2] = self.tracking_action.GetPolZ()
+        
+        wavelength = self.tracking_action.GetWavelength().astype(numpy.float32)
+        t0 = self.tracking_action.GetT0().astype(numpy.float32)
+        return { 'pos' : pos, 
+                 'dir' : dir, 
+                 'pol' : pol,
+                 'wavelength' : wavelength, 
+                 't0' : t0 }
+
+if __name__ == '__main__':
+    import time
+    import optics
+    gen = G4Generator(optics.water)
+    # prime things
+    gen.generate_photons('e-', 1, (0,0,0), (1,0,0))
+
+    start = time.time()   
+    n = 0
+    for i in xrange(100):
+        photons = gen.generate_photons('mu-', 700, (0,0,0), (1,0,0))
+        n += len(photons['t0'])
+        print photons['pos'][0].min(), photons['pos'][0].max()
+    stop = time.time()
+    print stop - start, 'sec'
+    print n / (stop-start), 'photons/sec'