[1337] | 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 | # ================================================================== |
---|
| 81 | class MyPrimaryGeneratorAction(G4VUserPrimaryGeneratorAction): |
---|
| 82 | "My Primary Generator Action" |
---|
| 83 | |
---|
| 84 | def __init__(self): |
---|
| 85 | G4VUserPrimaryGeneratorAction.__init__(self) |
---|
| 86 | self.particleGun= G4ParticleGun(1) |
---|
| 87 | |
---|
| 88 | def GeneratePrimaries(self, event): |
---|
| 89 | self.particleGun.GeneratePrimaryVertex(event) |
---|
| 90 | |
---|
| 91 | # ------------------------------------------------------------------ |
---|
| 92 | class MyRunAction(G4UserRunAction): |
---|
| 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 | # ------------------------------------------------------------------ |
---|
| 101 | class MyEventAction(G4UserEventAction): |
---|
| 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 | # ------------------------------------------------------------------ |
---|
| 112 | class MySteppingAction(G4UserSteppingAction): |
---|
| 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= demo_wp.MyPhysicsList() |
---|
| 170 | gRunManager.SetUserInitialization(myPL) |
---|
| 171 | |
---|
| 172 | # set user actions... |
---|
| 173 | myPGA= MyPrimaryGeneratorAction() |
---|
| 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 | |
---|