diff --git a/src/toast/todmap/conviqt.py b/src/toast/todmap/conviqt.py index b2fb4eb7d..ec3f46b57 100644 --- a/src/toast/todmap/conviqt.py +++ b/src/toast/todmap/conviqt.py @@ -507,49 +507,39 @@ def exec(self, data): beam_file = self._beam_file[det] except TypeError: beam_file = self._beam_file.replace("DETECTOR", det) - beam_file_i00 = beam_file.replace(".fits", "_I000.fits") - beam_file_0i0 = beam_file.replace(".fits", "_0I00.fits") - beam_file_00i = beam_file.replace(".fits", "_00I0.fits") + beam_file_i00= beam_file.replace('.fits', '_I000.fits') + beam_file_0i0= beam_file.replace('.fits', '_0I00.fits') + beam_file_00i= beam_file.replace('.fits', '_00I0.fits') - beamI00 = self.get_beam(beam_file_i00, det, verbose) - beam0I0 = self.get_beam(beam_file_0i0, det, verbose) - beam00I = self.get_beam(beam_file_00i, det, verbose) + + beamI00 =self.get_beam(beam_file_i00, det,verbose ) + beam0I0 =self.get_beam(beam_file_0i0, det,verbose ) + beam00I =self.get_beam(beam_file_00i, det,verbose ) detector = self.get_detector(det) theta, phi, psi = self.get_pointing(data, det, verbose) + psipol = self._get_psipol(data, det) + + hwpang = self._get_hwpangle(data) + if hwpang is None: + psitot = psipol + else: + psitot = psipol + 2 * hwpang + + # I-beam convolution + pnt = self.get_buffer(theta, phi, psi , det, verbose) + convolved_data = self.convolve(sky, beamI00, detector, pnt, det, verbose) + # Q-beam convolution + pnt = self.get_buffer(theta, phi, psi, det, verbose) + convolved_data += np.cos(2*psitot ) * self.convolve(sky, beam0I0, detector, pnt, det, verbose) + # U-beam convolution pnt = self.get_buffer(theta, phi, psi, det, verbose) + convolved_data += np.sin(2* psitot ) * self.convolve(sky, beam00I, detector, pnt, det, verbose) del theta, phi, psi - for obs in data.obs: - tod = obs["tod"] - focalplane = obs["focalplane"] - psipol = self._get_psipol(focalplane, det) - # import pdb - # pdb.set_trace() - hwpang = self._get_hwpangle(tod) - if hwpang is None: - psitot = 2 * psipol - else: - psitot = 2 * psipol + 4 * hwpang - - convolved_data = self.convolve( - sky, beamI00, detector, pnt, det, verbose - ) - theta, phi, psi = self.get_pointing(data, det, verbose) - pnt = self.get_buffer(theta, phi, psi, det, verbose) - del theta, phi, psi - convolved_data += np.cos(psitot) * self.convolve( - sky, beam0I0, detector, pnt, det, verbose - ) - theta, phi, psi = self.get_pointing(data, det, verbose) - pnt = self.get_buffer(theta, phi, psi, det, verbose) - del theta, phi, psi - convolved_data += np.sin(psitot) * self.convolve( - sky, beam00I, detector, pnt, det, verbose - ) - - self.calibrate(data, det, beamI00, convolved_data, verbose) - - self.cache(data, det, convolved_data, verbose) + + self.calibrate(data, det, beamI00, convolved_data, verbose) + + self.cache(data, det, convolved_data, verbose) del pnt, detector, beamI00, beam0I0, beam00I, sky @@ -558,6 +548,7 @@ def exec(self, data): return + def _get_detectors(self, data): """ Assemble a list of detectors across all processes and observations in `self._comm`. @@ -575,31 +566,45 @@ def _get_detectors(self, data): all_dets = self._comm.bcast(dets, root=0) return all_dets - def _get_hwpangle(self, tod): + + def _get_hwpangle (self, data ): """ Return the spinning HWP angle and the multiplicative factor to be applied when added to the signal polarization angle. """ hwpang = None - try: + all_hwpang =[] + for obs in data.obs: + tod= obs["tod"] + nsamp = tod.local_samples[1] hwpang = tod.local_hwp_angle() - except: - pass - return hwpang + if hwpang is not None : + all_hwpang .append(hwpang) + else : + return hwpang - def _get_psipol(self, focalplane, det): + return np.hstack (all_hwpang ) + + + def _get_psipol(self, data, det): """ Parse polarization angle in radians from the focalplane dictionary. """ - if det not in focalplane: - raise RuntimeError("focalplane does not include {}".format(det)) - if "pol_angle_deg" in focalplane[det]: - psipol = np.radians(focalplane[det]["pol_angle_deg"]) - elif "pol_angle_rad" in focalplane[det]: - psipol = focalplane[det]["pol_angle_rad"] - else: - raise RuntimeError("focalplane[{}] does not include psi".format(det)) - return psipol + all_psipol = [] + for obs in data.obs: + nsamp = obs["tod"].local_samples[1] + focalplane = obs["focalplane"] + if det not in focalplane: + raise RuntimeError("focalplane does not include {}".format(det)) + if "pol_angle_deg" in focalplane[det]: + psipol = np.radians(focalplane[det]["pol_angle_deg"]) + elif "pol_angle_rad" in focalplane[det]: + psipol = focalplane[det]["pol_angle_rad"] + else: + raise RuntimeError("focalplane[{}] does not include psi".format(det)) + + all_psipol .append( np.ones(nsamp) *psipol ) + return np.hstack(all_psipol ) def _get_epsilon(self, focalplane, det): """ Parse polarization leakage (epsilon) from the focalplane