NaI Geant4

import ROOT
from ROOT import TCanvas
import Geant4
from Geant4 import cm, mm, MeV, deg, rad
import g4py.EMSTDpl
import g4py.NISTmaterials
import g4py.ezgeom
import g4py.ParticleGun
from Geant4 import G4ThreeVector

import math
from math import sin, cos, atan2
import sys

############################################################
#                        CONSTANTS                         #
############################################################
PI = 3.1415926
OUTPUT_DATA = True
# Number of particles to emit:
PARTICLE_COUNT = 10000
# Number of emitted particles which will be slowly displayed
REVIEWED_CASES = 0
# Number of emitted particles which will be displayed
SHOWN_CASES = 0
# Histogram parameters
ENERGY_BIN_COUNT = 100
ANGLE_BIN_COUNT  = 18
ENERGY_BIAS = 20 * MeV
WORLD_SIZE = 200.0 * cm

############################################################
#                 DETECTOR CLASS DEFINITION                #
############################################################
hits = 0

# electron/gamma standard EM physics list
g4py.EMSTDpl.Construct()

class MySD(Geant4.G4VSensitiveDetector):
    "My sensitive detector"
    def __init__(self):
        Geant4.G4VSensitiveDetector.__init__(self, "MySD")
        self.eventEnergy = 0.0
    def ProcessHits(self, step, rohist):
        energy = step.GetTotalEnergyDeposit()/MeV
        self.eventEnergy += energy

        point = step.GetPreStepPoint().GetPosition()
        self.enX += energy * point.getX() / cm
        self.enY += energy * point.getY() / cm
        self.enZ += energy * point.getZ() / cm

class MyEventAction(Geant4.G4UserEventAction):
    def __init__(self, sd, hist, hist2):
        Geant4.G4UserEventAction.__init__(self)
        self.sd    = sd
        self.hist  = hist
        self.hist2 = hist2
        self.count = 0

        self.hist2.GetXaxis().SetTitle("Energy (MeV)")
        self.hist2.GetYaxis().SetTitle("Angle (degrees)")
    def BeginOfEventAction(self, event):
        self.sd.enX = 0
        self.sd.enY = 0
        self.sd.enZ = 0

        self.count += 1
        if OUTPUT_DATA == True:
            print '==== start of event #', self.count
        self.sd.eventEnergy = 0
    def EndOfEventAction(self, event):
        if self.sd.eventEnergy > ENERGY_BIAS:
            global hits
            hits += 1
        self.hist.Fill(self.sd.eventEnergy)
        self.theta = atan2((self.sd.enX**2 + self.sd.enY**2)**0.5, self.sd.enZ)
        if self.sd.negative == True and self.theta < PI/2.0:
            self.theta = PI - self.theta
        self.hist2.Fill(self.sd.eventEnergy, self.theta * 180.0/PI)
        if OUTPUT_DATA == True:
            print '==== end of event #', self.count, ', energy=', self.sd.eventEnergy

############################################################
#                  INITIALIZING WORLD                      #
############################################################

g4py.ezgeom.ResizeWorld(WORLD_SIZE, WORLD_SIZE, WORLD_SIZE)

g4py.NISTmaterials.Construct()
g4py.ezgeom.Construct()

# fill the world with air
air = Geant4.G4Material.GetMaterial("G4_AIR")
g4py.ezgeom.SetWorldMaterial(air)

# Sodium Iodide is the official name for NaI
NaI = Geant4.G4Material.GetMaterial("G4_SODIUM_IODIDE")
target = g4py.ezgeom.G4EzVolume("Sodium Iodide Sphere")

target.CreateShpereVolume(NaI, 25.0*cm, 85.0*cm, 0.0*deg, 360.0*deg, 18.0*deg, 144.0*deg)

target.PlaceIt(G4ThreeVector(0.0, 0.0, 0.0))

############################################################
#                INITIALIZING DETECTOR                     #
############################################################

detector = MySD()
target.SetSensitiveDetector(detector)

############################################################
#               INITIALIZING HISTOGRAMS                    #
############################################################

canvas = TCanvas("SomeName", "SomeTitle", 800, 600)
canvas.Divide(2,1)

hist  = ROOT.TH1D("hist", "Energy", ENERGY_BIN_COUNT, 0, 800.0)
hist2 = ROOT.TH2D("hist2", "Energy2", ENERGY_BIN_COUNT, 0, 800.0, ANGLE_BIN_COUNT, 0, 180.05)
uaction = MyEventAction(detector, hist, hist2)
Geant4.gRunManager.Initialize()
Geant4.gRunManager.SetUserAction(uaction)

Geant4.gApplyUICommand('/control/execute visinit.mac')
Geant4.gApplyUICommand('/vis/enable true')

pg = g4py.ParticleGun.Construct()

############################################################
#                     OTHER STUFF                          #
############################################################

pg.SetParticleByName("gamma")
pg.SetParticlePosition( Geant4.G4ThreeVector(0.0, 0.0, 0.0) )
pg.SetParticleEnergy(500.0*MeV)
for i in range(0,PARTICLE_COUNT):
    phi = ROOT.gRandom.Uniform(0.0, 2*PI)
    theta = math.acos(ROOT.gRandom.Uniform(-1.0, 1.0))

    if (theta >= PI/2.0):
        uaction.sd.negative = True
    else:
        uaction.sd.negative = False

    vx = sin(theta)*cos(phi)
    vy = sin(theta)*sin(phi)
    vz = cos(theta)
    pg.SetParticleMomentumDirection( Geant4.G4ThreeVector(vx, vy, vz) )
    Geant4.gRunManager.BeamOn(1)
    if i < REVIEWED_CASES:
        sys.stdin.readline()
    if i == SHOWN_CASES:
        Geant4.gApplyUICommand('/vis/enable false')

canvas.cd(1)
hist.Draw()
canvas.cd(2)
hist2.Draw("surf")
canvas.Update()
print 'Efficiency: ', (float(hits) / PARTICLE_COUNT) * 100.0, '%'
sys.stdin.readline()