diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 5e9c30c..692e3d9 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -1,6 +1,9 @@ v0.x.x ============= * Current version. +* Add dichroic detector (dichroic branch) +* Add subscans information +* Add beam ellipticity systematics v0.6.1 ============= diff --git a/s4cmb/input_sky.py b/s4cmb/input_sky.py index e1b8b97..5b569c8 100644 --- a/s4cmb/input_sky.py +++ b/s4cmb/input_sky.py @@ -21,9 +21,10 @@ class HealpixFitsMap to handle it easily. class HealpixFitsMap(): """ Class to handle fits file containing healpix maps """ def __init__(self, input_filename, - do_pol=True, verbose=False, fwhm_in=0.0, nside_in=16, - lmax=None, map_seed=53543, no_ileak=False, no_quleak=False, - compute_derivatives=False, ext_map_gal=False): + do_pol=True, verbose=False, fwhm_in=0.0, fwhm_in2=None, + nside_in=16, lmax=None, map_seed=53543, no_ileak=False, + no_quleak=False, compute_derivatives=False, + ext_map_gal=False): """ Parameters @@ -43,6 +44,10 @@ def __init__(self, input_filename, If input_filename is a CAMB lensed cl file, the generated maps will be convolved with a beam having this fwhm_in. In arcmin. No effect if you provide maps directly. + fwhm_in2 : float, optional + If provided, will generate another set of I, Q, U with this + resolution (useful for dichroic detectors). Default is None. + No effect if you provide maps directly. nside_in : int, optional If input_filename is a CAMB lensed cl file, the maps will be generated at a resolution nside_in. No effect if you provide maps @@ -77,6 +82,7 @@ def __init__(self, input_filename, self.no_quleak = no_quleak self.ext_map_gal = ext_map_gal self.fwhm_in = fwhm_in + self.fwhm_in2 = fwhm_in2 self.nside_in = nside_in if lmax is None: self.lmax = 2 * self.nside_in @@ -89,6 +95,10 @@ def __init__(self, input_filename, self.Q = None self.U = None + self.I2 = None + self.Q2 = None + self.U2 = None + fromalms = False if type(self.input_filename) == list: if self.verbose: @@ -120,6 +130,8 @@ def load_healpix_fits_map(self, force=False): """ Load from disk into memory a sky map. + Not updated for dichroic for the moment. + Parameters ---------- force : bool @@ -187,6 +199,13 @@ def load_healpix_fits_map_from_alms(self, force=False): But you can force it >>> hpmap.load_healpix_fits_map_from_alms(force=True) + + You can also generate 2 sets of maps with different resolution which + is useful for dichroic detectors + >>> hpmap = HealpixFitsMap(input_filename=filenames, fwhm_in=3.5, + ... fwhm_in2=1.8, nside_in=16,) + >>> hasattr(hpmap, 'Q2') + True """ if self.I is None or force: if self.do_pol: @@ -199,16 +218,35 @@ def load_healpix_fits_map_from_alms(self, force=False): nside=self.nside_in, pixwin=False, fwhm=self.fwhm_in / 60. * np.pi / 180., - sigma=None, pol=True, inplace=False, verbose=self.verbose) + sigma=None, pol=True, inplace=False, + verbose=self.verbose) + if self.fwhm_in2 is not None: + self.I2, self.Q2, self.U2 = hp.alm2map( + [tlm, elm, blm], + nside=self.nside_in, + pixwin=False, + fwhm=self.fwhm_in2 / 60. * np.pi / 180., + sigma=None, pol=True, inplace=False, + verbose=self.verbose) else: - self.I = hp.read_alm(self.input_filename[0]) + tlm = hp.read_alm(self.input_filename[0]) self.I = hp.alm2map( tlm, nside=self.nside_in, pixwin=False, fwhm=self.fwhm_in / 60. * np.pi / 180., - sigma=None, pol=False, inplace=False, verbose=self.verbose) + sigma=None, pol=False, inplace=False, + verbose=self.verbose) + if self.fwhm_in2 is not None: + self.I2 = hp.alm2map( + tlm, + nside=self.nside_in, + pixwin=False, + fwhm=self.fwhm_in2 / 60. * np.pi / 180., + sigma=None, pol=False, inplace=False, + verbose=self.verbose) + self.nside = hp.npix2nside(len(self.I)) else: print("External data already present in memory") @@ -221,8 +259,8 @@ def create_healpix_fits_map(self, force=False): Parameters ---------- force : bool - If true, force to recreate the maps in memory even if it is already - loaded. Default is False. + If true, force to recreate the maps in memory even + if it is already loaded. Default is False. Examples ---------- @@ -239,20 +277,44 @@ def create_healpix_fits_map(self, force=False): But you can force it >>> hpmap.create_healpix_fits_map(force=True) + + You can also load 2 sets of maps with different resolution, which + is useful for dichroic detectors + >>> filename = 's4cmb/data/test_data_set_lensedCls.dat' + >>> hpmap = HealpixFitsMap(input_filename=filename, fwhm_in=3.5, + ... fwhm_in2=1.8, nside_in=16, map_seed=489237) + >>> hasattr(hpmap, 'I2') + True """ if self.I is None or force: if self.do_pol: - self.I, self.Q, self.U = create_sky_map(self.input_filename, - nside=self.nside_in, - FWHM=self.fwhm_in, - seed=self.map_seed, - lmax=self.lmax) + self.I, self.Q, self.U = create_sky_map( + self.input_filename, + nside=self.nside_in, + FWHM=self.fwhm_in, + seed=self.map_seed, + lmax=self.lmax) + if self.fwhm_in2 is not None: + self.I2, self.Q2, self.U2 = create_sky_map( + self.input_filename, + nside=self.nside_in, + FWHM=self.fwhm_in2, + seed=self.map_seed, + lmax=self.lmax) else: - self.I = create_sky_map(self.input_filename, - nside=self.nside_in, - FWHM=self.fwhm_in, - seed=self.map_seed, - lmax=self.lmax) + self.I = create_sky_map( + self.input_filename, + nside=self.nside_in, + FWHM=self.fwhm_in, + seed=self.map_seed, + lmax=self.lmax) + if self.fwhm_in2 is not None: + self.I2 = create_sky_map( + self.input_filename, + nside=self.nside_in, + FWHM=self.fwhm_in2, + seed=self.map_seed, + lmax=self.lmax) self.nside = hp.npix2nside(len(self.I)) else: print("External data already present in memory") @@ -274,10 +336,20 @@ def set_leakage_to_zero(self): >>> hpmap = HealpixFitsMap('myfits_to_test_.fits', no_quleak=True) >>> print(hpmap.Q, hpmap.U) [ 0. 0. 0. ..., 0. 0. 0.] [ 0. 0. 0. ..., 0. 0. 0.] + + If you have two sets of maps, it will remove leakages from the two sets + >>> filename = 's4cmb/data/test_data_set_lensedCls.dat' + >>> hpmap = HealpixFitsMap(input_filename=filename, fwhm_in=3.5, + ... fwhm_in2=1.8, nside_in=16, map_seed=489237, + ... no_ileak=True, no_quleak=True) + >>> print(hpmap.I, hpmap.I2) + [ 0. 0. 0. ..., 0. 0. 0.] [ 0. 0. 0. ..., 0. 0. 0.] """ ## Set temperature to zero to avoid I->QU leakage if self.no_ileak: self.I[:] = 0.0 + if self.I2 is not None: + self.I2[:] = 0.0 ## Set polarisation to zero to avoid QU leakage if self.no_quleak: @@ -286,6 +358,11 @@ def set_leakage_to_zero(self): if self.U is not None: self.U[:] = 0.0 + if self.Q2 is not None: + self.Q2[:] = 0.0 + if self.U2 is not None: + self.U2[:] = 0.0 + def compute_intensity_derivatives(self, fromalm=False): """ Compute derivatives of the input temperature map (healpix). @@ -293,6 +370,8 @@ def compute_intensity_derivatives(self, fromalm=False): higher than 1 (see healpix for better treatment), so we use only an approximation. + Not updated for dichroic for the moment. + Parameters ---------- fromalm : bool, optional @@ -453,11 +532,12 @@ def create_sky_map(cl_fn, nside=16, FWHM=0.0, seed=548397, lmax=None): FWHM_rad = FWHM / 60. * np.pi / 180. np.random.seed(seed) - I, Q, U = hp.synfast([TT / llp, EE / llp, BB / llp, TE / llp], nside, - lmax=lmax, mmax=None, alm=False, - pol=True, pixwin=False, - fwhm=FWHM_rad, sigma=None, new=True, - verbose=False) + I, Q, U = hp.synfast( + [TT / llp, EE / llp, BB / llp, TE / llp], nside, + lmax=lmax, mmax=None, alm=False, + pol=True, pixwin=False, + fwhm=FWHM_rad, sigma=None, new=True, + verbose=False) return I, Q, U def write_healpix_cmbmap(output_filename, data, fits_IDL=False, @@ -504,13 +584,15 @@ def write_healpix_cmbmap(output_filename, data, fits_IDL=False, ## Need to introduce this workaround because the last ## version of healpy introduced non-backward compatibility... try: - hp.write_map(output_filename, data, fits_IDL=fits_IDL, - coord=coord, column_names=None, partial=partial, - extra_header=extra_header, overwrite=True) + hp.write_map( + output_filename, data, fits_IDL=fits_IDL, + coord=coord, column_names=None, partial=partial, + extra_header=extra_header, overwrite=True) except: - hp.write_map(output_filename, data, fits_IDL=fits_IDL, - coord=coord, column_names=None, partial=partial, - extra_header=extra_header) + hp.write_map( + output_filename, data, fits_IDL=fits_IDL, + coord=coord, column_names=None, partial=partial, + extra_header=extra_header) def write_dummy_map(filename='myfits_to_test_.fits', nside=16): """ diff --git a/s4cmb/instrument.py b/s4cmb/instrument.py index 48f7aca..bc77d2e 100644 --- a/s4cmb/instrument.py +++ b/s4cmb/instrument.py @@ -449,12 +449,22 @@ def __init__(self, Default is 58347. projected_fp_size : float, optional Size of the focal plane on the sky (in degree). This has to - do with the size of the mirror. Default = 3 degrees. + do with the size of the mirror and the beam size. + Default = 3 degrees. pm_name : string, optional The pointing model to load. Currently, only the five-parameter pointing model (Mangum 2001) is implemented (pm_name = 5params). verbose : boolean, optional If True, print out a number of useful comments for verboseging. + type_hwp : string, optional + Choose between CRWHP or stepped. + freq_hwp : float, optional + If type_hwp=CRHWP, then freq_hwp sets the spin frequency of the + HWP [Hz]. + angle_hwp : float, optional + If type_hwp=CRHWP, angle_hwp corresponds to the starting position. + If type_hwp=stepped, angle_hwp corresponds to the step size + which is updated after each scan (CES). [deg] Examples ---------- @@ -471,6 +481,48 @@ def __init__(self, self.half_wave_plate = HalfWavePlate(type_hwp, freq_hwp, angle_hwp) + def make_dichroic(self, fwhm=1.8, beam_seed=58347, projected_fp_size=3., + shift_angle=45): + """ + Add a layer of detectors on top of the existing ones + to have dichroic detectors. The new detectors will sit on top of the + existing ones (that is same number and same location on the focal + plane). The only difference is the beam size for the moment (and the + projected size of the fp on sky). + + Parameters + ---------- + fwhm : float, optional + Full Width Half Maximum of the beam (in arcmin). + Default = 3.5 arcmin. + beam_seed : int + Seed used to generate angle of rotation of beam axes. + Default is 58347. + projected_fp_size : float, optional + Size of the focal plane on the sky (in degree). This has to + do with the size of the mirror and the beam size. + Default = 3 degrees. + shift_angle : int, optional + Shift the polarisation angle with respect to the ones from the + first frequency channel. + + Examples + ---------- + >>> instrument = Hardware() + >>> instrument.make_dichroic() + >>> hasattr(instrument, 'focal_plane2') + True + """ + self.focal_plane2 = copy.copy(self.focal_plane) + + ## Shift the polarisation angles (and clip it btw 0 and 360 deg) + self.focal_plane2.bolo_polangle = ( + np.array(self.focal_plane2.bolo_polangle) + shift_angle) % 360 + + self.beam_model2 = BeamModel( + self.focal_plane2, fwhm, beam_seed, + projected_fp_size) + class FocalPlane(): """ Class to handle the focal plane of the instrument. """ def __init__(self, diff --git a/s4cmb/systematics.py b/s4cmb/systematics.py index b535b83..a4bc9c1 100644 --- a/s4cmb/systematics.py +++ b/s4cmb/systematics.py @@ -424,7 +424,8 @@ def modify_pointing_parameters(values, errors): values_mod = [p + err for p, err in zip(values, errors)] return values_mod -def step_function(nbolos, nsamples, mean=1, std=0.05, nbreaks=1, seed=0): +def step_function(nbolos, nsamples, mean=1, std=0.05, + nbreaks=1, sign='same', seed=0): """ Generate step functions for each bolometer from 1 to a values drawn from N(mean, std). The full timestream is broken into nbreaks @@ -445,6 +446,10 @@ def step_function(nbolos, nsamples, mean=1, std=0.05, nbreaks=1, seed=0): Default is 0.05 (that is \pm 5% with respect to a mean=1). nbreaks : int Number of break (number of retuning). + sign : string + If same, both detectors in the pair will have the same gain. + If opposite, detectors will have gains with opposite sign + (differential gain). Returns ---------- @@ -458,6 +463,9 @@ def step_function(nbolos, nsamples, mean=1, std=0.05, nbreaks=1, seed=0): >>> gains = step_function(nbolos, nsamples, nbreaks=1) >>> print(gains[0]) [ 1. 1. 1.08820262 1.08820262] + + >>> gains = step_function(nbolos, nsamples, nbreaks=1, sign='opposite') + >>> assert gains[0][0] == 2 - gains[1][0] """ ## Fix the seed state = np.random.RandomState(seed) @@ -483,9 +491,13 @@ def step_function(nbolos, nsamples, mean=1, std=0.05, nbreaks=1, seed=0): else: continue + if sign == 'opposite': + gains[1] = 2 - gains[1] + return gains -def step_function_gen(nsamples, mean=1, std=0.05, nbreaks=1, seed=0): +def step_function_gen(nsamples, mean=1, std=0.05, + nbreaks=1, sign='same', seed=0): """ Generator of step functions for each bolometer from 1 to a values drawn from N(mean, std). The full timestream is broken into nbreaks @@ -506,6 +518,10 @@ def step_function_gen(nsamples, mean=1, std=0.05, nbreaks=1, seed=0): Default is 0.05 (that is \pm 5% with respect to a mean=1). nbreaks : int Number of break (number of retuning). + sign : string + If same, both detectors in the pair will have the same gain. + If opposite, detectors will have gains with opposite sign + (differential gain). Returns ---------- @@ -518,6 +534,10 @@ def step_function_gen(nsamples, mean=1, std=0.05, nbreaks=1, seed=0): >>> gains_gen = step_function_gen(nsamples, nbreaks=1) >>> print(next(gains_gen)[0]) [ 1. 1. 1.08820262 1.08820262] + + >>> gains_gen = step_function_gen(nsamples, nbreaks=1, sign='opposite') + >>> g = next(gains_gen) + >>> assert g[0][0] == 2 - g[1][0] """ ## Fix the seed state = np.random.RandomState(seed) @@ -544,9 +564,13 @@ def step_function_gen(nsamples, mean=1, std=0.05, nbreaks=1, seed=0): else: continue + if sign == 'opposite': + gains[1] = 2 - gains[1] + yield gains -def linear_function(nbolos, nsamples, mean=1, std=0.05, nbreaks=1, seed=0): +def linear_function(nbolos, nsamples, mean=1, std=0.05, + nbreaks=1, sign='same', seed=0): """ Generate linear functions for each bolometer from 1 to a values drawn from N(mean, std). The full timestream is broken into nbreaks @@ -567,6 +591,10 @@ def linear_function(nbolos, nsamples, mean=1, std=0.05, nbreaks=1, seed=0): Default is 0.05 (that is \pm 5% with respect to a mean=1). nbreaks : int Number of break (number of retuning). + sign : string + If same, both detectors in the pair will have the same gain. + If opposite, detectors will have gains with opposite sign + (differential gain). Returns ---------- @@ -580,6 +608,9 @@ def linear_function(nbolos, nsamples, mean=1, std=0.05, nbreaks=1, seed=0): >>> gains = linear_function(nbolos, nsamples, nbreaks=1) >>> print(gains[0]) [ 1. 1.02940087 1.05880174 1.08820262] + + >>> gains = linear_function(nbolos, nsamples, nbreaks=1, sign='opposite') + >>> assert gains[0][0] == 2 - gains[1][0] """ ## Fix the seed state = np.random.RandomState(seed) @@ -606,9 +637,13 @@ def linear_function(nbolos, nsamples, mean=1, std=0.05, nbreaks=1, seed=0): else: continue + if sign == 'opposite': + gains[1] = 2 - gains[1] + return gains -def linear_function_gen(nsamples, mean=1, std=0.05, nbreaks=1, seed=0): +def linear_function_gen(nsamples, mean=1, std=0.05, + nbreaks=1, sign='same', seed=0): """ Generator returning linear functions for each bolometer from 1 to a values drawn from N(mean, std). The full timestream is broken into nbreaks @@ -629,6 +664,10 @@ def linear_function_gen(nsamples, mean=1, std=0.05, nbreaks=1, seed=0): Default is 0.05 (that is \pm 5% with respect to a mean=1). nbreaks : int Number of break (number of retuning). + sign : string + If same, both detectors in the pair will have the same gain. + If opposite, detectors will have gains with opposite sign + (differential gain). Returns ---------- @@ -641,6 +680,11 @@ def linear_function_gen(nsamples, mean=1, std=0.05, nbreaks=1, seed=0): >>> gains_gen = linear_function_gen(nsamples, nbreaks=1) >>> print(next(gains_gen)[0]) [ 1. 1.02940087 1.05880174 1.08820262] + + >>> gains_gen = linear_function_gen(nsamples, nbreaks=1, sign='opposite') + >>> g = next(gains_gen) + >>> assert g[0][0] == 2 - g[1][0] + """ ## Fix the seed state = np.random.RandomState(seed) @@ -668,6 +712,9 @@ def linear_function_gen(nsamples, mean=1, std=0.05, nbreaks=1, seed=0): else: continue + if sign == 'opposite': + gains[1] = 2 - gains[1] + yield gains def get_kernel_coefficients(beamprm, pairlist, nx=128, pix_size=None): diff --git a/s4cmb/tod.py b/s4cmb/tod.py index 520c6d4..ae29277 100644 --- a/s4cmb/tod.py +++ b/s4cmb/tod.py @@ -39,7 +39,8 @@ def __init__(self, hardware, scanning_strategy, HealpixFitsMap, nside_out=None, pixel_size=None, width=140., cut_pixels_outside=True, array_noise_level=None, array_noise_seed=487587, - mapping_perpair=False, verbose=False): + array_noise_level2=None, array_noise_seed2=56736, + mapping_perpair=False, mode='standard', verbose=False): """ C'est parti! @@ -94,6 +95,11 @@ def __init__(self, hardware, scanning_strategy, HealpixFitsMap, If True, assume that you want to process pairs of bolometers one-by-one, that is pairs are uncorrelated. Default is False (and should be False unless you know what you are doing). + mode : string, optional + Choose between `standard` (1 frequency band) and `dichroic` + (2 frequency bands). If `dichroic` is chosen, make sure your + hardware can handle it (see instrument.py) and HealpixFitsMap + should contain the inputs maps at different frequency. """ ## Initialise args self.verbose = verbose @@ -101,14 +107,30 @@ def __init__(self, hardware, scanning_strategy, HealpixFitsMap, self.scanning_strategy = scanning_strategy self.HealpixFitsMap = HealpixFitsMap self.mapping_perpair = mapping_perpair + + ## Check if you can run dichroic detectors + self.mode = mode + if self.mode == 'dichroic' and ( + not hasattr(self.HealpixFitsMap, 'I2')): + raise IOError("You need two sets of maps for dichroic detectors!") + if self.mode == 'dichroic' and ( + not hasattr(self.hardware, 'focal_plane2')): + raise IOError("You need two sets of det for dichroic detectors!") + if self.mode == 'dichroic' and ( + not hasattr(self.hardware, 'beam_model2')): + raise IOError("You need two sets of det for dichroic detectors!") + self.width = width self.cut_pixels_outside = cut_pixels_outside + + ## Check the projection self.projection = projection assert self.projection in ['healpix', 'flat'], \ ValueError("Projection <{}> for ".format(self.projection) + "the output map not understood! " + "Choose among ['healpix', 'flat'].") + ## Check if the number of CES is consistent with the scanning strategy self.CESnumber = CESnumber assert self.CESnumber < self.scanning_strategy.nces, \ ValueError("The scan index must be between 0 and {}.".format( @@ -182,6 +204,17 @@ def __init__(self, hardware, scanning_strategy, HealpixFitsMap, else: self.noise_generator = None + self.array_noise_level2 = array_noise_level2 + self.array_noise_seed2 = array_noise_seed2 + if self.array_noise_level2 is not None: + self.noise_generator2 = WhiteNoiseGenerator( + array_noise_level=self.array_noise_level2, + ndetectors=2*self.npair, + ntimesamples=self.nsamples, + array_noise_seed=self.array_noise_seed2) + else: + self.noise_generator2 = None + def get_angles(self): """ Retrieve polarisation angles: intrinsic (focal plane) and HWP angles, @@ -192,13 +225,22 @@ def get_angles(self): size=self.nsamples) self.intrinsic_polangle = self.hardware.focal_plane.bolo_polangle + self.intrinsic_polangle2 = None + if self.mode == 'dichroic': + self.intrinsic_polangle2 = self.hardware.focal_plane2.bolo_polangle ## Will contain the total polarisation angles for all bolometers ## That is PA + intrinsic + 2 * HWP if not self.mapping_perpair: self.pol_angs = np.zeros((self.npair, self.nsamples)) + self.pol_angs2 = None + if self.mode == 'dichroic': + self.pol_angs2 = np.zeros((self.npair, self.nsamples)) else: self.pol_angs = np.zeros((1, self.nsamples)) + self.pol_angs2 = None + if self.mode == 'dichroic': + self.pol_angs2 = np.zeros((1, self.nsamples)) def get_timestream_masks(self): """ @@ -318,7 +360,7 @@ def get_weights(self): else: return np.ones((2, 1), dtype=int) - def set_detector_gains(self, new_gains=None): + def set_detector_gains(self, new_gains=None, new_gains2=None): """ Set the gains of the detectors (unitless). Default is 1., that is perfectly calibrated. @@ -327,7 +369,10 @@ def set_detector_gains(self, new_gains=None): ---------- new_gains : 1d array Array containing the gain value for all detectors - (1 number per detector). + (1 number per detector) for the 1st frequency channel. + new_gains2 : 1d array + Array containing the gain value for all detectors + (1 number per detector) for the 2nd frequency channel. Examples ---------- @@ -341,6 +386,17 @@ def set_detector_gains(self, new_gains=None): >>> tod.set_detector_gains(new_gains=new_gains) >>> print(tod.gain[0]) 2.0 + + For dichroic + >>> inst, scan, sky_in = load_fake_instrument(fwhm_in2=1.8) + >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in, + ... mode='dichroic', CESnumber=1) + >>> assert tod.gain2 is not None + + Change the value of gains for both frequency channels + >>> new_gains = np.ones(2 * tod.npair) * 2. + >>> new_gains2 = np.ones(2 * tod.npair) * 4. + >>> tod.set_detector_gains(new_gains=new_gains, new_gains2=new_gains2) """ if new_gains is not None: assert len(new_gains) == 2 * self.npair, \ @@ -350,7 +406,18 @@ def set_detector_gains(self, new_gains=None): else: self.gain = np.ones(2 * self.npair) - def set_detector_gains_perpair(self, new_gains=None): + self.gain2 = None + if self.mode == 'dichroic': + if new_gains2 is not None: + assert len(new_gains2) == 2 * self.npair, \ + ValueError( + "You have to provide {} new gain values!".format( + 2 * self.npair)) + self.gain2 = new_gains2 + else: + self.gain2 = np.ones(2 * self.npair) + + def set_detector_gains_perpair(self, new_gains=None, new_gains2=None): """ Set the gains of all 2 pair detectors for each timestep (unitless). This is particularly useful to introduce drifts for example. @@ -360,7 +427,10 @@ def set_detector_gains_perpair(self, new_gains=None): ---------- new_gains : 2d array of size (2, nsamples) Array containing the gain value for all detectors - (nsamples number per detector). + (nsamples number per detector) for the 1st frequency channel. + new_gains2 : 2d array of size (2, nsamples) + Array containing the gain value for all detectors + (nsamples number per detector) for the 2nd frequency channel. Examples ---------- @@ -375,6 +445,19 @@ def set_detector_gains_perpair(self, new_gains=None): >>> tod.set_detector_gains_perpair(new_gains=new_gains) >>> print(tod.gain[0][0:4]) [ 2. 1. 2. 1.] + + Dichroic + >>> inst, scan, sky_in = load_fake_instrument(fwhm_in2=1.8) + >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in, + ... mode='dichroic', CESnumber=1) + + Change the value of gains every other sample. + >>> new_gains = np.ones((2, tod.nsamples)) + >>> new_gains[:, ::2] = 2. + >>> new_gains2 = np.ones((2, tod.nsamples)) + >>> new_gains2[:, ::2] = 4. + >>> tod.set_detector_gains_perpair( + ... new_gains=new_gains, new_gains2=new_gains2) """ if new_gains is not None: msg = "You have to provide ({}, {}) new gain values!" @@ -384,6 +467,16 @@ def set_detector_gains_perpair(self, new_gains=None): else: self.gain = np.ones((2, self.nsamples)) + self.gain2 = None + if self.mode == 'dichroic': + if new_gains2 is not None: + msg = "You have to provide ({}, {}) new gain values!" + assert new_gains2.shape == (2, self.nsamples), \ + ValueError(msg.format(2, self.nsamples)) + self.gain2 = new_gains2 + else: + self.gain2 = np.ones((2, self.nsamples)) + def get_boresightpointing(self): """ Initialise the boresight pointing for all the focal plane bolometers. @@ -455,10 +548,20 @@ def compute_simpolangle(self, ch, parallactic_angle, polangle_err=False): ---------- >>> inst, scan, sky_in = load_fake_instrument() >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in, CESnumber=0) - >>> angles = tod.compute_simpolangle(ch=0, + >>> angles, angles2 = tod.compute_simpolangle(ch=0, ... parallactic_angle=np.array([np.pi] * tod.nsamples)) >>> print(angles[:4]) [ 0. 0.31415927 0.62831853 0.9424778 ] + >>> assert angles2 is None + + If you have dichroic detectors, the routine will return angles for + the two frequency channels + >>> inst, scan, sky_in = load_fake_instrument(fwhm_in2=1.8) + >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in, + ... mode='dichroic', CESnumber=0) + >>> angles, angles2 = tod.compute_simpolangle(ch=0, + ... parallactic_angle=np.array([np.pi] * tod.nsamples)) + >>> assert angles2 is not None """ if not polangle_err: ang_pix = (90.0 - self.intrinsic_polangle[ch]) * d2r @@ -469,14 +572,25 @@ def compute_simpolangle(self, ch, parallactic_angle, polangle_err=False): pol_ang = parallactic_angle + ang_pix + 2.0 * self.hwpangle else: pol_ang = parallactic_angle - ang_pix - 2.0 * self.hwpangle + + pol_ang2 = None + if self.mode == 'dichroic': + ang_pix2 = (90.0 - self.intrinsic_polangle2[ch]) * d2r + + ## Demodulation or pair diff use different convention + ## for the definition of the angle. + if not hasattr(self, 'dm'): + pol_ang2 = parallactic_angle + ang_pix2 + 2.0*self.hwpangle + else: + pol_ang2 = parallactic_angle - ang_pix2 - 2.0*self.hwpangle else: print("This is where you call the systematic module!") sys.exit() pass - return pol_ang + return pol_ang, pol_ang2 - def return_parallactic_angle(self, ch): + def return_parallactic_angle(self, ch, frequency_channel=1): """ Return the parallactic angles (orientations of the pixel on sky) for one specific channel (full timestream). @@ -485,24 +599,52 @@ def return_parallactic_angle(self, ch): ---------- ch : int Index of the bolometer (within the focal plane). + frequency_channel : int, optional + When using dichroic detectors, you need to specify the channel. Returns ---------- pa : 1d array - Parallactic angles (timestream) for detector ch. + Parallactic angles (timestream) for detector ch. [radian] + + Examples + ---------- + Simple TES bolometers + >>> inst, scan, sky_in = load_fake_instrument() + >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in, CESnumber=1) + >>> pa = tod.return_parallactic_angle(0) + + Dichroic detectors + >>> inst, scan, sky_in = load_fake_instrument(fwhm_in2=1.8) + >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in, + ... mode='dichroic', CESnumber=1) + >>> pa1 = tod.return_parallactic_angle(0, frequency_channel=1) + >>> pa2 = tod.return_parallactic_angle(0, frequency_channel=2) """ + if frequency_channel == 1: + pa = self.pol_angs + elif frequency_channel == 2: + pa = self.pol_angs2 + if self.mapping_perpair is True: - ang_pix = (90.0 - self.intrinsic_polangle[ch]) * d2r + if frequency_channel == 1: + ang_pix = (90.0 - self.intrinsic_polangle[ch]) * d2r + elif frequency_channel == 2: + ang_pix = (90.0 - self.intrinsic_polangle2[ch]) * d2r + if not hasattr(self, 'dm'): - return self.pol_angs - ang_pix - 2*self.hwpangle + return pa - ang_pix - 2*self.hwpangle else: - return self.pol_angs + ang_pix + return pa + ang_pix else: - ang_pix = (90.0 - self.intrinsic_polangle[2*ch]) * d2r + if frequency_channel == 1: + ang_pix = (90.0 - self.intrinsic_polangle[2*ch]) * d2r + elif frequency_channel == 2: + ang_pix = (90.0 - self.intrinsic_polangle2[2*ch]) * d2r if not hasattr(self, 'dm'): - return self.pol_angs[ch] - ang_pix - 2*self.hwpangle + return pa[ch] - ang_pix - 2*self.hwpangle else: - return self.pol_angs[ch] + ang_pix + return pa[ch] + ang_pix def map2tod(self, ch): """ @@ -517,17 +659,35 @@ def map2tod(self, ch): Returns ---------- - ts : 1d array + ts : 1d array or array of 1d array The timestream for detector ch. If `self.HealpixFitsMap.do_pol` is True it returns intensity+polarisation, otherwise just intensity. + If dichroic, ts = np.array([ts_freq1, ts_freq2]). Examples ---------- + Simple TES bolometers >>> inst, scan, sky_in = load_fake_instrument() >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in, CESnumber=1) >>> d = tod.map2tod(0) - >>> print(round(d[0], 3)) #doctest: +NORMALIZE_WHITESPACE - -42.874 + >>> print(d.shape) + (115200,) + + Dichroic detectors + >>> inst, scan, sky_in = load_fake_instrument(fwhm_in2=1.8) + >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in, + ... mode='dichroic', CESnumber=1) + >>> d = tod.map2tod(0) + >>> print(d.shape) + (2, 115200) + + Dichroic detectors with different noise levels + >>> inst, scan, sky_in = load_fake_instrument(fwhm_in2=1.8) + >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in, + ... array_noise_level=2.5, array_noise_seed=487587, + ... array_noise_level2=25., array_noise_seed2=56736, + ... mode='dichroic', CESnumber=1) + >>> d = tod.map2tod(0) """ ## Use bolometer beam offsets. azd, eld = self.xpos[ch], self.ypos[ch] @@ -586,8 +746,15 @@ def map2tod(self, ch): else: noise = 0.0 + if self.noise_generator2 is not None: + noise2 = self.noise_generator2.simulate_noise_one_detector(ch) + else: + noise2 = 0.0 + if self.HealpixFitsMap.do_pol: - pol_ang = self.compute_simpolangle(ch, pa, polangle_err=False) + ## pol_ang2 is None if mode == 'standard' + pol_ang, pol_ang2 = self.compute_simpolangle( + ch, pa, polangle_err=False) ## For demodulation, HWP angles are not included at the level ## of the pointing matrix (convention). @@ -602,14 +769,47 @@ def map2tod(self, ch): elif ch % 2 == 0 and self.mapping_perpair: self.pol_angs[0] = pol_ang_out - return (self.HealpixFitsMap.I[index_global] + - self.HealpixFitsMap.Q[index_global] * np.cos(2 * pol_ang) + - sign * self.HealpixFitsMap.U[index_global] * - np.sin(2 * pol_ang) + noise) * norm + ts1 = ( + self.HealpixFitsMap.I[index_global] + + self.HealpixFitsMap.Q[index_global] * np.cos(2 * pol_ang) + + sign * self.HealpixFitsMap.U[index_global] * + np.sin(2 * pol_ang) + noise) * norm + + if self.mode == 'standard': + return ts1 + + elif self.mode == 'dichroic': + # For demodulation, HWP angles are not included at the level + ## of the pointing matrix (convention). + if hasattr(self, 'dm'): + pol_ang_out2 = pol_ang2 + 2.0 * self.hwpangle + else: + pol_ang_out2 = pol_ang2 + + ## Store list polangle only for top bolometers + if ch % 2 == 0 and not self.mapping_perpair: + self.pol_angs2[int(ch/2)] = pol_ang_out2 + elif ch % 2 == 0 and self.mapping_perpair: + self.pol_angs2[0] = pol_ang_out2 + + ts2 = ( + self.HealpixFitsMap.I2[index_global] + + self.HealpixFitsMap.Q2[index_global] * np.cos(2*pol_ang2) + + sign * self.HealpixFitsMap.U2[index_global] * + np.sin(2 * pol_ang2) + noise2) * norm + return np.array([ts1, ts2]) + else: - return norm * (self.HealpixFitsMap.I[index_global] + noise) + ts1 = norm * (self.HealpixFitsMap.I[index_global] + noise) + if self.mode == 'standard': + return ts1 + elif self.mode == 'dichroic': + ts2 = norm * (self.HealpixFitsMap.I2[index_global] + noise2) + return np.array([ts1, ts2]) - def tod2map(self, waferts, output_maps, gdeprojection=False): + def tod2map(self, waferts, output_maps, + gdeprojection=False, + frequency_channel=1): """ Project time-ordered data into sky maps for the whole array. Maps are updated on-the-fly. Massive speed-up thanks to the @@ -626,6 +826,9 @@ def tod2map(self, waferts, output_maps, gdeprojection=False): If True, perform a deprojection of a constant contribution in the polarisation timestream: d^{-} = G + Qcos(2*theta) + U*sin(2theta). Work only for pair difference. + frequency_channel : int, optional + If you are processing dichroic pixels, you need to specify the + index of the frequency channel (1 or 2). Default is 1. Examples ---------- @@ -672,7 +875,37 @@ def tod2map(self, waferts, output_maps, gdeprojection=False): ... npixsky=tod.npixsky, pixel_size=tod.pixel_size) >>> tod.tod2map(d, m) + HEALPIX + Dichroic + >>> inst, scan, sky_in = load_fake_instrument(fwhm_in2=1.8) + >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in, mode='dichroic', + ... CESnumber=0, projection='healpix', mapping_perpair=True) + >>> m1 = OutputSkyMap(projection=tod.projection, + ... nside=tod.nside_out, obspix=tod.obspix) + >>> m2 = OutputSkyMap(projection=tod.projection, + ... nside=tod.nside_out, obspix=tod.obspix) + >>> for pair in tod.pair_list: + ... d = np.array([tod.map2tod(det) for det in pair]) + ... d1 = d[:, 0] ## first frequency channel + ... d2 = d[:, 1] ## second frequency channel + ... tod.tod2map(d1, m1) + ... tod.tod2map(d2, m2) + + HEALPIX + deprojection + >>> inst, scan, sky_in = load_fake_instrument() + >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in, + ... CESnumber=0, projection='healpix', mapping_perpair=True) + >>> m = OutputSkyMapIGQU(projection=tod.projection, + ... nside=tod.nside_out, obspix=tod.obspix) + >>> for pair in tod.pair_list: + ... d = np.array([tod.map2tod(det) for det in pair]) + ... tod.tod2map(d, m, gdeprojection=True) + """ + if frequency_channel == 1: + pol_angs = self.pol_angs + elif frequency_channel == 2: + pol_angs = self.pol_angs2 + nbolofp = waferts.shape[0] npixfp = nbolofp / 2 nt = int(waferts.shape[-1]) @@ -691,14 +924,14 @@ def tod2map(self, waferts, output_maps, gdeprojection=False): assert npixfp == self.point_matrix.shape[0], msg assert nt == self.point_matrix.shape[1], msg - assert npixfp == self.pol_angs.shape[0], msg - assert nt == self.pol_angs.shape[1], msg + assert npixfp == pol_angs.shape[0], msg + assert nt == pol_angs.shape[1], msg assert npixfp == self.diff_weight.shape[0], msg assert npixfp == self.sum_weight.shape[0], msg point_matrix = self.point_matrix.flatten() - pol_angs = self.pol_angs.flatten() + pol_angs = pol_angs.flatten() waferts = waferts.flatten() diff_weight = self.diff_weight.flatten() sum_weight = self.sum_weight.flatten() @@ -739,7 +972,8 @@ def __init__(self, hardware, scanning_strategy, HealpixFitsMap, nside_out=None, pixel_size=None, width=140., cut_pixels_outside=True, array_noise_level=None, array_noise_seed=487587, - mapping_perpair=False, verbose=False): + array_noise_level2=None, array_noise_seed2=56736, + mapping_perpair=False, mode='standard', verbose=False): """ C'est parti! @@ -782,6 +1016,11 @@ def __init__(self, hardware, scanning_strategy, HealpixFitsMap, If True, assume that you want to process pairs of bolometers one-by-one, that is pairs are uncorrelated. Default is False (and should be False unless you know what you are doing). + mode : string, optional + Choose between `standard` (1 frequency band) and `dichroic` + (2 frequency bands). If `dichroic` is chosen, make sure your + hardware can handle it (see instrument.py) and HealpixFitsMap + should contain the inputs maps at different frequency. Examples ---------- @@ -794,6 +1033,22 @@ def __init__(self, hardware, scanning_strategy, HealpixFitsMap, >>> m = OutputSkyMap(projection=tod.projection, ... nside=tod.nside_out, obspix=tod.obspix, demodulation=True) >>> tod.tod2map(d, m) + + Same with dichroic detectors + >>> inst, scan, sky_in = load_fake_instrument(fwhm_in2=1.8) + >>> tod = TimeOrderedDataDemod(inst, scan, sky_in, mode='dichroic', + ... CESnumber=0, projection='healpix', mapping_perpair=True) + >>> d = np.array([tod.map2tod(det) for det in range(2)]) + >>> d1 = d[:, 0] + >>> d2 = d[:, 1] + >>> d1 = tod.demodulate_timestreams(d1) + >>> d2 = tod.demodulate_timestreams(d2) + >>> m1 = OutputSkyMap(projection=tod.projection, + ... nside=tod.nside_out, obspix=tod.obspix, demodulation=True) + >>> m2 = OutputSkyMap(projection=tod.projection, + ... nside=tod.nside_out, obspix=tod.obspix, demodulation=True) + >>> tod.tod2map(d1, m1, frequency_channel=1) + >>> tod.tod2map(d2, m2, frequency_channel=2) """ TimeOrderedDataPairDiff.__init__( self, hardware, scanning_strategy, HealpixFitsMap, @@ -802,7 +1057,10 @@ def __init__(self, hardware, scanning_strategy, HealpixFitsMap, cut_pixels_outside=cut_pixels_outside, array_noise_level=array_noise_level, array_noise_seed=array_noise_seed, + array_noise_level2=array_noise_level2, + array_noise_seed2=array_noise_seed2, mapping_perpair=mapping_perpair, + mode=mode, verbose=verbose) ## Prepare the demodulation of timestreams @@ -1951,14 +2209,14 @@ def set_goodpix(self): >>> m1.set_goodpix() """ ## We have 4 components: I, Q, U and G. - inonzero = [pix for pix in range(self.npixsky) if self.nhit[pix] > 4] + inonzero = [pix for pix in range(self.npixsky) if self.nhit[pix] > 10] self.goodpix = np.zeros((self.npixsky)) self.goodpix[inonzero] = 1 for ipix in inonzero: M = self.buildP(ipix) det = np.linalg.det(M) - if det == 0 or np.linalg.cond(M) > 1e8: + if det == 0 or np.linalg.cond(M) > 1e2: self.goodpix[ipix] = 0 def get_QU(self): @@ -2375,7 +2633,7 @@ def build_pointing_matrix(ra, dec, nside_in, nside_out=None, return index_global, index_local -def load_fake_instrument(nside=16, nsquid_per_mux=1, +def load_fake_instrument(nside=16, nsquid_per_mux=1, fwhm_in2=None, compute_derivatives=False): """ For test purposes. @@ -2402,7 +2660,7 @@ def load_fake_instrument(nside=16, nsquid_per_mux=1, ## Sky sky_in = HealpixFitsMap('s4cmb/data/test_data_set_lensedCls.dat', - do_pol=True, fwhm_in=0.0, + do_pol=True, fwhm_in=0.0, fwhm_in2=fwhm_in2, nside_in=nside, map_seed=48584937, compute_derivatives=compute_derivatives, verbose=False, no_ileak=False, no_quleak=False) @@ -2414,6 +2672,8 @@ def load_fake_instrument(nside=16, nsquid_per_mux=1, beam_seed=58347, projected_fp_size=3., pm_name='5params', type_hwp='CRHWP', freq_hwp=0.2, angle_hwp=0., verbose=False) + if fwhm_in2 is not None: + inst.make_dichroic(fwhm=fwhm_in2) ## Scanning strategy scan = ScanningStrategy(nces=2, start_date='2013/1/1 00:00:00',