The spatial sensitivity of an ultrasound transducer, which strongly influences its suitability for different applications, depends on the shape of the transducer surface. Accurate simulation of these spatial effects is important for transducer characterization and design, and for system response modelling in imaging applications. In optoacoustic imaging, broadband transducers are used to capitalize on the rich frequency content of the signals, but their usage makes highly accurate simulations with general wave equation solvers prohibitively memory- and time-intensive. Therefore, specialized tools for simulating the isolated spatial focusing properties described by the spatial impulse response (SIR) have been developed. However, the challenging numerics of the SIR and the necessity to convolve the SIR with the wave shape generated by the optoacoustic absorber to simulate the system response lead to numerical inaccuracies of SIR-based methods. In addition, the approximation error of these methods cannot be controlled a priori. To circumvent the problems associated with the explicit calculation of SIR, we propose directly computing the convolution of the required wave shape with the SIR, which we call the spatial pulse response (SPR). We demonstrate that by utilizing an h-adaptive cubature algorithm, SPR can be computed with significantly higher accuracy than an SIR-based reference method, and the approximation error can be controlled with a tolerance parameter. In addition, the integration of accurate SPR simulations into model-based optoacoustic image reconstruction is shown to improve image contrast and reduce noise artifacts. Precise system characterization and simulation leads to improved imaging performance, ultimately increasing the value of optoacoustic imaging systems for clinical applications.