How can I simulate a PHA data set in Sherpa, with/without added Poisson noise? 

The Sherpa fake_pha function is available for simulating a PHA data set, given a user-specified source model, integration time, and instrument response (ARF and RMF or RSP). The model parameters - which may have been previously fitted to an actual data set - are used to create the fake spectrum, which contains Poisson noise determined by the model parameters and integration time.

If you have already loaded in a PHA dataset along with a response and a source model then this can be done with a call like the following (for a 10 ks observation):

sherpa> fake_pha(1, arf=get_arf(), rmf=get_rmf(), exposure=10000)

New in CIAO 4.16 is the ability to change the method used to simulate the data from a model. The default is to use a Poisson routine (sherpa.utils.poisson_noise), but it can be set to a routine that just returns the model values, after converting to the nearest-integer value, rather than simulating noise:

sherpa> nonoise = lambda values, rng: np.rint(values)
sherpa> fake_pha(1, arf=get_arf(), rmf=get_rmf(), exposure=10000, method=nonoise)

For more information see the simulations section of the Sherpa threads.