Matching photon pathlength distribution to analytical form

[Melissa Wu]

Hi Qianqian,

I hope this message finds you well.

While using MCX recently I ran into a mismatch and was curious about your input. I am comparing the total pathlength distribution of photons traveling through a homogeneous box to the analytical form of P(s) in a semi-infinite medium described in Eq. 4.16 of David Boas' thesis:

​I ran a simulation in MCXLAB (using similar code as in demo_validation_homogeneous.m), re-simulated via replay to only extract the pathlengths of detected photons, and compared their distribution to the above analytical expression for the same optical configuration. The MCX paths are more right-shifted than the analytical solution:

​​

I tried making the detector size pretty small (0.1 mm), increasing the size of the cube medium, as well as changing the boundary reflection option, but none of them seemed to make that big of a difference. I've attached the code that I used to generate the above plot, but off the top of your head, would you happen to have any ideas of why I am observing this discrepancy?

Thank you so much in advance for any help you can provide. I really appreciate it.

Best,

Melissa

[Melissa Wu]

The original screenshots don't seem to be showing up so I've reattached them here and included them in the above Google drive link.

[Hirvi Pauliina]

Dear Dr. Fang and Melissa,

Thank you, Melissa, for sharing this analytical formula!

I believe that the primary reason for the mismatch between the analytical and MC distributions is that each photon trajectory recorded by MCX carries an associated probability, or weight, which is generally not equal to 1/(number of detected photons). Instead, one must compute the probabilities for each path length and utilize them to form the respective histogram. 

Adding the following to the end of your script:

% Compute the weight, or probability for each path length with:

% https://github.com/fangq/mcx/blob/master/utils/mcxdetweight.m .

cfg.unitinmm = 1; % assume 1 mm x 1 mm x 1 mm voxels

detw = mcxdetweight(detp2, cfg.prop, cfg.unitinmm); % array with detected photon weights

detw = detw./sum(detw); % normalize

% Form histogram for analytical model and extract bin information.

figure

h_analytical = histogram(pl_samples, 'Facealpha', 0.5, 'EdgeColor','none','Normalization','pdf');

h_NumBins = h_analytical.NumBins; 

h_BinEdges = h_analytical.BinEdges;

h_BinWidths = h_analytical.BinWidth;

% Form histogram for MC by summing normalized probabilities for each bin.

h_MCX = zeros(h_NumBins,1);

for b=1:1:h_NumBins

    b_rows = detp2.ppath>=h_BinEdges(b) & detp2.ppath<h_BinEdges(b+1); % find path lengths in current bin

    h_MCX(b,1) = sum(detw(b_rows,1)); % sum probabilities for paths

end

% Normalize with bin widths to match 'pdf' normalization setting above.

h_MCX = h_MCX./h_BinWidths;

hold on

histogram('BinEdges', h_BinEdges, 'BinCounts', h_MCX, 'Facealpha',0.5, 'EdgeColor','none')

xlabel('Photon pathlength (mm)')

ylabel('Probability density')

legend({'Analytical','MCX PH'})

I was able to produce a better match:

I believe that the remaining discrepancies can be explained with the limited number of detected photon packets in each bin.

If you search the mcx-users forum for "consider weights" , you might find additional useful discussions.

Hope this helps!

Best regards,

Pauliina Hirvi

---

From: mcx-...@googlegroups.com <mcx-...@googlegroups.com> on behalf of Melissa Wu <wu.mel...@gmail.com>
Sent: Tuesday, July 21, 2026 6:46 PM
To: mcx-users <mcx-...@googlegroups.com>
Subject: [mcx-users] Re: Matching photon pathlength distribution to analytical form

 

--
You received this message because you are subscribed to the Google Groups "mcx-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mcx-users+...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/mcx-users/b347a4fb-1e2c-4123-840c-431954c7f566n%40googlegroups.com.

[Qianqian Fang]

Pauliina's reply was spot-on! 

the pathlengths must be weighted by the photon's detected weight, as shown by Pauliina's sample code.

a similar quantity, TPSF at the detector, requires to perform similar accumulation/binning over time gates, and you can see a similar weighted summation is also performed (not in the length bins, but along the time-bins)

https://github.com/fangq/mcx/blob/master/utils/mcxdettpsf.m

Qianqian

To view this discussion visit https://groups.google.com/d/msgid/mcx-users/AS8P192MB132009BD670C455B66450AB8E9C22%40AS8P192MB1320.EURP192.PROD.OUTLOOK.COM.

[Melissa Wu]

Thanks so much, Qianqian and Pauliina! This was extremely helpful and insightful.