Implementing a point detector estimator

[jmatta1]

Hi All,
As a start to implementing a somewhat more general “Finite Detector Estimator” I am trying to implement a Point Detector Estimator. After some thinking about it I am going to break this post into three pieces. In the first piece I want to describe my understanding of the estimator both as a refresher to myself and for anyone who knows more to correct any misunderstandings I might have. Second, I wanted to describe my approach to the implementation in hopes of constructive criticism. Finally, I will ask a couple of questions that I have come up with as I have been working on this problem.

If this you just want to see the questions, please see the third post.

Part 1: Point Detector Estimator

A Point Detector Estimator is a method for estimating the flux of some particle at a point in space and can be written as:

This equation is evaluated for every source particle and collision event throughout the random walk.

The first term is the “source term” where:
r_i is the position of the particle source
W_i is the particle weight when it is created at the source
P_i is the probability that the particle is emitted in some solid angle towards the detector r_i

The second term is the “interaction term” which is summing across every collision j in the history of particle i.
r_j is the position of the interaction
W_j is the particle weight prior to it’s interaction
P_j is the probability that some particle is emitted in some solid angle towards the detector from r_j

For both terms:
r_d is the position of the point detector
Beta(x, y) is the number of attenuation lengths (optical path length) between x and y for the particle which for gamma-rays can be expressed as:

Where:
s is each different geometry segment between x and y.
d_s is the length of that geometry segment
lambda_s is the attenuation length of the material in that geometry segment

So you can estimate the scalar flux at a point detector by summing some “Direct term” of the source with a “Scattering term”. The direct term (how much of the source shines directly through to the detector without collision) is calculated using beam attenuation to the detector from the source point chosen for that event multiplied by a solid angle correction and multiplied by a probability of the particle being emitted in that direction from that point (for non-isotropic sources). The scattering term is (how much of the source reaches the detector at some reduced energy due to scattering from surroundings) is calculated using beam attenuation to the detector from the source point chosen for that event multiplied by a solid angle correction and multiplied by a probability of some process emitting a particle towards the detector.

If I have made any errors so far, subtle or egregious, in this, please let me know.

[jmatta1]

Part 2: Implementation Approach

The next part I of this is my approach to implementing this in Geant4. I accomplish this using a combination of the various user actions.

A useful thing I realized is that the geantino particle is practically made for probing the geometry in a straight line (in fact it genuinely might have been made for that purpose). For the calculations of optical path length in the above equations I simple throw a geantino to give me the geometry segments.

G4VUserPrimaryGeneratorAction
I select a point in my source (a large cylindrical volume) and an isotropic direction from that point and I throw a gamma-ray of the appropriate energy from that point in that volume. Additionally, I generate a second primary vertex from the particle gun, a geantino traveling towards the detector.

G4UserSteppingAction
This is the part of my approach that generates almost all the questions. Here I check if the particle for the step is a gamma, a positron, an electron, or a geantino, otherwise the particle is ignored.

If the particle is a gamma and the step ended because of an interaction in the material (and the step was not with the detector region) then I either create or make a note to create two or three geantino secondaries of the appropriate energy (one for each process that could produce a secondary gamma travelling towards the detector Rayleigh Scattering, Compton Scattering, and -above 1.022 MeV- Pair Conversion followed by annihilation) with track information objects that contain the starting weight of each of those geantino. The starting weight is just product of the solid angle correction to the detector and, per secondary, the probability that the corresponding process occurred and produced a secondary pointing towards the detector.

If the particle in the step action is a geantino then the step length is divided by the total attenuation length for the material in the prestep point and that is added to the optical path length stored for that geantino. Additionally, if the post step point of the geantino is within the designated detector volume, the geantino’s track is set to fStopAndKill and the geantino’s weights are queued to be sent to my output thread.

If the particle is a positron or an electron and the step was ended due a discrete process then the probability of a bremmstrahlung photon being emitted in the direction of the detector is calculated and a geantino of appropriate energy and weight is generated.

If the particle is a positron and the step was ended due to annihilation then the probability of the annihilation photon being directed towards the detector is calculated and a geantino of appropriate energy and weight is created.

Any comments or suggestions on this approach are quite welcome.

[jmatta1]

Part 3: Questions

In the course of working out the approach I came to several questions.

Question 1:
I read somewhere that the only modification to the track that we are really allowed to do from UserSteppingAction is to change the track status. If that is the case: Am I not allowed to modify the track information object from within my UserSteppingAction (as I would want to add the next bit of optical path length)? Is there some better way?

Question 2:
What is the correct way to spawn these secondary geantinos? Can I just append to the G4TrackVector obtained from the Stepping Manager in my user stepping action?

Question 3:
Is there a good way to get at the angular distributions of the underlying physical processes? I can do some of what I want using the basic G4EmCalculator, but I am struggling to find a way to calculate the correct probabilities for the secondary going in the direction of the detector. Is there a good way to do this?

Question 4:
Is there some good way of setting the track information for geantinos created as primary particles in the event generator? Right now I have some harebrained scheme of having a little registry of each one created and then a tracking action that extracts the weights from that registry to set User information object at the beginning of tracking for that geantino.

Thank you so much for reading so far. If you can answer my questions or comment on anything I have written, I welcome it.

[jmatta1]

Since people don’t seem to have input, I figure I will use this as a place to record progress and if people want to drop in with wisdom at some point I would welcome it.

I just stumbled across this statement in the AppDeveloper manual:

Only processes can change information of G4Track and add secondary tracks via ParticleChange . G4VProcess is a base class of all processes and it has 3 kinds of DoIt and GetPhysicalInteraction methods in order to describe interactions generically. If a user want to modify information of G4Track , he (or she) SHOULD create a special process for the purpose and register the process to the particle.

Given the position in the manual I really should have seen it sooner and of course I see that it invalidates a good deal of my previous approach. My new approach involves essentially cloning the EM physics processes that I care about and see what it takes to modify those to generate my geantino ray-traces.

Still need to work out if I am allowed to tweak track information from within UserSteppingAction, since being able to do that instead of make some sort of replacement TransportationAction for geantinos.

[rajn]

Hi jmatta1
I am interested too in looking at point detector scorer as envisioned by you but for optical photons. I think the strategy will be as you describe and is probably very similar to point detector scorer in GAMOS (in case you aren’t aware) except that in there variance reduction methods are used from what I understood.
However, that particular method is also for Gamma and neutron and does not help me much.
I wonder if someone else have tried variance reduction method for optical photons to look at point detectors (by definition these detectors receive very little signal on account of their size in comparison to the rest of the geometry so that you need many photons to get a reduced variance in detection hits
Thanks.