possible problem in bandpass_unit_conversion

[NicolettaK]

Hi, I'm working with bandpass integration. In order to check that results are correct I'm doing the following: I integrate a CMB map in two different bands and then check that the integrated maps in thermodynamic units are still equal. This is what I expect for CMB in thermodynamic units also after band integration.

The code I'm using is the following:

 def check_unit_conversion(nuc, deltanu):
    bandpass_frequencies = np.linspace(nuc-deltanu/2., nuc+deltanu/2., 51) * u.GHz
    weights = np.ones(len(bandpass_frequencies))
    sky = pysm3.Sky(nside=128, preset_strings=["c1"])
    CMB_rj_int = sky.get_emission(bandpass_frequencies, weights)
    CMB_thermo_int = CMB_rj_int*pysm3.bandpass_unit_conversion(bandpass_frequencies, weights, u.uK_CMB)
    return CMB_thermo_int

CMB_thermo_int_100 = check_unit_conversion(100, 50)
CMB_thermo_int_400 = check_unit_conversion(400, 50)

if then I plot CMB_thermo_int_100[0]-CMB_thermo_int_400[0] this is what I get:

while what I expect is a map of zeros.

I've done some calculations and I think the problem is in bandpass_unit_conversion, I can share them with you, but first we should agree that this is a good test to do and that I'm actually implementing it in the right way.

[zonca]

bandpass integration is in units of power, so I don't expect to get zeros. Still there could be a bug, we need more detailed test.

[NicolettaK]

Hi @zonca, I think that the CMB in thermodynamic units should always be flat, independently from the units in which you do the bandpass integration. See below:

Given that, I expect to get a map of zeros in the computation described in my previous message.

These equations are similar to the ones implemented by PySM, with the only difference that PySM gets a map in K_RJ and returns a map in K_RJ, meaning that what PySM does is:

which I think is correct!

To then convert the output map into K_CMB units, taking into account the frequency band, I guess that what should be done is the following:

Which is different from what bandpass_unit_conversion does.

[zonca]

isn't the default of pysm.Sky K_rj?

[NicolettaK]

yes it is, and the computation of pysm.Sky which returns the map in K_RJ I think is correct (third equation above). I think there might be a problem in bandpass_unit_conversion which should allow to convert the output of pysm.Sky into K_CMB, taking into account the band

[zonca]

do you understand where the function bandpass_unit_conversion is wrong and can provide a fix?

[NicolettaK]

The function I have implemented is this one:

def bandpass_unit_conversion(freqs, weights, output_unit):
    freqs = check_freq_input(freqs)
    weights_to_rj = (weights * u.uK_RJ).to_value(
            (u.Jy / u.sr), equivalencies=u.cmb_equivalencies(freqs * u.GHz)
        )
    weights_to_out = (weights * output_unit).to_value(
            (u.Jy / u.sr), equivalencies=u.cmb_equivalencies(freqs * u.GHz)
        )
    factor = np.trapz(weights_to_rj, freqs)/np.trapz(weights_to_out, freqs)
    return factor * u.Unit(u.Unit(output_unit) / u.uK_RJ)

It does what I believe is correct, but someone should double check!!

[zonca]

Great thanks, can you please implement that in a pull request and write a unit test?

[NicolettaK]

Ok, I can do that!

[zonca]

fixed in #60

[zonca]

More details on the problem and impact on existing simulations https://zonca.dev/2020/08/debug-bandpass-unit-conversion-pysm3.html

[zonca]

@keskitalo has an implementation of unit conversion for Planck derived from the IDL code. @keskitalo could you please provide CMB2MJysr and MJysr2KRJ for the frequency averaged bandpasses just for HFI?

See https://github.com/healpy/pysm/pull/60/files#diff-ffcf022db138d558b1ea4c56ed5338d4R44-R57 for what I have currently.

[zonca]

Trouble.

tod2flux agrees very well with the Planck published figures, PySM 2 and 3 give different results. I have identified that the difference is a linear weighting in frequency which is applied in tod2flux but not in PySM.

See this snippet from tod2flux:

      correction["mkcmb2mjysr"] = (                                                                          
          1e17                                                                                               
          * scipy.integrate.simps(db_dt * trans, freq)                                                       
          / scipy.integrate.simps(cfreq / freq * trans, freq)                                                
      ) 

I don't understand where cfreq/freq is coming from, @keskitalo do you know if this comes from the definition of the weights in HFI bandpasses?

If I remove this factor, tod2flux agrees with PySM 2 and 3 (they are very close).

Comparison tables for K_CMB to MJy/sr

In percent

The CSV files are available at https://gist.github.com/49f8b05a096d15b5bbd82c94f8c61897

[zonca]

also, the cfreq/freq factor is not normalized, here is its mean for all frequencies in order:

0.8174695583913084
0.8293499066922475
0.8674964664208977
0.9092738281395635
1.330435677498173
1.3648651889530898

[keskitalo]

The calculation of the HFI unit conversion coefficients is written up in https://arxiv.org/abs/1303.5070. The specific case of K_CMB->MJy/sr requires Eq:s (22), (23), (8) and (9). cfreq/freq follows from equating an IRAS spectrum source (slope=-1) with the CMB at the nominal central frequency, cfreq.

If you use as reference any other standard spectrum, you'll find a different unit conversion factor.

[zonca]

ok, thanks @keskitalo I think now it is clear and everything should be consistent in #61, I had to write myself a couple of notebooks to remind me how it worked (https://zonca.dev/2020/09/unit-conversion-broadband-cmb.html and https://zonca.dev/2020/09/unit-conversion-broadband-foregrounds.html) and also added plenty of explanation as comments into the code.