I am not an XPS expert at all, but if I understand the Shirley background (please correct me if I am wrong), the idea is that the larger peak of a doublet (say, 2p splitting) has a ghost peak (from some energy loss effect) that overlaps with the second peak in the doublet, and so you want to subtract from the second, smaller peak a "background" that is proportional to the area (or height) of the larger peak.
You might be able to define a `shirley` background function that included that "ghost peak" (If I recall, there is also a smoother background, perhaps mostly instrumental that would be included), and then constrain the amplitude (and perhaps centroid too?) of that "ghost peak" based on the parameters of the main peak.
That might be missing some of the important parts of the Shirley function, but maybe that is enough to get you started?
If that physics is sort of correct, I might suggest thinking about a model for the main peak that included that "ghost peak", instead of thinking of the peak as a simple-ish function with a background that is proportional to the main peak. That is, shouldn't there be a function that includes the full response from that electronic level? Maybe that is just a different lineshape? For XRF analysis, I use a variation of the so-called hypermet function, which I think might be slightly similar in that it has a low energy tail (mostly to account for detector incompleteness).
I hope that helps... but yes, I would be happy to add any lineshape functions, or do we can do here to make XPS analysis easier or better.
--Matt