ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
demo.py
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file demo.py
1 #!/usr/bin/python
2 # ==================================================================
3 # python script for Geant4Py test
4 #
5 # gtest01
6 # - check basic control flow
7 # ==================================================================
8 from Geant4 import *
9 import demo_wp
10 import g4py.MedicalBeam
11 import ROOT
12 
13 # ==================================================================
14 # ROOT PART #
15 # ==================================================================
16 
17 # ------------------------------------------------------------------
18 def init_root():
19 # ------------------------------------------------------------------
20  ROOT.gROOT.Reset()
21 
22  # plot style
23  ROOT.gStyle.SetTextFont(42)
24  ROOT.gStyle.SetTitleFont(42, "X")
25  ROOT.gStyle.SetLabelFont(42, "X")
26  ROOT.gStyle.SetTitleFont(42, "Y")
27  ROOT.gStyle.SetLabelFont(42, "Y")
28 
29  global gCanvas
30  gCanvas= ROOT.TCanvas("water_phantom_plots",
31  "Water Phantom Demo Plots",
32  620, 30, 800, 800)
33 
34 # ------------------------------------------------------------------
35 def hini():
36 # ------------------------------------------------------------------
37  global gPad1
38  gPad1= ROOT.TPad("2D", "2D", 0.02, 0.5, 0.98, 1.)
39  gPad1.Draw()
40  gPad1.cd()
41 
42  ROOT.gStyle.SetPalette(1);
43 
44  global hist_dose2d
45  hist_dose2d= ROOT.TH2D("2D Dose", "Dose Distribution",
46  200, 0., 400.,
47  81, -81., 81.)
48  hist_dose2d.SetXTitle("Z (mm)")
49  hist_dose2d.SetYTitle("X (mm)")
50  hist_dose2d.SetStats(0)
51  hist_dose2d.Draw("colz")
52 
53  gCanvas.cd()
54  global gPad2
55  gPad2= ROOT.TPad("Z", "Z", 0.02, 0., 0.98, 0.5)
56  gPad2.Draw()
57  gPad2.cd()
58 
59  global hist_dosez
60  hist_dosez= ROOT.TH1D("Z Dose", "Depth Dose", 200, 0., 400.)
61  hist_dosez.SetXTitle("(mm)")
62  hist_dosez.SetYTitle("Accumulated Dose (MeV)")
63  hist_dosez.Draw()
64 
65 # ------------------------------------------------------------------
66 def hshow():
67 # ------------------------------------------------------------------
68  gPad1.cd()
69  hist_dose2d.Draw("colz")
70  gPad2.cd()
71  hist_dosez.Draw()
72 
73 
74 # ==================================================================
75 # Geant4 PART #
76 # ==================================================================
77 
78 # ==================================================================
79 # user actions in python
80 # ==================================================================
82  "My Primary Generator Action"
83 
84  def __init__(self):
85  G4VUserPrimaryGeneratorAction.__init__(self)
87 
88  def GeneratePrimaries(self, event):
89  self.particleGun.GeneratePrimaryVertex(event)
90 
91 # ------------------------------------------------------------------
93  "My Run Action"
94 
95  def EndOfRunAction(self, run):
96  print "*** End of Run"
97  print "- Run sammary : (id= %d, #events= %d)" \
98  % (run.GetRunID(), run.GetNumberOfEventToBeProcessed())
99 
100 # ------------------------------------------------------------------
102  "My Event Action"
103 
104  def EndOfEventAction(self, event):
105  gPad1.Modified()
106  gPad1.Update()
107  gPad2.Modified()
108  gPad2.Update()
109  ROOT.gSystem.ProcessEvents()
110 
111 # ------------------------------------------------------------------
113  "My Stepping Action"
114 
115  def UserSteppingAction(self, step):
116  pass
117 
118 # ------------------------------------------------------------------
119 class ScoreSD(G4VSensitiveDetector):
120  "SD for score voxels"
121 
122  def __init__(self):
123  G4VSensitiveDetector.__init__(self, "ScoreVoxel")
124 
125  def ProcessHits(self, step, rohist):
126  preStepPoint= step.GetPreStepPoint()
127  if(preStepPoint.GetCharge() == 0):
128  return
129 
130  track= step.GetTrack()
131  touchable= track.GetTouchable()
132  voxel_id= touchable.GetReplicaNumber()
133  dedx= step.GetTotalEnergyDeposit()
134  xz= posXZ(voxel_id)
135  hist_dose2d.Fill(xz[1], xz[0], dedx/MeV)
136  if( abs(xz[0]) <= 100 ):
137  hist_dosez.Fill(xz[1], dedx/MeV)
138 
139 # ------------------------------------------------------------------
140 def posXZ(copyN):
141  dd= 2.*mm
142  nx= 81
143 
144  iz= copyN/nx
145  ix= copyN-iz*nx-nx/2
146 
147  x0= ix*dd
148  z0= (iz+0.5)*dd
149  return (x0,z0)
150 
151 
152 # ==================================================================
153 # main
154 # ==================================================================
155 # init ROOT...
156 init_root()
157 hini()
158 
159 # configure application
160 #app= demo_wp.MyApplication()
161 #app.Configure()
162 
163 myMaterials= demo_wp.MyMaterials()
164 myMaterials.Construct()
165 
166 myDC= demo_wp.MyDetectorConstruction()
167 gRunManager.SetUserInitialization(myDC)
168 
169 myPL= FTFP_BERT()
170 gRunManager.SetUserInitialization(myPL)
171 
172 # set user actions...
174 gRunManager.SetUserAction(myPGA)
175 
176 myRA= MyRunAction()
177 gRunManager.SetUserAction(myRA)
178 
179 myEA= MyEventAction()
180 gRunManager.SetUserAction(myEA)
181 
182 #mySA= MySteppingAction()
183 #gRunManager.SetUserAction(mySA)
184 
185 # set particle gun
186 #pg= myPGA.particleGun
187 #pg.SetParticleByName("proton")
188 #pg.SetParticleEnergy(230.*MeV)
189 #pg.SetParticleMomentumDirection(G4ThreeVector(0., 0., 1.))
190 #pg.SetParticlePosition(G4ThreeVector(0.,0.,-50.)*cm)
191 
192 # medical beam
193 beam= g4py.MedicalBeam.Construct()
194 beam.particle= "proton"
195 beam.kineticE= 230.*MeV
196 #beam.particle= "gamma"
197 #beam.kineticE= 1.77*MeV
198 beam.sourcePosition= G4ThreeVector(0.,0.,-100.*cm)
199 beam.SSD= 100.*cm
200 beam.fieldXY= [5.*cm, 5.*cm]
201 
202 # initialize
203 gRunManager.Initialize()
204 
205 # set SD (A SD should be set after geometry construction)
206 scoreSD= ScoreSD()
207 myDC.SetSDtoScoreVoxel(scoreSD)
208 
209 # visualization
210 gApplyUICommand("/control/execute vis.mac")
211 
212 # beamOn
213 gRunManager.BeamOn(100)
214 
215 #ROOT.gSystem.Run()
216