Hi there,
Not quite sure I understand what your code does.
Differential photometry using Poissonian/Gaussian profile and then using gaia DR2 to get the magnitude of the reference star?
How does the fitted result differ from what you would get using DAOPHOT?
Also, do you set a value that allows stars that are nearing saturation to be ignored?
Is the star to be measured denoted by hand or by RA/Dec somehow?
It does look like code people would find useful.