Epic restructuring of code to switch to a generator-based style of

event creation.

Now we have vertex generators (that produce initial particles), photon
generators (that create photons to propagate), and a standard data
structure using Python class containers and numpy arrays to hand
around the code.

Also cleaned up some naming of things before they become conventions.

1 files changed, 0 insertions, 165 deletions

diff --git a/g4gen.py b/g4gen.py
deleted file mode 100644
index 6766d42..0000000
--- a/g4gen.py
+++ /dev/null
@@ -1,165 +0,0 @@
-from Geant4 import *
-import g4py.ezgeom
-import g4py.NISTmaterials
-import g4py.ParticleGun
-import pyublas
-import numpy as np
-import pi0
-
-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, seed=None):
-        '''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, ... }.
-           seed: Random number generator seed for HepRandom. If None,
-                 generator is not seeded.
-        '''
-        if seed is not None:
-            HepRandom.setTheSeed(seed)
-
-        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_pi0(self, total_energy, position, direction):
-        '''Use GEANT4 to generate photons produced by a pi0 decay.
-
-        pi0 event will be generated by running two gammas with the appropriate
-        energy and angular distribution.
-
-               total_energy: Total energy of pi0 (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.
-        '''
-        direction = direction / np.linalg.norm(direction)
-        
-        cos_theta_rest = np.random.random_sample() * 2 - 1
-        theta_rest = np.arccos(cos_theta_rest)
-        phi_rest = np.random.random_sample() * 2 * np.pi
-
-        (gamma1_e, gamma1_dir), (gamma2_e, gamma2_dir) = pi0.pi0_decay(energy=total_energy,
-                                                                       direction=direction,
-                                                                       theta=theta_rest,
-                                                                       phi=phi_rest)
-
-        gamma1 = self.generate_photons('gamma', gamma1_e, position, gamma1_dir)
-        gamma2 = self.generate_photons('gamma', gamma2_e, position, gamma2_dir)
-
-        photons = {'subtracks' : [ {'name':'gamma', 'pos': position, 't0' : 0.0,
-                                    'dir' : gamma1_dir, 'total_e' : gamma1_e },
-                                   {'name':'gamma', 'pos': position, 't0' : 0.0,
-                                    'dir' : gamma2_dir, 'total_e' : gamma2_e } ]
-                    }
-        for key in gamma1.keys():
-            photons[key] = np.concatenate((gamma1[key], gamma2[key]))
-                                                                       
-        return photons
-
-    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 = np.zeros(shape=(n,3), dtype=np.float32)
-        pos[:,0] = self.tracking_action.GetX()
-        pos[:,1] = self.tracking_action.GetY()
-        pos[:,2] = self.tracking_action.GetZ()
-
-        dir = np.zeros(shape=(n,3), dtype=np.float32)
-        dir[:,0] = self.tracking_action.GetDirX()
-        dir[:,1] = self.tracking_action.GetDirY()
-        dir[:,2] = self.tracking_action.GetDirZ()
-
-        pol = np.zeros(shape=(n,3), dtype=np.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(np.float32)
-        t0 = self.tracking_action.GetT0().astype(np.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'