G4StepLimiterPhysics

[tobiascremer]

Hey it's me again. I'm still working with Geant4 in Julia and i have to say it's well done and very comfortable to work with.
I tried to implement a maximum stepsize into my simulation and noticed that it doesn't work although you use it in one of
your examples. Did i just oversee something or is the problem that the G4StepLimiterPhysics list isn't implemented yet?
Or could it be that the maximum stepsize just doesn't work for geantinos?

[peremato]

Hi. I didn't know that we need to set the User Limits and in addition to Register and instance of G4StepLimiterPhysics to the physics list. This class was not wrapped. I can add it and make a new release. Is anything else missing at this moment hat you noticed?

[tobiascremer]

Thanks for the fast answer a fast Implantation would bei great :). At least nothing i noticed ist Missing at the Moment but i will contact you again If i find Something.

[peremato]

I have made a new version (Geant4 v0.1.16) and is now registered in Julia. To use the new class G4StepLimiterPhysics you can define a custom physics list like this.

#----Physics---------------------------------------------------------------------------------------
struct B2PhysicsList <: G4VUserPhysicsList
  function B2PhysicsList(verbose)
      pl = FTFP_BERT(verbose)
      lp = G4StepLimiterPhysics()
      SetApplyToAll(lp, true)            # Apply to all particles (including neutrals)
      RegisterPhysics(pl, move!(lp))     # Register to the physics list
      return pl
  end 
end

The example B2a in https://github.com/JuliaHEP/G4Examples.jl/blob/main/basic/B2/B2a.jl has it now.

I would be interested if you could share your progress in your work since I have to prepare a contribution for JuliaCon and I could use it as example of use.

BTW, I did produce a set of tutorials notebooks here: https://github.com/peremato/Geant4.jl-tutorial

[tobiascremer]

Thanks for the quick response and implementation. I have tried using the maximum step size now for my code but it still doesn't work so I tried the B2a example with some modifications to look at the influence of the maximum step size. I added the following code after the inclusion of DetectorB2a.jl which basicly just opens a CSV-file to fill with data.

if isfile("test.csv")
	rm("test.csv")
end
const FILE_HANDLE = open("test.csv", "w")
write(FILE_HANDLE, "X,Y,Z,Edep\n")

I also added the following few lines into the processHits-function as a replacement for what was there before to write the exact position of every step and the energy deposition into my CSV-file.

  edep = step |> GetTotalEnergyDeposit
  pos = step |> GetPostStepPoint |> GetPosition
  position = [getX(pos), getY(pos), getZ(pos)]
  row = join(vcat(position, edep), ",")
  println(FILE_HANDLE, row)

(See fully modified code below)
So what i would expect to see now if i understood the maximum step size correctly is that (since we shoot in z-direction) the distance between two entries of my CSV-file in the z-coordinate is at max 2mm, but what I get as an output is the following.

X,Y,Z,Edep
-11.725465096139606,2.6111684519482248,-1659.66009865716,0.014426553714888745
-11.901514190642647,2.650268129196301,-1648.768496154882,0.004973045925526593
-12.428164063466637,2.7673874130911114,-1615.725968023137,0.017040183503037144
-12.947641904784359,2.8838450482893485,-1583.023461971727,0.011520887279207456
-13.242574556982547,2.947793670493517,-1564.5161921541794,0.00680453884670115

So in the 3rd column i would expect the difference between two lines to be less than 2mm but it shows that there are differences which are around 11mm instead. Did I just misunderstood the maximum step size or is still something not working for me? And if I misunderstood the step size how would I achieve the effect of limiting the distance traveled by a particle between each step?

Btw, i will ask my supervisors if I can provide some examples.

Full code of B2a.jl:

using Geant4
using Geant4.SystemOfUnits
using Printf, GeometryBasics

#---Define Detector Parameters struct--------------------------------------------------------------
include(joinpath(@__DIR__, "DetectorB2a.jl"))

if isfile("test.csv")
	rm("test.csv")
end
const FILE_HANDLE = open("test.csv", "w")
write(FILE_HANDLE, "X,Y,Z,Edep\n")

#---Define Tracker Hit-----------------------------------------------------------------------------
struct TrackerHit
  trackID::Int32
  chamberNb::Int32
  edep::Float64
  pos::Point3{Float64}
end
function Base.show(io::IO, hit::TrackerHit)
  (;trackID, chamberNb, edep, pos) = hit
  @printf(io, "\ntrackID: %3d chamberNb: %2d Edep: %.3f MeV Position: (%3f, %3f, %3f)", trackID, chamberNb, edep/MeV, pos...) 
end

#--------------------------------------------------------------------------------------------------
#---Define Sensitive Detector----------------------------------------------------------------------
#--------------------------------------------------------------------------------------------------
#---SD collected data------------------------------------------------------------------------------
struct B2aSDData <: G4JLSDData
  trackerHits::Vector{TrackerHit}
  B2aSDData() = new([])
end
#---Initialize method------------------------------------------------------------------------------
function _initialize(::G4HCofThisEvent, data::B2aSDData)::Nothing
  empty!(data.trackerHits)
  return
