1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
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'
|
