I am trying to calcuate CLs upper limit of the signal cross-section by fitting the signal to a gaussian and the background to a polynomial using the scikit-hep/hepstats package.
There are several sources of uncertainities to the the generated signal events. The background is obtained by fitting to data itself (similar to H -> gamma gamma analysis).
Can you please tell me how to include the effect of systematics while calculating the upper limits?
I think this is a very broad question that maybe boils down to the question of "how to incorporate systematic uncertainties into the likelihood.
hepstats uses the models and the likelihood to infer the limits, but they're all given by the former, therefore they have to be encoded in there.
How to encode them depends strongly on your use-case and the type of fits you're doing. For example, in templated, binned fits, BB or BB-lite is often used, morphing for some templates. Other possibilities are including constraints on parameters. But ultimately, it all depends on the model that you're building and what is in the likelihood.
Update I highly assume here that you're using zfit. Then, I would recommend you to use a free parameter for your systematic and use a constraint on that, an example can be found here.
If you have more specific, physics related questions, it could also be worth to contact the authors directly.