end
#---End of Event method----------------------------------------------------------------------------
function _endOfEvent(::G4HCofThisEvent, data::B2aSDData)::Nothing
  return
end
#---Process Hit method-----------------------------------------------------------------------------
function _processHits(step::G4Step, ::G4TouchableHistory, data::B2aSDData)::Bool
	#edep <  0. && return false
	#pos = step |> GetPostStepPoint |> GetPosition
	#push!(data.trackerHits, TrackerHit(step |> GetTrack |> GetTrackID,
	#                                   step |> GetPreStepPoint |> GetTouchable |> GetCopyNumber,
	#                                   edep,
	#                                   Point3{Float64}(x(pos),y(pos),z(pos))))
  edep = step |> GetTotalEnergyDeposit
  pos = step |> GetPostStepPoint |> GetPosition
  position = [getX(pos), getY(pos), getZ(pos)]
  row = join(vcat(position, edep), ",")
  println(FILE_HANDLE, row)
  return true
end
#---Create SD instance-----------------------------------------------------------------------------
chamber_SD = G4JLSensitiveDetector("Chamber_SD", B2aSDData();           # SD name an associated data are mandatory
                                    processhits_method=_processHits,    # process hist method (also mandatory)
                                    initialize_method=_initialize,      # intialize method
                                    endofevent_method=_endOfEvent)      # end of event method

#---End Event Action-------------------------------------------------------------------------------
function endeventaction(evt::G4Event, app::G4JLApplication)
  hits = getSDdata(app, "Chamber_SD").trackerHits
  eventID = evt |> GetEventID
  if eventID < 10 || eventID % 1000 == 0
    G4JL_println("Event: $eventID with $(length(hits)) hits stored in this event")
  end
  return
end

#--------------------------------------------------------------------------------------------------
#---Particle Gun initialization--------------------------------------------------------------------
#--------------------------------------------------------------------------------------------------
particlegun = G4JLGunGenerator(particle = "proton", 
                               energy = 3GeV, 
                               direction = G4ThreeVector(0,0,1), 
                               position = G4ThreeVector(0,0,-2940.0))

#----Physics---------------------------------------------------------------------------------------
struct B2PhysicsList <: G4VUserPhysicsList
  function B2PhysicsList(verbose)
      pl = FTFP_BERT(verbose)
      lp = G4StepLimiterPhysics()
      SetApplyToAll(lp, true)            # Apply to all particles
      RegisterPhysics(pl, move!(lp))     # Register to the physics list
      return pl
  end 
end
#--------------------------------------------------------------------------------------------------
#---Create the Application-------------------------------------------------------------------------
#--------------------------------------------------------------------------------------------------
app = G4JLApplication(;detector = B2aDetector(nChambers=5),          # detector with parameters
                       generator = particlegun,                      # primary particle generator
                       nthreads = 0,                                 # # of threads (0 = no MT)
                       physics_type = B2PhysicsList,                 # what physics list to instantiate
                       endeventaction_method = endeventaction,       # end event action
                       sdetectors = ["Chamber_LV+" => chamber_SD]    # mapping of LVs to SDs (+ means multiple LVs with same name)
                      )
              
#---Configure, Initialize and Run------------------------------------------------------------------                      
configure(app)
initialize(app)

#ui`/tracking/verbose 1`
beamOn(app,100)
close(FILE_HANDLE)

[peremato]

I did not check carefully the actual step length, but I noticed that by adding these lines in the script I saw that new steps are created because of the limits.

ui`/tracking/verbose 1`
beamOn(app, 1)

Can you try with geantino as particle type?

[tobiascremer]

Ok so I did some testing with different particles and different step size for one particle shot for each test.

particle_type, no_maximum_stepsize, 0.2mm, 20mm, 0.002mm

proton, 83 entries, 146 entries, 88 entries, 17 entries

gamma, 1603 entries, 1505 entries, 1395 entries, 1419 entries 

geantino, 6 entries, 6 entries, 6 entries, 6 entries

What I noticed from the testing is that the computation time goes up for smaller step sizes which makes sense but the number of entries generated seems totally random to me. I did wonder if my scorer doesn't score what I think it does, but when i checked the logic it should score every step with the according position and energy deposition. So technically i should see more entries for a smaller maximum step size.

[peremato]

Did you try with?

ui`/tracking/verbose 1`

It could be a discrepancy between the number of steps and the number of hits. In the processHit function there is

  edep = step |> GetTotalEnergyDeposit
  edep <  0. && return false

Perhaps is worth to put the question to the Geant4 Forum

[tobiascremer]

Oh my bad no i didn't. But with that I atleast found out my stupid mistake while using the example. I had my step size limit applied to the tracker logical volume while my sensitive detector was applied to the chamber logical volumes. Thats why I couldn't see the additional steps caused by the StepLimiter inside of my CSV-file. I also found the mistake in my own simulation. Sorry for the trouble and thanks for the great help :)

[peremato]

Welcome!