[Gate-users] Unable to use phase space file from IAEA with Gate 6.2

[Stathis Kamperis]

Ok, I "fixed" it for me.

There are many problems involved here.

1) The document from IAEA descriping the format of the phase space
files is not consistent with the phase space files that they provide
(specifically with the ELEKTA_PRECISE_6mv_part1.IAEAphsp that I
checked).

The document implies that the record of a particle starts with its
X,Y,Z coordinates, which is not the case. The data in the Elekta phase
space file are of the form:

<particle type> 1 byte   (sign of particle designates z direction)
<energy>         4 bytes (sign of energy designates history)
<x coord>        4 bytes
<y coord>        4 bytes
<u cosine>       4 bytes
<v cosine>       4 bytes
<weight>          4 bytes
<extra longs>   8 bytes

At a total sum = 33 bytes (just like what the header says). You could
easily see for yourself, by looking at the data when you adopt the
above format, i.e.:

* The 1st field is always 1, 2 or 3 which looks like particle type.
* If you plot the 2nd field for particles with type=1, you get a nice
x-ray energy spectrum.
* If you look at fields 5 and 6 you realize that the values are
bounded to [-1, 1], so these must be the direction cosines.
* The 7th field is always 0.05, which is the statistical weight (just
like what the header says).

2) The ELEKTA_PRECISE_6mv_part1.IAEAheader has 1 for W (not to be
confused with weight), but W does not occupy a field in its own in the
binary phase space file. The reason why this doesn't matter is because
the code doesn't even check for W and acts as if W was 0.

http://www.opengatecollaboration.org/lxr/source/6.2.0/source/general/src/GateIAEARecord.cc#L156

So far, 1) and 2) don't affect the runtime behaviour of Gate (they
only confuse the reader), because the code and the data files are
consistent with each other (but inconsistent with the IAEA document).

3) The problem is that even though Z is constant (just like the header
file says), for some reason Gate assumes it isn't and tries to read it
from the binary stream, which messes things up big time. How do I
know? Even though:

* The header file says that Z is 0.
* Gate succesfully reads that value from the header (I checked).
http://www.opengatecollaboration.org/lxr/source/6.2.0/source/general/src/GateIAEAHeader.cc#L60
* In get_record_contents() 'iz' is 0 (I checked).
* When read_particle() is called 'iz' is 1 (I checked).
http://www.opengatecollaboration.org/lxr/source/6.2.0/source/general/src/GateIAEARecord.cc#L156

I hard-coded 'iz' for me to be 0 and suddenly Gate started
interpreting the binary stream correctly and the whole simulation
worked fine. I don't have time to fully debug it and supply a patch,
because I operate under very tight time constraints, but I hope that
this mail will help others.

Best regards,
Stathis

On Sat, Mar 8, 2014 at 10:34 PM, Stathis Kamperis <ekamperi at gmail.com> wrote:
> On Sat, Mar 8, 2014 at 10:13 PM, Stathis Kamperis <ekamperi at gmail.com> wrote:
>> I'm not sure why Gate can't parse the phase space file. I tried to
>> extract the information with Mathematica and it seems that the phase
>> space file is ok (not sure about the header though, especially where
>> the record contents are declared).
>>
>> For example I get only 1,2,3 (gammas, e-, e+) for particles:
>>
>> file =
>>   BinaryReadList[
>>    "~/Dropbox/foo",
>>    Join[
>>     Flatten@
>>      {"Byte", PadLeft[{}, 6, ToString at Real32], "Integer32", "Integer32"}],
>>    ByteOrdering -> -1];
>>
Join[] is redundant here. A left over from previous experiments.

>> h = HistogramList[file[[All, 1]]]
>> {{1, 2, 3, 4}, {1543049, 8306, 161}}
>>
>> And the ratios seem to be in alignment with what is declared at the header:
>>
>> First at Last@h/Total at Last@h // N
>> 0.994543
>>
>> 125806703/126491813 // N
>> 0.994584
>>
>> I've only checked 50 MB worth of data but I doubt that there's a wrong
>> particle number in the rest (will check though).
>>
>
> The whole phase space file is ok. I checked it with the following code snippet:
>
> [stathis at fedora]~/Downloads% cat read.c
>                                                      22:32:29
> #include <assert.h>
> #include <stdio.h>
>
> struct rec {
>         char particle;
>         char pad[32];
> } __attribute__((packed));
>
> int
> main(int argc, char *argv[])
> {
>         FILE *fp;
>         struct rec r;
>         int cnt, rv;
>
>         fp = fopen(argv[1], "r");
>         assert(fp);
>
>         cnt = 0;
>         do {
>                 ++cnt;
>
>                 rv = fread(&r, sizeof(struct rec), 1, fp);
>                 assert(rv > 0 || (rv == 0 && feof(fp)));
>
>                 if (r.particle != 1 && r.particle != 2 && r.particle != 3) {
>                         printf("Wrong particle = %d at record = %d\n",
>                             r.particle, cnt);
>                 }
>         } while(!feof(fp) && !ferror(fp));
>
>         fclose(fp);
>
>         return 0;
> }
>
> [stathis at fedora]~/Downloads% gcc -Wall -ansi -pedantic read.c
> [stathis at fedora]~/Downloads% ./a.out ELEKTA_PRECISE_6mv_part1.IAEAphsp
> [stathis at fedora]~/Downloads%
> [stathis at fedora]~/Downloads% echo $?
> 0
> [stathis at fedora]~/Downloads%
>
>> So... has anyone ever successfully used the IAEA phase space files
>> with Gate 6.2 ?
>>
>> Best regards,
>> Stathis
>>
>>
>>
>>
>> On Fri, Mar 7, 2014 at 11:40 PM, Stathis Kamperis <ekamperi at gmail.com> wrote:
>>> Relevant:
>>>
>>> http://comments.gmane.org/gmane.comp.science.opengate.user/3115
>>>
>>> Best regards,
>>> Stathis
>>>
>>> On Fri, Mar 7, 2014 at 11:26 PM, Stathis Kamperis <ekamperi at gmail.com> wrote:
>>>> Greetings everyone,
>>>>
>>>> I downloaded it from here:
>>>> https://www-nds.iaea.org/phsp/photon1/
>>>>
>>>> ELEKTA_PRECISE_6mv_part1.IAEAphsp and its associated header
>>>> ELEKTA_PRECISE_6mv_part1.IAEAheader
>>>>
>>>> When I try to use it from within Gate with:
>>>>
>>>> /gate/source/addSource mybeam phaseSpace
>>>> /gate/source/mybeam/addPhaseSpaceFile
>>>> /home/stathis/Downloads/ELEKTA_PRECISE_6mv_part1.IAEAphsp
>>>>
>>>> I get the following error:
>>>> /home/stathis/Downloads/gate_v6.2/source/physics/src/GateSourcePhaseSpace.cc
>>>> (l.252): Source phase space: particle not available in IAEA phase
>>>> space format.
>>>>
>>>> I checked the source code:
>>>> http://www.opengatecollaboration.org/lxr/source/6.2.0/source/physics/src/GateSourcePhaseSpace.cc#L231
>>>>
>>>> and it appears as if Gate can't recognize some particle type that the
>>>> header file is including. But,
>>>> https://www-nds.iaea.org/phsp/photon1/ELEKTA_PRECISE_6mv_part1.IAEAheader
>>>> has only photons, electrons, positrons that Gate should recognize.
>>>>
>>>> Any thoughts?
>>>>
>>>> Best regards,
>>>> Stathis