Task4mine

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

n = int(sys.argv[1])
print n, ' iterations'
offset = 0.001
energy = 6.0
absorbedGamma = 0
#classes for registering of detector responce
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.hist1 = hist1
        self.sd2 = sd2
        self.hist2 = hist2
        self.hist3 = hist3
        self.count = 0
    def BeginOfEventAction(self, event):
        self.count += 1
        print '==== start of event #', self.count
        self.sd1.eventEnergy = 0
        self.sd2.eventEnergy = 0
    def EndOfEventAction (self, event):
        self.hist1.Fill(self.sd1.eventEnergy)
        self.hist2.Fill(self.sd2.eventEnergy)
        self.hist3.Fill(self.sd1.eventEnergy / (self.sd1.eventEnergy + self.sd2.eventEnergy))
        global absorbedGamma
        if (self.sd2.eventEnergy >= 0.99 * energy):
            absorbedGamma += 1

        print '==== end of event # ', self.count, ', energy1, energy2 = ', self.sd1.eventEnergy, self.sd2.eventEnergy


#draw
art = ROOT.TCanvas("art")
art.Clear()
art.Divide(2,2)
art.cd(1)

#initialization
g4py.EMSTDpl.Construct()
g4py.NISTmaterials.Construct()
g4py.ezgeom.Construct()
#create outer space
air = Geant4.G4Material.GetMaterial("G4_AIR")
g4py.ezgeom.SetWorldMaterial(air)
g4py.ezgeom.ResizeWorld(50.0 * cm, 50.0 * cm, 70.0 * cm)
#detector
Ge = Geant4.G4Material.GetMaterial("G4_Ge")
target1 = g4py.ezgeom.G4EzVolume("Ge Detector")
target1.CreateBoxVolume(Ge, 5.0 * cm, 5.0 * cm, 5.0 * cm)
target1.PlaceIt(Geant4.G4ThreeVector(0.0 * cm, 0.0 * cm, 7.5 * cm))
#absorber
Pb = Geant4.G4Material.GetMaterial("G4_Pb")
target2 = g4py.ezgeom.G4EzVolume("Pb Absorber")
target2.CreateBoxVolume(Pb, 20.0 * cm, 20.0 * cm, 20.0 * cm)
target2.PlaceIt(Geant4.G4ThreeVector(0.0 * cm, 0.0 * cm, 20.0 * cm))
#create sensitive detectors
detector = MySD()
target1.SetSensitiveDetector(detector)

absorber = MySD()
target2.SetSensitiveDetector(absorber)
#histagrammsi
hist1 = ROOT.TH1D("first hist", "Emitted energy in detector", 200, -0.1, 6.1)
hist2 = ROOT.TH1D("second hist", "Emitted energy in absorber", 200, -0.1, 6.1)
hist3 = ROOT.TH1D("third hist", "Ratio", 200, -0.1, 1.1)
#action
uaction = MyEventAction(detector, absorber, hist1, hist2, hist3)
Geant4.gRunManager.Initialize()
Geant4.gRunManager.SetUserAction(uaction)

Geant4.gApplyUICommand("/control/execute visinit.mac")
Geant4.gApplyUICommand("/vis/enable true")
#create particle gun
pg = g4py.ParticleGun.Construct()
pg.SetParticleByName("gamma")
pg.SetParticlePosition(Geant4.G4ThreeVector(0.0 * cm, 0.0 * cm, 0.0 * cm))
pg.SetParticleEnergy(energy * MeV)
#cycle shoots
for i in range(0,n):
    if (i < 10):
        raw_input()
    if (i == 10):
        Geant4.gApplyUICommand("/vis/enable false")
    phi = ROOT.gRandom.Uniform(0.0, 2.0 * pi)
    theta = ROOT.gRandom.Gaus(0.0, offset * pi / 180)
    vx = sin(theta) * cos(phi)
    vy = sin(theta) * sin(phi)
    vz = cos(theta)
    pg.SetParticleMomentumDirection(Geant4.G4ThreeVector(vx, vy, vz))
    Geant4.gRunManager.BeamOn(1)




print 'absorbed photons part: ', 1.0 * absorbedGamma / n

hist1.Draw()
art.cd(2)
raw_input()
hist2.Draw()
art.cd(3)
raw_input()
hist3.Draw()
raw_input()