DataMedSci / pymchelper

Python toolkit for SHIELD-HIT12A and FLUKA
http://datamedsci.github.io/pymchelper/
15 stars 7 forks source link

values/pages from SH12A bdo file are not pr particle, bug? #400

Open kkorsgaard opened 3 years ago

kkorsgaard commented 3 years ago

Hi After I should compare some different versions of Shieldhit12A I did notice that different identical setup had different values for dose depending on how many simulated particles I did use.

The problem

I did find out that the reason for the different values are that pymchelper does output value from Shieldhit12A multiplied with how many particles I simulate pr simulation.

This means that even though I combine 200 simulations each with 100k primaries it only gives the total for ONE simulation. Very unintuitive if it is a feature and not a bug. If it should give the total I would like it to be the TOTAL for all simulations (primaries) I ran and independent of how I do arrive to this number.

Goal

But anyway I would rather have it to output value (e.g. dose) pr simulated particle as I thought it did. This is what Shieldhit12A do if I output as a txt file.

Reproduce the problem

I can easy reproduce the problem. Just by running a simulation and having both a bdo file and txt file saved. This is also the reason that I do say pymchelper multiply all values with number of simulated primaries as the txt and bdo files that way get identical (I multiply the txt values with number of primaries). Furthermore I am using the method fromfile(). Convertmc image gives the same plots/ values. I have attached a txt file and bdo file which are from the same run and setup. files.zip

Tested values

I have tested dose, Nkerma, Energy and materiel. It should at least not multiply material and zone at it does not depend on the number of simulated particles, so it indicates that it is a bug.

Version

I am using pymchelper version 1.5.1 and Shieldhit version: v0.9.1-6-gc0599d8.

kkorsgaard commented 3 years ago

As mentioned in post #395 when doing parallelization not all jobs will finish for a given time limit. So I guess which could lead to some potential errors - if pymchelper in 1/10 of the cases gives the total dose for 100k particles and in 9/10 gives the total doses for 200k which would lead to some large uncentainties. However I don't think it is something I have seen, but the question remains how does pymchelper handle this?

I cannot think of a very intuitive way (I would say it did something wrong if it e.g. multiplied all values with 200k even if some only have run with 100k primaries). So dose pr simulated particle would give more sense in these cases.

nbassler commented 3 years ago

@grzanka it might be helpful for pymchelper to add a normalization tag on SH side. SH_NORM_NONE SH_NORM_PER_PARTICLE ... ?

nbassler commented 3 years ago

Not directrly related, but SHBDO_PAG_NORMALIZE tag was not set previously on SH12A.

I propose to use it, then only SH12A needs to keep track of which detectors need to be normalized and which one do not. This is the MR which activates the tag: https://gitlab.au.dk/shieldhit/shieldhit/-/merge_requests/259

nbassler commented 3 years ago

@grzanka following our discussion, I realize we have 3 cases how to deal with results in .bdo.

so, we need these averaging cases:

1) SUM:          X = sum_j x_j                    # this is e.g for COUNT, COUNTPOINT and other artificial stuff.
2) NORMALIZED:   X = sum_j x_j / sum_j I_j        # this changes unit into "per primary". Dose, Fluence etc.
3) AVERAGED:     X = (sum_j x_j * sum_j I_j) / sum_j I_j      #  unit is unchanged. dLET, tLET, ...z2/beta2

I do not see any difference in fluence or dose averaged here, since the feld is identical for the voxel of the constiuent compoents, so fluence = kf * I_j and dose = kd * I_j

kf and kd are constants and independent of j (assuming sufficent statistics). Therefore they cancels out in case 3) if you do dose or fluence averaged, respectively.

SHBDO_PAG_NORMALIZE is a char, which can carry the modus.

nbassler commented 3 years ago

We can even add a 4th case, where just the first result is used.

0) X = x_0    # for GEOMAP scorers, they are always the same for all j.
nbassler commented 3 years ago

https://gitlab.au.dk/shieldhit/shieldhit/-/merge_requests/259

new tag behaviour (was not set at all before)

    SHBDO_PAG_NORMALIZE,             /* Flags for page->postproc on how to postprocess the data in SHBDO_PAG_DATA

                                        Given:
                                        - the data in page->data as x_j
                                        - for j instances of this simulation
                                        - which was done with I_j number of paritcles.

                                        The resulting data will be termed X and has the units given by SHBDO_PAG_DATA_UNIT

                                        0: X = x_1                                for GEOMAP type scorers
                                        1: X = sum_j x_j                          COUNT, ...
                                        2: X = (sum_j x_j) / (sum_j I_j)          NORMCOUNT, ...
                                        3: X = (sum_j x_j * I_j) / (sum_j I_j)    LET, ...

                                        (Defines are in score/sh_scoredef.h)
                                      */
grzanka commented 2 years ago

Part of the scorers (code 2) will be fixed in PR #521

nbassler commented 1 year ago

related: #628 which adds another averaging tag.

reviewpad[bot] commented 1 year ago

AI-Generated Summary: This issue describes a problem with the output values of pymchelper when processing Shieldhit12A simulations. The user reports that different identical setups produce varied dose values depending on the number of simulated particles. It seems that pymchelper multiplies the output values by the number of simulated particles, which leads to incorrect results. The goal is to output the value (e.g., dose) per simulated particle as Shieldhit12A does when outputting as a text file. The issue can be reproduced by running a simulation and having both a BDO and TXT file saved. The user has tested this for dose, Nkerma, Energy, and material parameters. The issue is observed in pymchelper version 1.5.1 and Shieldhit version v0.9.1-6-gc0599d8.

grzanka commented 5 months ago

Should be fixed in PR #723