Predicting rolling contact fatigue crack hot spots or regions with increased local driving forces in rails is challenging due to the wide range of factors that influence crack initiation. Rail sections experience fluctuating creepage conditions, contact positions, and loads throughout their lifespan, influencing the development and location of fatigue cracks. A new computational method is proposed that predicts the orientation and regions prone to rolling contact fatigue cracks under realistic service loading. It combines multi-body simulations, finite element analysis, and critical plane approaches. A novel multi-variable sampling technique simplifies loading spectra into representative traction profiles, which are then analyzed using finite element analysis and the Smith–Watson–Topper damage indicator parameter (DIP<inf>SWT</inf>). The maximum DIP<inf>SWT</inf> value identifies the critical plane and potential crack orientation. A case study on the Swedish heavy haul train line (Malmbanan) considers measured traffic and loading conditions, analyzing the wheel load spectrum for a 384 m long section of a R = 450 m curve. Results show that the DIP<inf>SWT</inf> is highest for the locomotive with a loaded payload configuration, with a maximum value of 3.84 × 10⁻⁸ located at 38.59 mm from the lower gauge face corner. The DIP<inf>SWT</inf> critical plane aligns with experimental measurements of RCF cracks orientations near the gauge corner. This computational method, when combined with other predictive tools, can efficiently identify conditions that lead to RCF cracks and determine their possible locations and orientations in railway tracks.