hassan Let me just verify quickly if that is still a problem, it does seems to be that it changes quite a lot from one run to another (which had different qualities of the centroid)

Updated after fix introduced in PIPE2D-726. Residuals are now on the order of 3% which I presume is on the order of what I would expect. I will investigate if that is true.

I meant to comment that the centroiding errors would mess up the normalisation; sorry.

Are the 3% errors coming from the normalisation? We should be able to get the right at a level that doesn't affect sky subtraction; or is the 3% coming from the state of your PSF modelling? Or the photon noise?

Ok, the ``effect'' is due to two factors, first one being dominant.

1. The dominant effect: The pipeline (here: https://github.com/Subaru-PFS/drp_stella/blob/32139567407954633cfc6d0a9bac1a32776640db/python/pfs/drp/stella/subtractSky2d.py#L268) normalizes by:
a) computing sum(psf_model**2)
b) computing sum(psf_model*image)
c) computing renormalizing value as sum(psf_model**2)/sum(psf_model*image)
d) psf_model_renormalized=psf_model * value 

As we can see, this does not ensure that sum(psf_model_renormalized)==sum(image), as I expected (unless psf_model_renormalized = constant * image). 

When the difference between pfs_model and image is larger, due to bad centering or bad modeling, the ``effect'' increases. 

2. Less important: I have been reporting values from 14x14 array around the psf cut, while the pipeline normalizes in square given by size of input psf, i.e., 21x21 pixels. However, that is a subdominant effect to the first one.

The reason for this normalization is now understood. This is an optimal solution for the images dominated by noise errors (e.g., faint galaxies in HSC). In the future we might investigate if this is an optimal solution for this use case. At the moment closing this with Won't Fix and possibly reopening this discussion in the future.