We have introduced a new Monte Carlo approach to calculate the titration isotherms in a canonical ensemble. Unlike conventional constant pH (cpH) simulation methods, which inherently operate in the grand canonical ensemble and assess the protonation state of molecules solely with respect to the pH within the reservoir, our canonical titration method works directly with the isolated system. The simulation method employs a reactive Monte Carlo algorithm to determine the protonation state of macromolecules in relation to the total number of protons present within the canonical simulation cell. To compute the pH of a fully equilibrated system, we have developed a new surface insertion Widom algorithm, which effectively circumvents the challenges associated with the bulk Widom particle insertion, particularly for extremely low hydronium ion concentrations. To accurately account for the long-range Coulomb forces, we have adopted the Ewald summation method, highlighting the significance of the Bethe potential in the precise calculations of pH of canonical systems. Although the present simulation method was developed within the framework of the primitive model, in which water is treated as a dielectric continuum, there is no conceptual difficulty of extending it to more realistic atomistic simulations. In this respect, the canonical approach is much easier to implement than the alternative grand canonical methods, which require a simultaneous insertion of an anion together with a protonation move, in order to preserve the overall charge neutrality . Clearly in a dense system with atomistic water, most of insertion attempts will be rejected. On the other hand, one can easily combine the present canonical approach with a molecular dynamics (MD) simulation -so that the evolution of the system is performed using standard MD algorithms with a suitable thermostat -combined with titration moves in which one of the hydronium ions is transformed into a water molecule, with a proton transferred to one of the polyelectrolyte sites. The acceptance probability for such a protonation move will be given by Eq. ( ), and similarly for a deprotonation move. The pH calculation introduced in this paper can also be combined with standard MD algorithms, so that the sampling needed to perform the average in Eq. ( ) can be performed using a MD simulation. The implementation of the present approach to atomistic system will be the subject of the future work.