Reflections in small scale problematic? -> 1mm scale

[mahnen]

Dear Sebastian,

great work, installed out of the box no issues. I wrote the following scenery quickly, to test it on a 10mmx10mmx10mm block of supposed LYSO:

{
    "functions": [
        {
            "name": "unity",
            "values": [
                [
                    2e-7,
                    1
                ],
                [
                    0.0000012,
                    1
                ]
            ]
        },
        {
            "name": "refraction_lyso",
            "values": [
                [
                    2e-7,
                    1.81
                ],
                [
                    0.0000012,
                    1.81
                ]
            ]
        },
        {
            "name": "refraction_light_guide",
            "values": [
                [
                    2e-7,
                    1.49
                ],
                [
                    0.0000012,
                    1.49
                ]
            ]
        },
        {
            "name": "+infinity",
            "values": [
                [
                    2e-7,
                    9000
                ],
                [
                    0.0000012,
                    9000
                ]
            ]
        },
        {
            "name": "zero",
            "values": [
                [
                    2e-7,
                    0
                ],
                [
                    0.0000012,
                    0
                ]
            ]
        }
    ],
    "colors": [
        {
            "rgb": [
                22,
                9,
                255
            ],
            "name": "blue"
        },
        {
            "rgb": [
                255,
                255,
                0
            ],
            "name": "yellow"
        },
        {
            "rgb": [
                255,
                91,
                49
            ],
            "name": "red"
        },
        {
            "rgb": [
                16,
                255,
                0
            ],
            "name": "green"
        },
        {
            "rgb": [
                23,
                23,
                23
            ],
            "name": "grey"
        }
    ],
    "media": [
        {
            "name": "lyso",
            "refraction": "refraction_lyso",
            "absorbtion": "+infinity"
        },
        {
            "name": "light_guide",
            "refraction": "refraction_light_guide",
            "absorbtion": "+infinity"
        }
    ],
    "default_medium": "lyso",
    "surfaces": [
        {
            "name": "lyso_face_inside",
            "material": "Transparent",
            "specular_reflection": "unity",
            "diffuse_reflection": "unity",
            "color": "yellow"
        },
        {
            "name": "lyso_face_outside",
            "material": "Transparent",
            "specular_reflection": "unity",
            "diffuse_reflection": "unity",
            "color": "blue"
        },
        {
            "name": "esr_inside",
            "material": "Phong",
            "specular_reflection": "unity",
            "diffuse_reflection": "zero",
            "color": "red"
        },
        {
            "name": "esr_outside",
            "material": "Phong",
            "specular_reflection": "unity",
            "diffuse_reflection": "zero",
            "color": "green"
        },
        {
            "name": "perfect_absorber",
            "material": "Phong",
            "specular_reflection": "zero",
            "diffuse_reflection": "zero",
            "color": "grey"
        }
    ],
    "children": [
        {
            "type": "Mesh",
            "id": 1,
            "pos": [0, 0, 0],
            "rot": {"repr": "tait_bryan", "xyz": [0, 0, 0]},
            "surface": {
              "inner": {"medium": "lyso", "surface": "esr_inside"},
              "outer": {"medium": "light_guide", "surface": "esr_outside"}
            },
            "vertices": [
                [0, 0, 0],
                [0.01, 0, 0],
                [0, 0.01, 0],
                [0.01, 0.01, 0],
                [0, 0, 0.01],
                [0.01, 0, 0.01],
                [0, 0.01, 0.01],                
                [0.01, 0.01, 0.01]
            ],
            "faces": [
              [0,5,4],
              [0,1,5],
              [1,7,5],
              [1,3,7],
              [3,6,7],
              [3,2,6],
              [2,4,6],
              [2,0,4],
              [4,7,6],
              [4,5,7]
            ]
        },
        {
            "type": "Mesh",
            "id": 2,
            "pos": [0, 0, 0],
            "rot": {"repr": "tait_bryan", "xyz": [0, 0, 0]},
            "surface": {
              "inner": {"medium": "lyso", "surface": "lyso_face_inside"},
              "outer": {"medium": "light_guide", "surface": "lyso_face_outside"}
            },
            "vertices": [
                [0, 0, 0],
                [0.01, 0, 0],
                [0, 0.01, 0],
                [0.01, 0.01, 0]
            ],
            "faces": [
              [1,0,2],
              [1,2,3]
            ]
        },
        {
            "type": "Disc",
            "id": 3,
            "pos": [0.005, 0.005, -0.003],
            "rot": {"repr": "tait_bryan", "xyz": [0, 0, 0]},
            "radius": 0.05,
            "surface": {
                "inner": {"medium": "light_guide", "surface": "perfect_absorber"},
                "outer": {"medium": "light_guide", "surface": "perfect_absorber"}
            }
        }
    ]
}

