Setting accolinearity with back to back source

[MaxTousss]

Greetings,

If one looks into the mailing list or in the Gate code, one can find that there is an option to simulate accolinearity from back to back source. The commands to use are as follow:
/gate/source/NAME/setAccolinearityFlag true
/gate/source/NAME/setAccoValue 0.5 deg

In short, is this feature validated and still supported?

In the following, I detail what I found from my search on that feature.

This feature is not mentioned in the current openGate wiki (or at least, it was not found with the wiki search tool). I did a quick "cntrl+f" search of "setAccolinearityFlag" for the userguides available from v6.0 to the last available and I did not found a single mention of that keyword. However, if we look at the mailing list, we can find the following:
A possible bug was found in 2011 for which the resulting thread does not provide a clear conclusion about if the bug was corrected or if it was a real bug [1].
Around two years later, the feature was proposed to a user [4].
This feature seems to be still used currently [2], [3]. In [3], the user asked for some tips on the feature which were not answered.
From Gate v5.0 to v6.0, a bug was found and corrected [5], [6], [7].
I did try to look into the available release note of Gate and the Gate github but I did not find a mention about that feature.

Currently, the option to test Gate with accolinearity validation does not work out-of-the-box with that feature. When I tried to activate the accollinearity error graph creation, the resulting graph indicated pure back-to-back events. I have built some tests by hand[*] and it showed that a Gaussian angular error was indeed simulated which is good. However, when I compared the accolinearity error resulting from that feature to when using a positron source[**], the profile were clearly different. The latter being Gaussian-like, as expected. It would seems that something unexpected happens when the angular error is applied to the gamma pairs momentum with that feature. While I did try to make sense of the code in the method "GenerateVertex" of GateBacktoBack.cc, I must admit I fail to comprehend all the steps it does.

Thus, it seems that this feature might have some problems or that it is not completely supported. If it is indeed the case, I will propose a modification to the wiki such that it mention the existence of that feature but warns the users that it might not work as intended. If it should work, I will try to make a MWE of my tests and ask for validation.

Best,
Maxime Toussaint

[*] Outputting the "dev" variable from the file GateBacktoBack.cc
[**] Outputting the momentum from the file GateBacktoBack.cc and GatePositronAnnihilation.cc, and extracting the photon pairs accolinearity.

[dsarrut]

Hi, thank a lot! We will study and probably contact you in the future to resolve this issue.

[MaxTousss]

Greetings,

The following is a kind of MWE to observe the issue. At the end, you will obtain an image that should be simular to
which indicate that the acolinearity applied by the acolinearity flag for back to back source does something wrong. The figures show the histogram of the acolinearity amplitude (absolute value) in radian extracted from the emission direction of the gamma pair generated in a simulation. The left figure is obtained from a source of F18 positron while the right figure is obtained from a back-to-back source using the acolinearity flag with the value set at 0.5 Deg.

All the file needed are inside the following: mwe.tar.gz

To generate the data, you need to modify the Gate code. See modifToCode.txt. Note that this test was made with Gate v8.0. Do not forget to compile Gate after the modification.

Then, launch the simulations of F18 positron (acoF18.mac) and BTB source (acoBTB.mac). From there, you should only need to launch python3 resultExtractoy.py to obtain the previous image.

[dsarrut]

Doc has been updated. Accolinearity does not work in the current version

[dsarrut]

Mirjam is on it ;)

[m-dupont]

@MaxTousss If I redo the same as you in GateBacktoBack, ie compute angle between particle1 and particle, I find the right distribution.
See :

I think the problem in your python code, is that you do not divide |v1||v2| by the product of the norm of v1 and v2 before taking the arcosinus.

I used

x1 = dataAcoBtb[:, :3]
x2 = dataAcoBtb[:, 3:]
num = np.sum(x1*x2, axis=1)
denom = np.sqrt(np.sum(x1*x1, axis=1)* np.sum(x2*x2, axis=1))
plt.hist(np.arccos(num/denom) - np.pi, bins=50);

[MaxTousss]

@wrzof Yup, it seems I made some false assumptions.

For Positron sources, the vectors "Direction" and "DirectionPhoton" are both normalized to 1.0.
For BTB sources, the vector "gammaMom" is normalized to 0.511 while the norm of the vector "DirectionPhoton" as a distribution similar, in shape, to the image on the right shown in the post of Jan. 7.

