Best-fit template not agreeing with the best-fit photometry

When fitting some simulated HST+JWST fluxes, comparing the original EAZY to eazy-py, where I've fit the same data set with the same templates (and template error function), without using a prior. The minimum chi-square values (and the chi-square surfaces) agrees quite well between EAZY and eazy-py, and the resulting redshifts are similar (well, the redshifts corresponding to the minimum chi-square).

However, when using self.show_fit(), I've discovered that many of the targets have fits that do not seem to agree with the best-fit templates. For instance:

The F150W and F200W template fluxes are significantly higher than the best-fit template. I've gone and pulled this template out by examining how self.show_fit() calculates the template from the coefficients, and passed it through the F150W and F200W filters by hand:

When I fit this object with the original EAZY (and the same parameters, including the same templates), I get this fit (in green):

Notice that the best-fit photometry is almost identical between this fit and the eazy-py fit, but the template combination here leads to stronger line emission accounting for the high F150W and F200W fluxes. I don't know enough about how the coefficients are calculated within eazy-py to understand why this discrepancy occurs. The templates that I am using are the standard set of EAZY templates along with some generated from fsps:

# Template definition file
#
# No blank lines allowed (for now).
#
# Column definitions:
#   1. Template number
#   2. Template file name
#   3. Lambda_conv (multiplicative factor to correct wavelength units)
#   4. Age of template model in Gyr (0 means template is always used)
#   5. Template error amplitude (for INDIVIDUAL template fits)
#   6. Comma/space separated list of template numbers to be combined
#      with current template if combined fits are enabled.
#
# Sample entry:
# 1 [path_to_file]/template1.sed 1.0 14.7 0.2 2,3,5
1   templates/eazy_v1.1_sed7.sed    1.0 0 1.0   
2   templates/eazy_v1.1_sed6.sed    1.0 0 1.0   
3   templates/eazy_v1.1_sed5.sed    1.0 0 1.0   
4   templates/eazy_v1.1_sed4.sed    1.0 0 1.0   
5   templates/eazy_v1.1_sed3.sed    1.0 0 1.0   
6   templates/eazy_v1.1_sed2.sed    1.0 0 1.0   
7   templates/eazy_v1.1_sed1.sed    1.0 0 1.0   
8   templates/ssp_25Myr_z008_withem.sed    1.0 0 1.0   
9   templates/ssp_5Myr_z008_withem.sed    1.0 0 1.0   
10  templates/c09_del_8.6_z_0.019_chab_age09.40_av2.0.sed  1.0 0 1.0   
11 templates/erb2010_highEW.sed   1.0 0 1.0    
12  templates/tau_0.01_age_0.1_dust2_0.0_fsps_model.sed    1.0 0 1.0 
13  templates/tau_0.0398_age_0.3162_dust2_0.0_fsps_model.sed    1.0 0 1.0 
14  templates/tau_1.0_age_11.749_dust2_0.0_fsps_model.sed    1.0 0 1.0 
15  templates/tau_10.0_age_0.01_dust2_0.0_fsps_model.sed    1.0 0 1.0 
16  templates/tau_10.0_age_0.01_dust2_0.6_fsps_model.sed    1.0 0 1.0

And the simulated data that I'm using is:

# id z_spec f_HST_F435W e_HST_F435W f_HST_F606W e_HST_F606W f_HST_F775W e_HST_F775W f_HST_F814W e_HST_F814W f_HST_F850LP e_HST_F850LP f_NRC_F090W e_NRC_F090W f_NRC_F115W e_NRC_F115W f_NRC_F150W e_NRC_F150W f_NRC_F200W e_NRC_F200W f_NRC_F277W e_NRC_F277W f_NRC_F335M e_NRC_F335M f_NRC_F356W e_NRC_F356W f_NRC_F410M e_NRC_F410M f_NRC_F444W e_NRC_F444W  
# id z_spec F233 E233 F236 E236 F238 E238 F239 E239 F240 E240 F363 E363 F364 E364 F365 E365 F366 E366 F375 E375 F381 E381 F376 E376 F383 E383 F377 E377
4116 -9999.0 2.8299999237060547 7.302999973297119 3.0799999237060547 8.663000106811523 3.2699999809265137 6.857999801635742 3.3399999141693115 4.50600004196167 3.569999933242798 1.1139999628067017 5.039000034332275 1.0789999961853027 6.1539998054504395 0.7509999871253967 7.366000175476074 0.8029999732971191 7.269999980926514 0.7829999923706055 4.236999988555908 0.5210000276565552 3.7190001010894775 0.8180000185966492 4.986999988555908 0.5659999847412109 4.456999778747559 0.8560000061988831 4.303999900817871 0.7170000076293945  
4152 -9999.0 19.8799991607666 11.027000427246094 16.579999923706055 5.099999904632568 15.0600004196167 1.437000036239624 14.9399995803833 7.045000076293945 14.529999732971191 1.8279999494552612 17.166000366210938 1.0880000591278076 18.32900047302246 0.7570000290870667 26.575000762939453 0.8180000185966492 28.052000045776367 0.8059999942779541 28.492000579833984 0.6110000014305115 24.31399917602539 0.9589999914169312 24.016000747680664 0.652999997138977 24.506000518798828 0.9789999723434448 22.983999252319336 0.8100000023841858  
4157 -9999.0 4.400000095367432 9.006999969482422 3.9000000953674316 3.565000057220459 4.070000171661377 1.6019999980926514 4.360000133514404 9.281999588012695 5.329999923706055 2.114000082015991 3.8940000534057617 1.0770000219345093 7.982999801635742 0.7509999871253967 7.0980000495910645 0.796999990940094 7.150000095367432 0.7799999713897705 4.627999782562256 0.5460000038146973 3.75600004196167 0.8529999852180481 4.074999809265137 0.5860000252723694 4.761000156402588 0.8930000066757202 4.1579999923706055 0.7450000047683716  
4212 -9999.0 13.289999961853027 13.012999534606934 7.869999885559082 7.86299991607666 10.039999961853027 9.72700023651123 2.759999990463257 2.447999954223633 0.07999999821186066 8.678000450134277 1.6820000410079956 1.055999994277954 4.639999866485596 0.734000027179718 4.000999927520752 0.7799999713897705 3.6710000038146973 0.7599999904632568 2.308000087738037 0.5400000214576721 2.8239998817443848 0.8550000190734863 2.609999895095825 0.5870000123977661 2.877000093460083 0.8949999809265137 3.8559999465942383 0.7590000033378601

