Strange results scoring Phase space in transformed geometry

[jw3126]

When use the new egs_phsp_scoring on transformed geometries, I sometimes get strange results. Namely scored particles are not at the boundary of the scoring geometry. Here is an example:

:start geometry definition:
    :start geometry:
        name = world
        box size = 100.0 100.0 100.0
        :start media input:
            media = AIR521ICRU
            set medium = 0 0
        :stop media input:
        library = egs_box
        type = EGS_Box
    :stop geometry:
    :start geometry:
        name = box
        box size = 90.0 90.0 10.4
        :start media input:
            media = H2O521ICRU
            set medium = 0 0
        :stop media input:
        library = egs_box
        type = EGS_Box
    :stop geometry:
    :start geometry:
        name = tbox
        :start transformation:
            translation = 0.0 0.0 -0.2
        :stop transformation:
        my geometry = box
        library = egs_gtransformed
    :stop geometry:
    :start geometry:
        name = simgeo
        base geometry = world
        inscribed geometries = tbox
        library = egs_genvelope
    :stop geometry:
    simulation geometry = simgeo
:stop geometry definition:
:start run control:
    ncase = 10000
    nbatch = 10
:stop run control:
:start rng definition:
    initial seeds = 97 42
    type = ranmar
:stop rng definition:
:start source definition:
    :start source:
        name = source
        position = 0.0 0.0 10.0
        :start spectrum:
            energy = 1.0
            type = monoenergetic
        :stop spectrum:
        charge = 0
        library = egs_point_source
    :stop source:
    simulation source = source
:stop source definition:
:start ausgab object definition:
    :start ausgab object:
        name = hello
        phase space geometry = tbox
        output format = IAEA
        particle type = all
        score particles on = entry
        output directory = /tmp
        library = egs_phsp_scoring
    :stop ausgab object:
:stop ausgab object definition:

These are the positions of the first few particles of the phsp. I would expect that almost all particles enter from the top e.g. have z=5.0. However a lot of the particles are somewhere in the middle of the volume instead.

 [2.66136, -4.31234, 4.83054]  
 [-3.11729, -1.51948, 4.21642] 
 [-4.44435, -6.25954, 4.92993] 
 [-9.57176, -7.03649, 2.36184] 
 [5.18321, -10.0571, 2.56581]  
 [39.2378, 15.4158, 4.29871]   
 [19.0599, -2.30777, 1.36036]  
 [5.58726, -0.19545, -0.48804] 
 [-26.4626, 13.8884, -3.52076] 
 [12.9989, 7.46059, 0.126464]  
 [-21.8424, 17.4936, 2.96826]  
 [7.96049, -4.13504, 5.0]      
 [4.24452, -20.6891, 3.98471]

[blakewalters]

Hi @jw3126: The geom->isInside(x) function that phase space scoring using the "phase space geometry" method relies on is, unfortunately, subject to roundoff errors, and this is what's happening in your transformed box. For now, I would recommend using the "from regions", "to regions" method of scoring:

:start ausgab object:
    name = hello
    output format = IAEA
    from regions = 0
    to regions = 1
    particle type = all
    output directory = /tmp
    library = egs_phsp_scoring
:stop ausgab object:

Of course, this relies on knowledge of the region numbering scheme in the geometry. This is simple in your case above, but can be determined using egs_view in more complex geometries.

[jw3126]

Thanks for investigating this! Do I understand these correctly:

• The question if a particle hits a boundary is decided by howfar. This function is consistent with isInside only up to rounding errors. So when hitting a boundary, region numbers are updated correctly, but isInside could be true for both/none or one of the geometries.

• The issue is not transformed geometry specific, but can easily happen for most kinds of geometry?

[blakewalters]

Hi @jw3126: Regarding your first point: Yes, the distance to a region boundary is determined by howfar. isInside is an unrelated function and is only used to answer the question of whether the current position is inside the geometry from which it is called. It does not query the region number, which has been set based on howfar. As for your second point: Also, yes, this roundoff error could theoretically happen with any geometry, although this is the first time I've seen it. I think that the inverse transforms required to establish particle position/direction with respect to a transformed geometry may be introducing the errors in this case. Interestingly enough, when tbox is translated in the +Z direction, I don't see the issue. Still, the "phase space geometry" method was intended to work with any geometry (compound, transformed, etc) and so, until an accurate method of determining whether a particle is inside a geometry or not (using an as-yet-uncoded isInside method based on regions, perhaps), I would recommend using the "from/to regions" method just to be on the safe side.

[jw3126]

Thanks for answering these!