OBJECTIVE Dose-averaged linear energy transfer (LETD) is one of the factors which determines relative biological effectiveness (RBE) for treatment planning in proton therapy. It is usually determined from Monte Carlo (MC) simulation. However, no standard simulation protocols were established for sampling of LETD. Simulation parameters like maximum step length and range cut will affect secondary electrons production and have an impact on the accuracy of dose distribution and LETD. We aim to show how different combinations of step length and range cut in GEANT4 will affect the result in sampling of LETD using different MC scoring methods. METHODS In this work, different step length and range cut value in a clinically relevant voxel geometry were used for comparison. Different LETD scoring methods were established and the concept of covariance between energy deposition per step and step length is used to explain the differences between them. RESULTS We recommend a maximum step length of 0.05 mm and a range cut of 0.01 mm in MC simulation as this yields the most consistent LETD value across different scoring methods. Different LETD scoring methods are also compared and variation up to 200% can be observed at the plateau of 80 MeV proton beam. Scoring Method one has one of the lowest percentage differences compared across all simulation parameters. CONCLUSIONS We have determined a set of maximum step length and range cut parameters to be used for LETD scoring in a 1 mm voxelized geometry. LETD scoring method should also be clearly defined and standardized to facilitate cross-institutional studies. CONCLUSIONS Establishing a standard simulation protocol for sampling LETD would reduce the discrepancy when comparing data across different centres, and this can improve the calculation for RBE.