with fluxes in nJy. I can provide other files, including the param file, the extra templates file, and/or the h5 file, if necessary.

As an update, this also happens if you run the fitting on Object 4116 above using these templates only:

# Template definition file
#
# No blank lines allowed (for now).
#
# Column definitions:
#   1. Template number
#   2. Template file name
#   3. Lambda_conv (multiplicative factor to correct wavelength units)
#   4. Age of template model in Gyr (0 means template is always used)
#   5. Template error amplitude (for INDIVIDUAL template fits)
#   6. Comma/space separated list of template numbers to be combined
#      with current template if combined fits are enabled.
#
# Sample entry:
# 1 [path_to_file]/template1.sed 1.0 14.7 0.2 2,3,5
1   templates/eazy_v1.1_sed7.sed    1.0 0 1.0   
2   templates/eazy_v1.1_sed6.sed    1.0 0 1.0   
3   templates/eazy_v1.1_sed5.sed    1.0 0 1.0   
4   templates/eazy_v1.1_sed4.sed    1.0 0 1.0   
5   templates/eazy_v1.1_sed3.sed    1.0 0 1.0   
6   templates/eazy_v1.1_sed2.sed    1.0 0 1.0   
7   templates/eazy_v1.1_sed1.sed    1.0 0 1.0   
8   templates/ssp_25Myr_z008_withem.sed    1.0 0 1.0   
9   templates/ssp_5Myr_z008_withem.sed    1.0 0 1.0   
10  templates/c09_del_8.6_z_0.019_chab_age09.40_av2.0.sed  1.0 0 1.0   
11  templates/erb2010_highEW.sed   1.0 0 1.0

Does it have anything to do with combining too many templates at once?

I doubt anyone is looking at this, but I think that I've discovered the issue - it has to do with the order of the template files specific, and to understand we have to look at the show_fit function:

        for i in range(self.NTEMP):
            zargs = {'z':z, 'redshift_type':TEMPLATE_REDSHIFT_TYPE}
            fnu = self.templates[i].flux_fnu(**zargs)*self.tempfilt.scale[i]
            try:
                tempflux[i, :] = fnu
            except:
                tempflux[i, :] = np.interp(templ.wave,
                                           self.templates[i].wave, fnu)

Here, in order to sum the fluxes for the plot (and for the array templf that's output as part of show_fit), the script checks to see if it can cram the template fluxes into the tempflux array, which is implicitly making an assumption that the wavelength arrays are the same, which is why they have the same lengths. If it can't do that, then the code runs an interpolation to templ.wave, which is set a few lines earlier:

templ = self.templates[0]

As the wavelength of the first template in the list specified by the TEMPLATES_FILE. Once things are interpolated, then the code can actually sum things. However, if your first template in the list has an especially sparse wavelength array, then this has ramifications for the final output SED, so I recommend that the first template file in the list have a pretty fine wavelength array. Notice, in the comment dated June 30th, the first template I specified was eazy_v1.1_sed7.sed (https://github.com/gbrammer/eazy-photoz/blob/master/templates/EAZY_v1.1_lines/eazy_v1.1_sed7.dat), which has only 468 wavelength values, while the others have 2818 wavelength values. By interpolating to the wavelength array from eazy_v1.1_sed7.dat, I was inadvertently lowering the resolution of the final output SED.

Thank you, Kelvin!!! I met with the same problem and your post saved my life :-D !!!

Thanks for your detective work on this @kevinhainline and followup @LittleLin1999 . There are a lot of hacks throughout that should be cleaned up, particularly everywhere there is crudetry/except clause that doesn't catch specific exceptions. Would you be interested in adding a fix and PR for this particular issue with show_fit? A few ideas:

The TemplateGrid object could have a property that is automatically an array of the unique wavelength values of all of the templates in the list
Use Template.resample to interpolate to this grid, either automatically on initialization reading the templates or just inside show_fit, rather than the simple and wrong np.interp interpolation.

Hi @gbrammer Thank you for your reply! I did a test using Hainline+23 templates (https://zenodo.org/records/8092529) and met with a new problem. There are three kinds of wavelength lengths: 468, 2818, and 6900. If I put the 2818 one as the first in the template list, the problem gets resolved. But if I put the 6900 one as the first, the same issue appears, showing extremely large emodel values. So I guess the wavelength grid of the first template cannot be too fine? Is this caused by the same problem in show_fit, too?

Here I attach the examples with the first template being the 6900 grid one and the 2800 one below. test_6900_first test_2800_first

gbrammer / eazy-py

Best-fit template not agreeing with the best-fit photometry #17