Thus, my code, which assumed that the variables were normalized, gave the wrong conclusion and acolinearity seems to work with BTB sources!

However, is it possible that the difference in vector norm between the two type of sources might cause other problems?

[dsarrut]

Good catch, thanks a lot @wrzof ! I still keep that issue open for the meeting dedicated to PET e+ range (21/10/2021) for reference, but I will probably close the issue then.
@MaxTousss could you please add that information (about norms) in the documentation? Also, could you please send us the macro + python code so that we can prepare the benchmark? It will be very helpful I think.

[MaxTousss]

@dsarrut The macro + python code are already in the mwe.tar.gz file. Is there something missing or you simply want me to update the python code with wrzof correction?

For the documentation, you want me to add a warning about the back-to-back acolinearity norm behavior? If yes, I will do it sometime this week.

[mirjamlenz]

@dsarrut I am attaching updated files (for code modification and a python evaluation script), the GATE macros are still the ones from Maxime. The file for code modification (modifToCode.cc) includes the missing normalisation of gammaMom before the rotation DirectionPhoton.rotateUz(gammaMom);.
mwe_271021.zip

I double-checked the angle calculations from momentum vectors by making use of the G4ThreeVector function angle(), therefore the evaluation is done two times in the script. This results in four plots, which are created by the updated python script. Furthermore, it fits a Gaussian distribution to the histograms and outputs the fit parameters.
This script is helpful for testing the influence of the vector normalization, but for the benchmark we can narrow it down and only make use of the angle() function.

Some aspects we could add to the documentation:

• If the AccoValue is specified in degrees, it will automatically be converted to radians.

• The AccoValue corresponds to the FWHM of the distribution, not sigma.

• Specifying the acollinearity flag as true is not sufficient, since m_source->GetAccoValue() returns 0 when no value is specified in the macro. Can we add a default value somewhere?
  Currently it has to be specified:

/gate/source/addSource src
/gate/source/src/setType backtoback
/gate/source/src/setAccolinearityFlag True
/gate/source/src/setAccoValue 0.5887 deg

@MaxTousss
About the different vector norms: In the positron annihilation source code, the (normalised) momentum direction vector is being called in G4DynamicParticle(), which takes the momentum direction as an argument along with particle type and energy:

// GatePositronAnnihilation.cc

fParticleChange.AddSecondary( new G4DynamicParticle (G4Gamma::Gamma(), DirectionPhoton, E1) );
fParticleChange.AddSecondary( new G4DynamicParticle (G4Gamma::Gamma(), -Direction, E2) );

In the back-to-back source, the direction vector is parsed as the full momentum vector of a G4PrimaryParticle - so I assume it has to be 0.511 instead of 1 for this reason:

// GateBackToBack.cc

G4PrimaryParticle* particle = aEvent->GetPrimaryVertex( 0 )->GetPrimary( 0 );
G4PrimaryParticle* particle1 = aEvent->GetPrimaryVertex( 0 )->GetPrimary( 1 );

[…]

particle1->SetMomentum( DirectionPhoton.x(), DirectionPhoton.y(), DirectionPhoton.z() );
particle->SetMomentum( -gammaMom.x(), -gammaMom.y(), -gammaMom.z() );

[dsarrut]

Hi, should be ok now: see this example: https://github.com/OpenGATE/GateBenchmarks/tree/master/t19_acollinearity

[MaxTousss]

Currently, the example use the command "/gate/source/setDebugPositronAnnihilationFlag True" which seems to not be available in Gate9.2*. As such, it can't be tested with the currently available stable versions of Gate.

However, I have no objection to claim the issue closed since the benchmark exist.

---

• Gate9.2 is March 2022 vs May 2022 for the command

[dsarrut]

right, the debug flag has been added after, it will be in 9.2.1 ;)
I close the issue, thanks !

[MaxTousss]

Greetings,

I tried using the t19_acollinearity benchmark from GateBenchmarks/master (pulled today) with Gate/develop (pulled last week) and it produced the following figure:

As you can see, the current version of the benchmark seems to only show one of the four figures and the measured sigma is 0.22 compared to the previously shown 0.25.

Is this discrepancy normal? The problem might be due to the version (i.e. branch) of code I use.