I am trying to use ParticleHP in my model to create secondary neutrons from a reaction, but then go in and edit the energy and direction of the neutrons from either a stacking action or tracking action before they are transported.
All of the user actions I can find take a const G4Track* which does not allow me to change any of the properties of the track. Is there a simple way to edit the information of an existing track or should I attempt to implement this by defining my own physics process to generate the secondary?
Canonically, you should write your own process, which you set as StronglyForced (so it gets called no matter what) and put last in the PostStep indexing. In that process, you can get the energy and direction the track is “going to have” from G4Step::PostStepPoint(). Then you can populate your process’ ParticleChange with the new values you want. That allows the Geant4 tracking management system to handle updating the G4Track kinematics properly.
Would it be possible to “edit” an existing process’ ParticleChange method using a wrapper similar to what is done in cross-section biasing? I’m wondering if that might be simpler than writing a whole new process if I want to continue using the interaction length from an existing process.
You could certainly do a wrapper process, but the neutron physics involves a number of processes. Are you confident that you have just a specific process you want to modify?
Writing a forced process to just do ParticleChange is essentially no more complicated than writing a wrapper process. You don’t need to compute an interaction length (just return DBL_MAX and set the forceCondition). You don’t need a fancy IsApplicable() – just attach the process to the specific particle type you want. In PostStepDoIt(), check the G4Step’s process name and either do something or nothing.
One other question: I only want this process called once per particle. Does returning DBL_MAX for the interaction length combined with setting the condition to forced insure that the process is called once and then won’t be called a second time?
Just the opposite – the process will get called at every step (that’s the Forced or StronglyForced condition). In your PostStepDoIt, you would check the state of the G4Step to see whether you really want to do anything or not. If you decide not, then
It sounds like what you want to do is modify new tracks (are you changing the channel production or something?). Is that right? In that case, your process should catch the parent particle, and modify the secondaries that got attached to the G4Step by the parent’s interaction during that step. Again, you test the G4Step to see whether the preceding interaction is really what happened or not.
Yep exactly, I want to change new tracks generated during the simulation. Its my understanding that when a secondary is created, the track of the particle which created it is killed, so I wouldn’t be able to catch the parent particle after its creates the secondaries unless I write a new process for creating the secondaries myself.
The current approach I’m trying is as you said earlier, to write a process which catches the neutrons right after they’ve been created and then updates their energy and direction. However, I’m realizing I’m not sure I’ll be able to get the relevant information about the parent track, specifically energy and direction at the time of secondary particle creation.
Also, I’m new to writing physics processes, do you know if any of the examples show how to write a custom physics process?
This part isn’t true! If your new process runs after everything else, then the seconaries created elsewhere will already be stored in the G4Step. They aren’t fully-fledged tracks yet, and the parent still exists (it’s the G4Track passed into your process). So you can access those new secondaries from the G4Step and mess with them.
That may be in some of the extended examples, and is certainly in some of the advanced examples. Let’s see…
$ cd <my path to G4>/examples/
$ find . -name '*.hh' | xargs grep ': public G4VDiscreteProcess'
./advanced/hadrontherapy/include/HadrontherapyStepMax.hh:class HadrontherapyStepMax : public G4VDiscreteProcess
[... lots of others, mostly all StepMax ...]
Okay I did not realize that! I think all I need to do then is write a process which grabs the parent particle, checks for the correct secondaries, and then edits the tracks of those secondaries. In this case is it correct to initialize a ParticleChange on the tracks of the secondaries as opposed to the track of the parent?
And thank you for the pointer to the advanced examples I’ve found a few that seem to do similar PostStep process and they’ve been very helpful!
You do need to initialize a ParticleChange, but you don’t need to populate it with anything. The secondaries will already have been created and attached to the step by the preceding process’ own ParticleChange.
Since the secondaries have not yet been “tracked” (i.e., they have not been processed by the G4 tracking manager to become new, numbered tracks), you can edit them by doing a const_cast<G4Track*>. As you know, this is generally considered evil, but whatever…
Okay I see - what do you think about the idea of grabbing the existing secondaries, set those tracks to be killed, and then repopulating the particle change with new secondaries? Would that cause problems for the tracking manager?
Oh, now that’s clean! There’s a small CPU and memory fragmentation penalty involved with the new/delete cycles, but your process is definitely better than const_cast editing. It won’t cause any problems with the tracking manager, because from its perspectivem the secondaries don’t even exist yet.
You could test this for yourself, if you already have an SD or SteppingAction where you can access the G4Step. Apply the same criteria your process would use to pick the “right” G4Step. Go through the secondaries and print out GetTrackID() for each one. You should see all zeros, because the stacking manager hasn’t assigned them IDs yet, or put them onto the stack.
I have run into a small issue implementing this, because the GetSecondary() method returns a const G4Track for which I cannot set the status. I could use const_cast to allow me to set the status but then we’re sort of back were we started, unless there’s another way to clear the secondaries associated with a step.
Update: it appears I can just cast the output of the GetSecondary() method as a G4Track* instead of a const G4Track*. I believe I’m still committing the sin of const_cast-ing the track but hopefully just changing the track status won’t cause too many problems?