We present a Monte Carlo approach for performing titration simulations in the canonical ensemble. The standard constant pH (cpH) simulation methods are intrinsically grand canonical, allowing us to study the protonation state of molecules only as a function of pH in the reservoir. Due to the Donnan potential between a system and an (implicit) reservoir of a semi-grand canonical simulation, the pH of the reservoir can be significantly different from that of an isolated system, for an identical protonation state. The new titration method avoids this difficulty by using the canonical reactive Monte Carlo algorithm to calculate the protonation state of macromolecules as a function of the total number of protons present inside the simulation cell. The pH of an equilibrated system is then calculated using a new surface insertion Widom algorithm, which bypasses the difficulties associated with the bulk Widom particle insertion for intermediate and high pH values. To properly treat the long range Coulomb force, we use the Ewald summation method, showing the importance of the Bethe potential for calculating the pH of canonical systems.