The HEGRA system of 4 Imaging Atmospheric Cherenkov Telescopes (IACTs) has been used to determine the flux and the spectrum of cosmic ray protons over a limited energy range around 1.5 TeV. Although the IACT system is designed for the detection of gamma-rays with energies above 500 GeV, it has also a large detection area of $ \simeq 10^6 m^2 \cdot 3 msr$ for primary protons of energies above 1 TeV and the capability to reconstruct the primary proton energy with a reasonable accuracy $Δ$E/E of 50% near this threshold. Furthermore, the principle of stereoscopic detection of air showers permits the effective suppression of air showers induced by heavier primaries already on the trigger level, and in addition on the software level by analysis of the stereoscopic images. The combination of both capabilities permits a determination of the proton spectrum almost independently of the cosmic ray chemical composition. The accuracy of our estimate of the spectral index at 1.5 TeV is limited by systematic uncertainties and is comparable to the accuracy achieved with recent balloon and space borne experiments. In this paper we describe in detail the analysis tools, namely the detailed Monte Carlo simulation, the analysis procedure and the results. We determine the local (i.e. in the range of 1.5 to 3 TeV) differential spectral index to be $γ_p = 2.72\pm 0.02_{stat.}{\pm 0.15}_{syst.}$ and obtain an integral flux above 1.5 TeV of $F(> 1.5 TeV) =3.1\pm 0.6_{stat.}\pm 1.2_{syst.} 10^{-2}/s sr m^2$.

39 pages, 13 figures, 10 tables, prepared with REVtex. accepted by Physical Review D, in press