geant4koryak

import 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


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
        self.psum = Geant4.G4ThreeVector()
    def ProcessHits(self, step, rohist):
        energy = step.GetTotalEnergyDeposit()/MeV
        point = step.GetPreStepPoint().GetPosition()
        self.psum += energy * point
        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
        self.sd1.psum = Geant4.G4ThreeVector()
    def EndOfEventAction(self, event):
        self.hist1.Fill(self.sd1.eventEnergy)
        self.hist2.Fill(self.sd2.eventEnergy)
        #print "==== end of event #", self.count
        if self.sd1.eventEnergy > 0. :
                    sumP = (1. / self.sd1.eventEnergy) * self.sd1.psum
                    self.hist3.Fill(sumP.x/cm,sumP.y/cm)
                    #print "HIT!"
                    #print "X="+str(sumP.x)
                    #print "Y="+str(sumP.y)

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(500.0*cm,500.0*cm,500.0*cm)

Au = Geant4.G4Material.GetMaterial("G4_Au")
target1 = g4py.ezgeom.G4EzVolume("Plast")
target1.CreateBoxVolume(Au, 10*cm, 10*cm, 1.0*mm)
target1.PlaceIt(Geant4.G4ThreeVector(0.0*cm, 0.0*cm, 50*cm+0.5*mm))

Pb = Geant4.G4Material.GetMaterial("G4_Pb")
target2 = g4py.ezgeom.G4EzVolume("Cub")
target2.CreateBoxVolume(Pb, 100.0*cm, 100.0*cm, 100.0*cm)
target2.PlaceIt(Geant4.G4ThreeVector(0.0, 0.0, 0.0))



detector1 = MySD()
target1.SetSensitiveDetector(detector1)

detector2 = MySD()
target2.SetSensitiveDetector(detector2)


hist1 = ROOT.TH1D("hist1", "EnergyPlast", 1000, 0, 20)
hist2 = ROOT.TH1D("hist2", "EnergyCub", 1000, 0, 20)
hist3 = ROOT.TH2F("hist3", "average X Y", 100, -10. , 10., 100 , -10. , 10.)

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

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


pg = g4py.ParticleGun.Construct()
pg.SetParticleByName("e-")
#pg.SetParticlePosition( Geant4.G4ThreeVector(0.0*cm,50*cm+1*mm+100*cm,0.0*cm) )
pg.SetParticlePosition( Geant4.G4ThreeVector(0.0*cm,0.0*cm,50*cm+1*mm+100*cm) )

pg.SetParticleEnergy(10.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, 5.0*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, -1*vz) )
    Geant4.gRunManager.BeamOn(1)
    if(i<=10):
        raw_input('--> ')

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

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