Long-term, we need to adapt QRules so that it can handle an initial state with two particles (https://github.com/ComPWA/qrules/issues/29). For now though, you'll have to create a virtual particle that represents the state $p\gamma$. An example for $e^+e^-$ can be found here:
https://qrules.readthedocs.io/0.10.x/usage/particle/#adding-custom-particle-definitions-through-python

I'm not entirely sure what the quantum numbers of the virtual particle for $p\gamma$ would be. You can try out some combinations, like creating a Particle with isonspin=None. Be careful though that the combination of quantum numbers for the Particle complies with Gell-Mann–Nishijima.

For the record:

allowed reaction types: "strong", "EM", and "weak":

allowed reaction types: "strong", "EM":

allowed reaction types: "strong":

Awesome 👌