Results for various values of λ
The Monte Carlo method is a numerical technique that using random numbers and probability to solve problems. It represents an attempt to model nature through direct simulation for any possible results, by substituting a range of values (a probability distribution) for any factor that has inherent uncertainty. The method is named after the city in the Monaco principality, because of roulette, a simple random number generator. The name and the systematic development of Monte Carlo method dates from about 1944.The name “Monte Carlo” refers to the Monte Carlo Casino in Monaco because of the similarity of statistical simulation to games of chance and was coined by Metropolis during the Manhattan Project of World War II, .
Monte Carlo is now used routinely in many fields, such as radiation transport in the Nuclear Engineering, Dosimetry in Medical Physics field, Risk Analysis, Economics… in all the applications the physical process of the solution is simulated directly based on the major components of a Monte Carlo algorithm that must be available during the simulations. The primary components of a Monte Carlo simulation are:
Probability density functions (pdf’s) the physical system must be described by a set of pdf’s;
Random number generator: a source of random numbers uniformly distributed on the unit interval must be available;
Sampling rule, a prescription for sampling from the specified pdf’s;
Scoring: the outcomes must be accumulated into overall tallies or scores for the quantities of interest;
Error estimation: an estimate of the statistical error (variance) as a function of the number of trials and other quantities must be determined;
Variance reduction techniques: methods for reducing the variance in the estimated solution to reduce the computational time for Monte Carlo simulation;
Parallelization and vectorization algorithms to allow.
In the field of nuclear engineering, deterministic and stochastic (Monte Carlo) methods are used to solve radiation transport problems. Deterministic methods solve the transport equation for the average particle behavior and also contain uncertainties associated with the discretization of the independent variables such as space, energy and angle of the transport equation and can admit solutions that exhibit non-physical features. Although the physics of photon and electron interactions in matter is well understood, in general it is impossible to develop an analytic expression to describe particle transport in a medium. This is because the electrons can create both photons (e.g., as bremsstrahlung) and secondary or knock-on electrons (δ-rays) and conversely, photons can produce both electrons and positrons. The Monte Carlo (MC) method obtains results by simulating individual particles and recording some aspects of their average behavior. The average behavior of particles in the physical system is then inferred from the average behavior of the simulated particles. This method also enables detailed, explicit geometric, energy, and angular representations and hence is considered the most accurate method presently available for solving complex radiation transport problems. For example the most important role of Monte Carlo in radiotherapy is to obtain the dosimetric parameters with high spatial resolution. As the cost of computing in the last decades continues to decrease, applications of Monte Carlo radiation transport techniques have proliferated dramatically. On the other hand, Monte Carlo techniques have become widely used because of the availability of powerful code such as BEAM, EGSnrc, PENELOPE and ETRAN/ITS/MCNP on personal computers. These codes able to accommodate complex 3-D geometries, inclusion of flexible physics models that provide coupled electron-photon and neutron-photon transport, and the availability of extensive continuous-energy cross section libraries derived from evaluated nuclear data files.
It should be noted that these codes are general purpose, and are therefore not optimized for any particular application and are strongly depended on the solution subject. One of the difficulties associated with Monte Carlo calculations is the amount of computer time required to generate sufficient precision in the simulations. Despite substantial advancements in computational hardware performance and widespread availability of parallel computers, the computer time required for analog MC is still considered exorbitant and prohibitive for the design and analysis of many relevant real-world nuclear applications especially for the problems with complex and large geometry. But there are many ways (other than increasing simulation time) in the Monte Carlo method that users can improve the precision of the calculations. These ways known as Variance Reduction techniques and are required enabling the Monte Carlo calculation of the quantities of interest with the desired statistical uncertainty. Without the use of variance reduction techniques in complex problems, Monte Carlo code should run the problem continuously for weeks and still not obtain statistically significant reliable results. The goal of Variance Reduction techniques is to produce more accurate and precise estimate of the expected value than could be obtained in analog calculation with the same computational efforts. Variance reduction parameters are vary with problem types so iterative steps must be repeated to determine VR parameters for different problems.
2. Conceptual role of the Monte Carlo simulation
The conceptual role of the Monte Carlo simulations is to create a model similar to the real system based on known probabilities of occurrence with random sampling of the PDFs.
This method is used to evaluate the average or expected behavior of a system by simulating a large number of events responsible for its behavior and observing the outcomes. Based on our experience concerning the distribution of events that occur in the system; almost any complex system can be modeled. Increasing the number of individual events (histories) improve the reported average behavior of the system.
In many applications of Monte Carlo the physical process is simulated directly and there is no need to even write down the differential equations that describe the behavior of the system. The only requirement is that the physical or mathematical system be described by probability density functions (PDF). Once the probability density functions are known the Monte Carlo simulation can proceed by random sampling from the probability density functions. Many simulations are then performed multiple trials or histories and the desired result is taken as an average over the number of observations. In many practical applications one can predict the statistical error the variance in this average result and hence an estimate of the number of Monte Carlo trials that are needed to achieve a given error
3. Accuracy, precision and relative error in Monte Carlo simulation
The first component of a Monte Carlo calculation is the numerical sampling of random variables with specified PDFs. Each random variable defines as a real number that is assigned to an event. It is random because the event is random and also is variable because the assignment of the value varies over the real values. In principle, a random number is simply a particular value taken on by a random variable.
When the random number generator is used on a computer, random number sequence is not totally random. Real random numbers are hard to obtain. A logarithm function made the random number and the function repeats itself over time. When the sequence walked through, it will start from the beginning. The typical production of random numbers is in the range between 0 and 1.
A sequence of real random numbers is unpredictable and therefore un-reproducible. A random physical process, for example radioactive decay, cosmic ray arrival times, nuclear interactions, and etc, can only generate these kinds of sequences. If such a physical process is used to generate the random numbers for a Monte Carlo calculation, there is no theoretical problem. The randomness of the sequence is therefore not totally random; this phenomenon is called pseudorandom.  Pseudo Random numbers look nearly random however when algorithm is not known and may be good enough for our purposes. Pesudo random numbers are generated according to a strict mathematical formula and therefore reproducible and not at all random in the mathematical sense but are supposed to be indistinguishable from a sequence generated truly randomly. That is, someone who does not know the formula is not supposed to be able to tell that a formula was used rather than a physical process. When using Monte Carlo simulation, it is desirable to have any variable depending on a uniform distributed variable, ρ. The probability,
The probability is in the range between 0 and 1, can be rewritten as a cumulative distribution:
The left side of equation (2) is the uniform distribution between 0 and 1 and
Monte Carlo results are obtained by simulating particle histories and assigning a score
By assuming a scalar value for each Monte Carlo simulation output, the Monte Carlo sample mean of the first
Where is the sample variance,
The standard deviation of scores has been obtained by the square root of the variance (, which is estimated via Monte Carlo method as
The standard deviation is obtained by the following equation:
To define the confidence interval in Monte Carlo estimation two statistical theorems are used: the law of large number and the central limit theorem.
The law of large number provides an estimate of the uncertainty in the estimate without any idea concerning the quantity of
To define confidence interval for the precision of a Monte Carlo result, the Central Limit Theorem of probability is used as follow:
And for large
λ is the number of standard deviation, from the mean, over which the unit normal is integrated to obtain the confidence coefficient. Results for various values of λ are shown in Table1. So to have confidence level the estimation for
For example for λ=1, the interval,  has a 68% chance of containing the true mean.
Eq. (8) shows that the deviation of the sample mean from the true mean approaches zero as n , and the quantity of present a measured of the deviation of the sample mean from the population mean by using
To construct a confidence interval for sample mean, that has a specified probability of the containing the true unknown mean, the sample standard deviation
The sample variance of is then given by:
It should be noted that the confidence intervals are valid only if the physical phase space is adequately sampled by the Monte Carlo calculation. The uncertainty of the Monte Carlo sampled physical phase space represents the precision of the simulation. There are several factors that can affect the precision such as tally type, variance reduction techniques and the number of histories simulated. Generally uncertainty or error caused by the statistical fluctuations of the
On the other hand the difference between the true mean and the sample mean is called the systematic error. To estimate the relative error at the 1δ level which represents the statistical precision Eq. (10) is used. 
In terms of Central Limit Theorem, the estimated relative error squared R2 should be proportional to . So as each history will take on average, the same amount of computer time and the used computer time, T, in a Monte Carlo calculation should be directly proportional to n (the number of histories); therefore R2T should be approximately constant. Thus, the metric of efficiency for a given tally, called the figure-of-merit (FOM), includes computer time as well and define as:
The FOM is also a tally reliability indicator in the sense that if the tally is well behaved, the FOM should be approximately constant (with the possible exception of statistical fluctuations early in the problem), and is thus an important and useful parameter to assess the quality (statistical behavior) of a tally bin. If the FOM is not approximately constant, the confidence intervals may not include the expected score value
Considering the following form of the previous relation can show the significant of the actual value of the FOM. 
The above expression shows a direct relationship between computer time and the value of the FOM. Increasing the FOM for a given tally will subsequently reduce the amount of computer time required to reach a desired level of precision. Thus the FOM can be used to measure the efficiency when the variance reduction techniques have been used. The ratio of FOMs before and after using the variance reduction techniques, gives the factor of improvement.
Another use of FOM is to investigate the improvement of the new version of a Monte Carlo code.
The ratio of the FOMs for identical sample problems gives the factor of improvement. When the FOM is not a constant as a function of
Another additional use of the FOM is to estimate the required computer time to reach a desired precision by considering: T~1/(R2 x FOM).
4. Variance reduction
The uncertainty of Monte Carlo simulation can be decreased by implementing some accurate physical models but this leads to longer calculation times. On the other hand, the accuracy of Monte Carlo dose calculation  is mainly restricted by the statistical noise, because the influence of Monte Carlo method approximations should be much smaller. This statistical noise can be decreased by a larger number of histories leading to longer calculation times as well. However, there are a variety of techniques to decrease the statistical fluctuations of Monte Carlo calculations without increasing the number of particle histories. These techniques are known as variance reduction. Variance Reduction techniques are often possible to substantially decrease the relative error,
Decreasing the standard deviation,
The main goal of all the variance reduction techniques is to increase precision and decrease the relative error. The precision of the calculation is increased by increasing the number of particle histories but needs a large amount of computer running time so to accelerate the Monte Carlo simulation and reduce the computing time these techniques are applied.[16-17]
Monte Carlo variance reduction techniques can be divided into four classes:
The truncation method like geometry truncation, time and energy cut off;
The population control method like Russian roulette, geometry splitting and weight windows;
The modified sampling method (source biasing); and
The partially deterministic method like point detectors.
5. Popular variance reduction techniques
Several of the more widely used variance reduction techniques are summarized as follow:
Geometric Splitting/Russian roulette is one of the oldest and most widely used variance reduction techniques, and when used properly, can significantly reduce the computational time of a Monte Carlo calculation.
Approximately 50% of CPU time is consumed to track secondary and higher-order photons and the electrons they set in motion. It is possible to remove a part of these photons by Russian roulette. 
Generally when particles move from a region of importance
In case of moving to a less important region Russian roulette is played and the particle is killed with probability
It means that the particles are killed to prevent wasting time on them.(Figure 3)
The objective of these techniques is to spend more time sampling important spatial cells and less time sampling unimportant spatial cells. This is done by dividing the problem geometry into cells and assigning each cell
Energy splitting/roulette are similar to geometric splitting/roulette except that energy splitting/roulette is performed on the energy domain rather that on the spatial domain.
Russian roulette can be shown that the weights of all particle tracks are the same in a cell no matter which geometrical path the tracks have taken to get to the cell, assuming that no other biasing techniques, e.g. implicit capture, are used. In the simulations if a track's energy drops through a prescribed energy level, the roulette game (based on the input value of the survival probability) is played. If the game is won, the track's history is continued, but its weight is increased by the reciprocal of the survival probability to conserve weight. Russian roulette is frequently, if not always, used in radiation transport problems and can be applied at any time during the life of a particle, usually after an interaction has taken place. Russian roulette always increases variance since it cuts off histories that could still contribute to the detector, but it also always reduces the simulation time in compare with an implicit capture (which will explained later) scheme without weight thresholds.
Generally, in a deep penetration shielding problem the number of particles diminishes to almost nothing in an analog simulation, but splitting helps keep the numbers built up. To have accurate and precise results it is recommended to keep the population of tracks traveling in the desired direction more or less constant that is, approximately equal to the number of particles started from the source. Particles are killed immediately upon entering a zero importance cell, acting as a geometry cutoff. Geometry splitting/Russian roulette works well only in problems without any extreme angular dependence. In the extreme case, if no particles ever enter an important cell where the particles can be split, the Splitting/Russian roulette is useless. Energy splitting and roulette typically are used together. Energy Splitting/roulette is independent of spatial cell. If the problem has space-energy dependence, the space-energy dependent weight window is normally a better choice.
Splitting and roulette are very common techniques in Monte Carlo simulation; not only because of their simplicity but also since they only deal with variance reduction via population control and do not modify pdfs, they can be used in addition to most other techniques to have more effect.
5.2. Energy cut off
A Monte Carlo simulation can be made much faster, by stopping a particle once its energy drops below certain threshold energy (cutoff energy). According to the particle energy and the material that the particle is travelling through, the travelling path length of the particle can estimate. If this path length is below the required spatial resolution, particles are terminated and assume their energy is absorbed locally. This can be done by energy cut off that terminate tracks and thus decrease the time per history.  Because low-energy particles can produce high energy particles, the energy cutoff can be used only when it is known that low-energy particles are either of zero or almost zero importance at the specific region (low energy particles have zero importance in some regions and high importance in others). In the Monte Carlo simulations Ecut is the photon energy cut-off parameter. It means, if a scattered photon is created with energy less than Ecut the photon will not be transported and the energy deposited locally. According to the above explanation seems the smaller Ecut the more accurate are the results but there are two criteria that should be considered in the simulations for selecting Ecut: (a) the mean free path (MFP) of photons with energy equal or less than Ecut should be small in compared with the voxel sizes or (b) the energy fraction carried by photons with energy less than Ecut is negligible compared with the energy fraction deposited. In terms of efficiency selecting the higher Ecut results in decreasing the CPU time, but on the other hand, selecting a higher value for Ecut can makes it a source of additional statistical fluctuations if it becomes comparable or even bigger than the average energy deposited by electrons. In this case the answer will be biased (low) if the energy cutoff is killing particles that might otherwise have contributed in the process even if N →∞. 
5.3. Time cutoff
A time cutoff is like a Russian roulette, with zero survival probability. The time cutoff terminates tracks and thus decreases the computer time per history. Particles are terminated when their time exceeds the time cutoff. The time cutoff can only be used in time-dependent problems where the last time bin will be earlier than the cutoff. The energy cutoff and time cutoffs are similar; but more caution must be considered with the energy cutoff because low energy particles can produce high-energy particles, whereas a late time particle cannot produce an early time particle. 
6. Weight window technique/weight window generator
The weight window technique administers the splitting and rouletting of particles based on space and energy dependent importance. This technique is one of the most used and effective variance reduction methods that deals with both the direct decrease of variance via a large number of samples (through splitting) and the decrease of simulation time via Russian roulette, and is therefore a very effective variance reduction technique. On the other hand this technique combines Russian roulette and splitting.
To apply this variance reduction technique a lower weight bound and the width of the weight window for each energy interval of each spatial cell should be considered.
If a particle's weight is below the lower weight bound, Russian roulette is performed, and the particle's weight is either increased to be within the weight window or the particle is terminated. On the other hand, if the particle's weight is above the upper weight bound, the particle is split such that the split particles all have weights within the weight window. If the particle's weight falls within the weight window, no adjustment is performed.
As shown in Figure 4, if a particle has a weight equal to
One problem that may arise when using weight windows is that over-splitting might occur when a particle enters a region or it is generate in a region with higher weight than the upper limit of the weight window in that region. This can usually be solved by modifying some of the weight window parameters. 
The weight windows generator is used to determine weight windows for the simulations. When generating weight windows, it is easy to generate unwanted zeros. Zero weight windows in a region are either due to particles not entering that region or due to particles that did enter the region but did not add to the tally score. To increase the number of particles that enter or generate in a region of the system, a uniformly distributed volumetric photon source that covers the whole system is used.
The weight window generator calculates the importance of each cell in the problem. This is done by noting that the importance of a particle at a point in phase space is equal to the expected score a unit weight particle would generate. Thus, the cell's importance can be estimated as follow:
As both the weight window and geometry splitting use the Russian roulette, there is a question concerning the difference between these two methods. The main differences are:
The weight windows are space–energy dependent whereas geometry splitting is only dependent on space;
The geometry splitting, splits the particles despite the weight of the particle but the weight window works completely in opposite way, it means before roulette is played and the particles split, the weight of the particle is checked against the weight window;
The geometry splitting is only applied at surfaces, but the weight window method is applied at surfaces and collision sites or both;
The geometry splitting method is based on the ration of importances across the surface but weight windows utilizes absolute weights;
As the geometry splitting is weight independent, will preserve any weight fluctuation but weight window can control weight fluctuation by other variance reduction techniques to force all particles in a cell to have a weight within a weight window.
The weight windows can be generated via the weight window generator but it requires considerable user understanding and intervention to work correctly and effectively.
7. Implicit capture
Implicit capture, survival biasing, and absorption by weight reduction are synonymous. Implicit capture is a very popular variance reduction technique in radiation transport MC simulations. The implicit capture technique involves launching and tracing packets of particles instead of one by one. At launch, each packet is assigned an initial weight
All of the variance reduction techniques vary the physical laws of radiation transport to sample more particles in regions of the phase-space that contribute to the objective.
To compensate for this departure from the physical laws of radiation transport, the concept of particle weight,
When implicit capture is used rather than sampling for absorption with probability
8. Forced collisions
The forced collision method is a variance reduction scheme that increases sampling of collisions in specified regions.
If the number of mean-free paths (MFP) to the next photon interaction be larger than the simulated phantom thickness, the photon will leave the region of interest without any interacting or depositing energy. Prediction of this event is not possible and to have precision results the photon behaviour must be traced through the interest region until it escapes from the region. In this case the computing time spent on the transport of escaping photons is then wasted and if the fraction of escaping photons is very large, the simulation will be very inefficient. The forced collision technique improves the efficiency by considering only the fraction of photons that interact in the phantom.
As shown in Figure 5, when a specified particle enters a region defined as the forced collision region, the incident particle splits into un-collided and collided particles. The un-collided particle passes through the current cell without collision and is stored in the bank until later when its track is continued at the cell boundary. The collided particle is forced to collide within the specified cell.
It means that the “passing-through case” and the “collision case” are analysed due to a single incident particle simulation. These split particles are weighted by the following Equation:
The probability of colliding within a distance
The particle’s weight is then reduced appropriately and Russian roulette is used to ensure that the calculation time is not increased significantly by the technique. 
9. Exponential transformation
The exponential transform also called path length stretching is a variance reduction technique designed to enhance efficiency for deep penetration problems (e.g. shielding calculations) or surface problems (e.g. build-up in photon beams). It is often used for neutron Monte Carlo simulation and is directly applicable to photons as well.
In applying the exponential transformation a stretching parameter is used to increase distances travelled in directions of interest, while in the use of Russian roulette and splitting other parameters are introduced to increase the death probabilities in regions of low importance and the number of independent particles found in regions of high importance. 
Exponential transformation samples the distance to collision from a non-analog probability density function. Specifically, it involves stretching the distance between collisions in the direction of interest and reducing the distance between collisions in directions of little interest by modifying the total macroscopic cross section as follow:
Where is the modified total cross section, is the true total cross section, p is the exponential transform parameter used to vary the degree of biasing, < 1, and is the cosine of the angle between the preferred direction and the particle’s direction.
It should be mentioned that the exponential transformation technique can produce large weight fluctuations and subsequently produce unreliable mean and variance estimates. Exponential transformation generally decreases the variance but increases the time per history.
Also it should be noted that due to the large weight fluctuations that can be produced by the exponential transform the exponential transform should be used accompanied by weight control.
In this chapter we discussed the accuracy, precision, relative error & figure of merit in the Monte Carlo simulation, the Monte Carlo method is considered to be the most accurate method presently available for solving radiation transport problems in nuclear engineering field but it is extremely expensive computationally, because this method should simulate individual particles and simulating the average behaviour of the particles in the system. Also for complex problems that the probability is small that a particle will contribute to the interest region, some form of variance reduction must be applied to reduce the required computer time to obtain the results with sufficient precision.
The disadvantage associated with Monte Carlo codes is that they require long calculation times to obtain well converged results, especially when dealing with complex systems. Computation time constitutes an important and a problematic parameter in Monte Carlo simulations, which is inversely proportional to the statistical errors so there comes the idea to use the variance reduction techniques.
In other word to shorten the calculation time and to decrease the error of the results obtained with Monte Carlo methods, i.e. to improve the efficiency of a Monte Carlo calculation, variance reduction techniques must be used so in this chapter we also discussed various variance reduction techniques. These techniques play an important role in reducing uncertainties and improving the statistical results. As already mentioned, it is possible to reduce the estimated variance of a sample for a given number of histories using techniques which increase the positive sampling efficiency and reduce the intrinsic dispersion of the contributions to the response, biasing the natural occurrence of more important paths and reducing the natural occurrence of less important, and especially ‘nearly zero’ importance paths.
In this chapter, several variance reduction strategies have been described with the aim of reducing CPU time.
In general assign the variance reduction techniques results in several factors that should be noted in the simulations:
Highly reliable systems
Financial engineering/quantitative finance
Lots of alternatives
Noisy gradient estimates
Real-time control using simulation
The usage of variance reduction methods in Monte Carlo simulations is not straightforward. There is a problem that these methods may be used as a black box leading to results not being analysed correctly. Note, all of the variance reduction techniques have a certain amount of over-head associated with their use, and in many situations the cost of this over-head outweighs the advantages of the technique.