How to Find a Matching PSF Template

Suppose we want to model the PSF for the WFC observations of XMMJ2235 (the galaxy cluster at a redshift of 1.4) and the image in question is the F775W image taken at August 16, 2005 (DATASET ID=J9BA51E6Q). The file name of the flat-fielded image is j9ba51e6q_flt.fits. We need to py-drizzle the image to correct for geometric distortion (with the Lanczos3 kernel and an output pixel scale of 0.05" to match the scheme we used for the PSF library). The exposure time is 1285 seconds and the image contains numerous cosmic rays.
In principle, one can median-stack all the images for the cluster and remove cosmic rays on this file before we proceed any further. However, we still have ~10 high S/N stars that are not affected by cosmic rays. Therefore, we do not attempt to delete cosmic rays in this example. Many algorithms exist for detecting objects in astronomical images (e.g., daofind, imcat, SExtractor, etc.) and we believe that any of these performs well enough to find high S/N stars. Useful parameters for selecting such stars include flux error, FWHM, and ellipticity. We find 13 stars from the above image that satisfy our conditions for the "goodness". Their ellipticities are shown below.

The sticks are aligned to the direction of the elongation of the stars and each stick's length is proportional to its ellipticity. People whose expertise is not lensing may find the program quadrupole.pro helpful to compute the ellipticity. The table below summarizes the above stars' location, quadrupole moments, and their uncertainties. 
ID       X         Y        Q11     Q22      Q12   Q11_err Q22_err Q12_err
--------------------------------------------------------------------------
115   3365.63   3740.99   0.9312  0.8880  -0.0354  0.0264  0.0264  0.0167
155   2948.41   3590.92   1.0847  1.0050   0.0219  0.0044  0.0045  0.0029
219   1831.53   3351.89   0.9144  0.8712   0.0337  0.0037  0.0037  0.0023
469    809.06   2692.92   0.8837  0.9465   0.0210  0.0005  0.0005  0.0003
601   2322.94   2267.02   0.9248  0.9630   0.0774  0.0255  0.0252  0.0161
687    688.21   1912.16   0.8433  0.8695   0.0304  0.0341  0.0337  0.0213
695   1440.27   1865.78   0.8824  0.9423   0.0444  0.0011  0.0010  0.0007
868   3905.48   1235.64   0.9249  0.9928  -0.0060  0.0225  0.0221  0.0142
872    216.58   1219.80   0.8533  0.9334  -0.0289  0.0063  0.0061  0.0039
901   2130.83   1116.11   0.8377  0.9143   0.0006  0.0078  0.0076  0.0049
991   2652.06    806.25   0.9642  0.9464   0.0022  0.0018  0.0018  0.0011

Now with the table above we are ready to search our PSF library for the template that closely matches the pattern. As stated in the paper, the best matching PSF template minimizes the following chi^2: , where Qij and Q'ij are the measurement of the stars in the above image and the model values at the same location, respectively. We find that the stellar observation field for NGC104 taken on May 30, 2006 at 23:08:22 UT (shown below) resembles the pattern most closely among the used templates.

Of course, because the template field contains many stars, we can utilize them to model a PSF for any given object in the j9ba51e6q field.