galsci / pysm

PySM 3: Sky emission simulations for Cosmic Microwave Background experiments
https://pysm3.readthedocs.io/
BSD 3-Clause "New" or "Revised" License
36 stars 23 forks source link

CMB only with different scaling in integrated bands #49

Closed jchamilton75 closed 4 years ago

jchamilton75 commented 4 years ago

I am simulating pure CMB and because I want to update my cosmology easily, I use the pysm.CMBMap(nside, filename) with a set of I,Q,U maps I have simulated with synfast from my model (Note that I don't care about having the non-gaussian features from taylens here, so doing what I do is sufficient).

I do it like that:

cmbmap = pysm.CMBMap(nside, map_IQU=filename)
sky.add_component(cmbmap)

and then when I want to get this map from PySM, I do:

themaps_iqu = sky.get_emission([nus_edge[i], nus_edge[i + 1]] * u.GHz)

where nus_edges defin my bandpass, I then convert to uK_CMB with:

mysky[i, :, :] = np.array(themaps_iqu.to(u.uK_CMB, equivalencies=u.cmb_equivalencies(nus_in[i] * u.GHz))).T

where nus_in[i] is the center of my bandpass

Doing so gives pixels with a few percent difference from one band to another and they don;'t match the initial CMB maps. So I probably do something wrong here...

I could get consistent results by using single frequency instead of integration but I'm afraid this is not so good when I will be including something else but CMB...

zonca commented 4 years ago

bandpasses are assumed to be in power units, MJ/sr, so it is expected that bandpass-integrated emissions have a discrepancy with a delta-bandpass at the center.

jchamilton75 commented 4 years ago

OK thanks I understand, so I will weight them accordingly. Thanks a lot

jchamilton75 commented 4 years ago

Sorry Andrea, me again... I'm afraid I could not find a way to weight my band integration properly in order to get exactly the input CMB in each band... Could you help me with this ?

I did the following:

        nfreqinteg = 5
        nus = np.linspace(nus_edge[i], nus_edge[i + 1], nfreqinteg)
        freqs = utils.check_freq_input(nus)

        convert_to_uK_RJ = (np.ones(len(freqs), dtype=np.double) * u.uK_CMB).to_value(
        u.uK_RJ, equivalencies=u.cmb_equivalencies(freqs))
        print('Convert_to_uK_RJ :',convert_to_uK_RJ)

        weights = np.ones(nfreqinteg) * convert_to_uK_RJ

        ### Integrate through band using filter shape defined in weights
        themaps_iqu = self.sky.get_emission(nus * u.GHz, weights=weights)
        sky[i, :, :] = np.array(themaps_iqu.to(u.uK_CMB, equivalencies=u.cmb_equivalencies(nus_in[i] * u.GHz))).T
        ratio = np.mean(self.input_cmb_maps[0,:]/sky[i,:,0])
        print('Ratio to initial: ',ratio)

Which leads to a 0.5% difference bewteen PySM CMB maps and thos used as an input:

Ratio to initial: 1.0051389388185639

zonca commented 4 years ago

I don't have time to do user support, try to do a uniform weighting in MJ/sr