Igor's Geant4

import ROOT
import sys
import math
import ROOT
import Geant4
from Geant4 import cm, mm, MeV
import g4py.EMSTDpl
import g4py.NISTmaterials
import g4py.ezgeom
import g4py.ParticleGun

totalNumb = 10000
absorbNumb = 0
gammaEnergy = 6.0

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

class MyEventAction(Geant4.G4UserEventAction):
    def __init__(self, sd1, sd2, hist1, hist2, hist3):
        Geant4.G4UserEventAction.__init__(self)
        self.sd1 = sd1
        self.sd2 = sd2
        self.hist1 = hist1
        self.hist2 = hist2
        self.hist3 = hist3
    def BeginOfEventAction(self, event):

        self.sd1.eventEnergy = 0
        self.sd2.eventEnergy = 0
    def EndOfEventAction(self, event):
        self.hist1.Fill(self.sd1.eventEnergy)
        self.hist2.Fill(self.sd2.eventEnergy)
        global absorbNumb
        if (self.sd2.eventEnergy > gammaEnergy * 0.99):
            absorbNumb = absorbNumb + 1
        if(self.sd1.eventEnergy or self.sd2.eventEnergy):
            self.hist3.Fill(self.sd1.eventEnergy / (self.sd1.eventEnergy + self.sd2.eventEnergy))

c= ROOT.TCanvas("c")
c.Clear()
c.Divide(2,2)
c.cd(1)
# electron/gamma standard EM physics list
g4py.EMSTDpl.Construct()
g4py.NISTmaterials.Construct()
#print Geant4.gMaterialTable

g4py.ezgeom.Construct()
# fill the world with air
air = Geant4.G4Material.GetMaterial("G4_AIR")
g4py.ezgeom.SetWorldMaterial(air)
g4py.ezgeom.ResizeWorld(100.0*cm,100.0*cm,100.0*cm)

Ge = Geant4.G4Material.GetMaterial("G4_Ge")
target1 = g4py.ezgeom.G4EzVolume("Detector")
target1.CreateBoxVolume(Ge, 5.0*cm, 5.0*cm, 5.0*cm)
target1.PlaceIt(Geant4.G4ThreeVector(0.0*cm, 0.0*cm, 5*cm + 2.5*cm))

Pb = Geant4.G4Material.GetMaterial("G4_Pb")
target2 = g4py.ezgeom.G4EzVolume("Absorber")
target2.CreateBoxVolume(Pb, 20.0*cm, 20.0*cm, 20.0*cm)
target2.PlaceIt(Geant4.G4ThreeVector(0.0, 0.0, 5.*cm + 5.*cm + 10.*cm))



detector1 = MySD()
target1.SetSensitiveDetector(detector1)

detector2 = MySD()
target2.SetSensitiveDetector(detector2)


hist1 = ROOT.TH1D("hist1", "Detector", 500, -1, 7)
hist2 = ROOT.TH1D("hist2", "Absorber", 500, -1, 7)
hist3 = ROOT.TH1D("hist3", "Proprtion", 500, -0.5, 1.5)

uaction = MyEventAction(detector1,detector2, hist1,hist2,hist3)
Geant4.gRunManager.Initialize()
Geant4.gRunManager.SetUserAction(uaction)

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


pg = g4py.ParticleGun.Construct()
pg.SetParticleByName("gamma")
pg.SetParticlePosition( Geant4.G4ThreeVector(0.0*cm,0.0*cm,0.*cm) )

pg.SetParticleEnergy(6.0*MeV)


for i in range(0,10000):
    if(i==10):
        Geant4.gApplyUICommand("/vis/enable false")
    phi   = ROOT.gRandom.Uniform(0.0, 2*math.pi)
    theta = ROOT.gRandom.Gaus(0.0, 0.001*math.pi/180)
    vx = math.sin(theta)*math.cos(phi)
    vy = math.sin(theta)*math.sin(phi)
    vz = math.cos(theta)
    pg.SetParticleMomentumDirection( Geant4.G4ThreeVector(vx, vy, vz) )
    Geant4.gRunManager.BeamOn(1)
    if(i<=10):
        print i
        raw_input('--> ')

hist1.Draw()
c.cd(2)
hist2.Draw()
c.cd(3)
hist3.Draw()

print "absorption efficiency: ", float(absorbNumb) / totalNumb

print "=============="
print "=============="
print "=============="
print "=============="
print "=============="
raw_input('--> ')
raw_input('--> ')
raw_input('--> ')
raw_input('--> ')
raw_input('--> ')