Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

few edits on conviqt.py #346

Closed
wants to merge 2 commits into from
Closed

few edits on conviqt.py #346

Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
3c25a30
few edits on conviqt.py
May 8, 2020
fab460f
avoiding loops over all observations in the convolution method
May 11, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
107 changes: 57 additions & 50 deletions src/toast/todmap/conviqt.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -507,47 +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')
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why these changes? black will revert them back immediately and they clutter the review. It would be helpful to format the source before committing.

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 )
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same comment about formatting.

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)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Instead of get_buffer(theta, phi, psi, ..) try get_buffer(theta, phi, psi - psipol, ...)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

sorry, not psipol but psitot to include the HWP

convolved_data = self.convolve(sky, beamI00, detector, pnt, det, verbose)
# Q-beam convolution
pnt = self.get_buffer(theta, phi, psi, det, verbose)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

psi -> psi - psipol

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)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

psi -> psi-psipol

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)
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

Expand All @@ -556,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`.
Expand All @@ -573,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):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This implementation will fail if the detector in question is not available in all observations. Please see how get_pointing avoids this issue.

""" 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
Expand Down