If I render this, it seems to work (at least the scaled up version to 1m is fine:

-> The disk is supposed to absorb and measure the position of the photons exiting the mesh indicated in green (this is "esr_outside)".

However, when scanning the data and looking at it, it seems that there is no actual reflections going on on the sides. I attach an illustration, where the black disk is at z=0, i.e. directly at the LYSO touching it, but still there is photons interacting with the disk outside of the red indicated crystal.

I attach the script that produced the above plot.

import json
import numpy as np
import tempfile
import os
import merlict_c89_wrapper
import photon_source
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import matplotlib as mpl

with open("20200526_001_single_lyso.json", "r") as my_file:
    SINGLE_LYSO_SCENERY = json.load(my_file)

RANDOM_SEED = 1
MAX_INTERACTIONS = 100

# photon src at center of lyso cube
N_PHOTONS = 17000
CENTER_POINT_SRC = photon_source.isotrop_point_source(
    emission_position=[0.005, 0.005, 0.005],
    num_photons=N_PHOTONS  # close to saint gobain data sheet for 511keV LYSO
)

def run(
    scenery=SINGLE_LYSO_SCENERY,
    object_idx=3,  # the disk beneath
    photons=CENTER_POINT_SRC,
    random_seed=RANDOM_SEED,
    max_interactions=MAX_INTERACTIONS,
):
    with tempfile.TemporaryDirectory(prefix="pet_tutorial_") as tmp_dir:

        with open(os.path.join(tmp_dir, "scenery.json"), "wt") as f:
            f.write(json.dumps(scenery, indent=4))

        merlict_c89_wrapper.wrapper.json_scenery_to_merlict_scenery(
            json_in_path=os.path.join(tmp_dir, "scenery.json"),
            merlict_out_path=os.path.join(tmp_dir, "scenery.mli"),
            octree_out_path=os.path.join(tmp_dir, "scenery.octree.mli")
        )

        intersection_table = merlict_c89_wrapper.wrapper.propagate(
            scenery_path=os.path.join(tmp_dir, "scenery.mli"),
            octree_path=os.path.join(tmp_dir, "scenery.octree.mli"),
            object_idx=object_idx,
            photons=photons,
            random_seed=random_seed,
            max_interactions=max_interactions
        )

        return intersection_table

steps = 21
z_start = -1e-9
z_stop = 0.01
count = 0

n_bins = 200
x_start = -0.05
x_stop = np.fabs(x_start)
binsx=np.linspace(x_start,x_stop,n_bins+1)
binsy=binsx

x_crystal_max = 0.01
y_crystal_max = 0.01
z_crystal_max = 0.01

N_GAMMA_RAYS = 3
plot_z_max = 8000/10*N_GAMMA_RAYS

# loop over multiple z position
for z_position in np.linspace(z_start,z_stop,steps):
    np.random.seed(RANDOM_SEED)
    data = []
    
    for gamma_ray_number in range(N_GAMMA_RAYS):
        SINGLE_LYSO_SCENERY["children"][2]["pos"][2] = z_position
        x_emit = x_crystal_max/2
        y_emit = y_crystal_max/2
        z_emit = z_crystal_max/2
        #x_emit = np.random.rand()*x_crystal_max
        #y_emit = np.random.rand()*y_crystal_max
        #z_emit = np.random.rand()*z_crystal_max
        CENTER_POINT_SRC = photon_source.isotrop_point_source(emission_position=[x_emit,y_emit,z_emit],num_photons=N_PHOTONS)

        data.append(run(
            scenery=SINGLE_LYSO_SCENERY,
            object_idx=3,
            photons=CENTER_POINT_SRC,
            random_seed=RANDOM_SEED,
            max_interactions=MAX_INTERACTIONS))

    data = np.vstack(data)

    # plot some
    plt.clf()
    plt.hist2d(data[:,0],data[:,1],bins=[binsx,binsy],norm=mpl.colors.LogNorm(1,plot_z_max))
    #plt.hist2d(data[:,0],data[:,1],bins=[binsx,binsy])
    ax = plt.gca()
    ax.set(aspect='equal')
    ax.add_patch(
        patches.Rectangle(
            xy=(-x_crystal_max/2, -y_crystal_max/2),  # point of origin.
            width=x_crystal_max,
            height=y_crystal_max,
            linewidth=1,
            color='red',
            fill=False
        )
    )
    plt.title("crystal:10x10x10mm^3, guide_z:%02.01f mm"%(z_position*1000))
    plt.colorbar()
    plt.tight_layout()
    plt.savefig("%02d.png"%count)
    count+=1

[relleums]

Clearly this is not expected. So it seems some photons leak out of the boxy crystal. This calls for a unittest in merlict_c89 itself.

[relleums]

So there were two issues in the young merlict_c89

1. Indeed structures in the 1e-3m scale were partly ignored when setting up the octree. I reduced this to 1e-9m for now.

2. The surface-normals of triangles were not normalized which ruined the reflections within the crystal.

I implemented unit-tests in merlict_c89, fixed it, and updated the pet_tutorial.

When testing again, you must wipe out the old installation of pet_tutorial before installing again. To be sure remove the automatically created pet_tutorial/pet_tutorial/merlict_c89_wrapper/wrapper.c

Thank you very much for pointing this out.

[relleums]

I adjusted your z_stop to run from 0 to ```-1e-2```` so that the sensor-disc moves down away from the crystal.