# 2017 Laboratory C: Monte Carlo Simulations: Error Analysis

### From VirtualPhotonics

**GOALS: Gain familiarity with the error characteristics of Monte Carlo simulations**

Open VTS_MATLAB folder and click on file short_course_monte_carlo_lab.m.

**I. Dependence of Fluence Predictions on Number of Photons Simulated**

Goal: This exercise explores how fluence estimates change with the number of photons simulated.

- Go to the section Example 1a. This example sets up a Monte Carlo simulation for a system with a collimated source impinging on a slab of tissue. The tissue optical properties are specified in the "tissueInput.LayerRegions" struct under "RegionOP".
tissueInput.LayerRegions = struct(... 'ZRange', ... {... [-Inf, 0], ... % air "z" range [0, 100], ... % tissue "z" range [100, +Inf] ... % air "z" range }, ... 'RegionOP', ... {... [0.0, 1e-10, 1.0, 1.0], ... % air optical properties [0.01, 1.0, 0.8, 1.4], ... % tissue OPs [ua, us', g, n] [0.0, 1e-10, 1.0, 1.0] ... % air optical properties } ... );

The vector with comment "tissue OPs" specifies μ

_{a}mm^{-1}, μ_{s}' mm^{-1}, g and n of the tissue slab. Fluence is calculated in a cylindrical coordinate system (ρ,*z*) for increasing number of photons launched from the source, N = 10, 100, 1000, 10000. - Configure an Analog simulation by setting "options.AbsorptionWeightingType" to 'Analog' (Line 17).
options.AbsorptionWeightingType = 'Analog';

- Set the tissue absorption to μ
_{a}=0.01 mm^{-1}(Line 46). - Execute Example 1a by right clicking on the cell so that it turns a pale yellow color and selecting "Evaluate Current Cell" or typing Ctrl+Enter.
- How do the fluence results change as N is increased?
- The relative error (RE) for the fluence is determined by dividing the standard deviation by the mean. For N=10,000, why does the relative error increase with distance from the source?
- Change the Monte Carlo estimator to Discrete Absorption Weighting (DAW) by setting "options.AbsorptionWeightingType" to 'Discrete'.
options.AbsorptionWeightingType = 'Discrete';

Execute the Example 1a cell with this new setting.

- How do the results using Analog and DAW compare for a given number of launched photons N?
- For N=10000 which estimator provides fluence results with the least RE?
- Now move down to Example 1b. This cell plots the difference in the relative error (RE) between Analog and DAW estimators.
- Execute this cell. In what spatial regions does DAW show smaller RE than Analog? Can you explain why?

**II. Spatially-Resolved Reflectance Predictions for Analog versus CAW Simulations**

Goal: This exercise compares error estimates of spatially-resolved reflectance using Analog versus Continuous Absorption Weighting (CAW) simulations.

- Example 2 sets up two Monte Carlo simulations using: 1) Analog, and 2) Continuous Absorption Weighting (CAW). The source is collimated and the tissue has optical properties [μ
_{a}, μ'_{s}, g, n] = [0.01, 1.0, 0.8, 1.4]. Spatially-resolved reflectance is calculated in ρ bins. The number of photons launched for each simulation is N=10,000. - Set the tissue absorption to μ
_{a}=0.01 mm^{-1}[0.01, 1.0, 0.8, 1.4], ... % tissue OPs [ua, us', g, n]

and execute this cell.

- What span of source-detector separations (ρ) does Analog have smaller standard deviations (SDs) than CAW? Where does CAW have smaller SDs? Can you explain why?
- Note the time each took to execute in the Matlab command window. Why does Analog run so much faster?
- Determine the Analog relative error in the last ρ bin (ρ=9.9 mm) by typing "d1SD(100)/d1.Mean(100)" in the matlab command window. Do the same for CAW "d2SD(100)/d2.Mean(100)". Which method, Analog or CAW, has the smaller relative error?
- The efficiency of an estimator is calculated as Efficiency =1/(R^2 T) where R is the relative error and T is the time it took for the simulation to execute. Which estimator has the larger efficiency for detector estimates at ρ=9.9 mm?
- Set the tissue absorption to μ
_{a}=0.1 mm^{-1}[0.1, 1.0, 0.8, 1.4], ... % tissue OPs [ua, us', g, n]

and execute the cell again.

- Looking at the results for μ
_{a}=0.01 mm^{-1}with μ_{a}=0.1 mm^{-1}, how do the two methods RE compare as absorption is increased? How does the efficiency compare?