Problem with generation of pointing in case of more than one detector

[paganol]

When I generate pointing for more than one detector, the array returned by get_pointings seems strange. For example, in the plot, in blue the field [0,:,0] of lbs.scanning.get_pointings(...) in case of 1 detector (only Top in a pixel), in orange the same field in case of two detectors generated (Top and Bottom) for the same pixel.

The two curves should be overlapped. But two detectors case shows something strange, it seems that the same theta angle is repeated two consecutive times. Probably this is realated to how the output of get_pointings gets filled when more than one detector is considered.
Last info, I'm using the dipole branch, I didn't check if this problem is there also for the master.

[ziotom78]

Hi @paganol , might you provide here the simplest script that reproduces the issue?

[paganol]

import litebird_sim as lbs
import numpy as np
import matplotlib.pylab as plt

imo = lbs.Imo()
sim = lbs.Simulation(base_path="./tut01",
                     start_time=0,
                     duration_s=1.5,
                     )

mft_file = sim.imo.query("/releases/v1.0/satellite/MFT/instrument_info")
sim.generate_spin2ecl_quaternions(imo_url="/releases/v1.0/satellite/scanning_parameters/")
instr = lbs.InstrumentInfo(name="MFT", 
    boresight_rotangle_rad=np.deg2rad(mft_file.metadata['boresight_rotangle_deg']),
    spin_boresight_angle_rad=np.deg2rad(mft_file.metadata['spin_boresight_angle_deg']),
    spin_rotangle_rad=np.deg2rad(mft_file.metadata['spin_rotangle_deg']),)
channel = lbs.FreqChannelInfo.from_imo(url='/releases/v1.0/satellite/MFT/M1-140/channel_info',imo=imo)

dets=[]
detquats=[]
detnames=[]
#Boresight
for detname in channel.detector_names[122:123]:
    det=lbs.DetectorInfo.from_imo(url="/releases/v1.0/satellite/MFT/M1-140/"+detname+"/detector_info",imo=imo)
    dets.append(det)
    detquats.append(det.quat)
    detnames.append({"name":det.name})


obs, = sim.create_observations(detectors=dets)

pointings1det = lbs.scanning.get_pointings(obs,
                                       spin2ecliptic_quats=sim.spin2ecliptic_quats,
                                       detector_quats=detquats,
                                       bore2spin_quat=instr.bore2spin_quat,
                                      )

dets=[]
detquats=[]
detnames=[]
#two detectors
for detname in channel.detector_names[122:124]:
    det=lbs.DetectorInfo.from_imo(url="/releases/v1.0/satellite/MFT/M1-140/"+detname+"/detector_info",imo=imo)
    dets.append(det)
    detquats.append(det.quat)
    detnames.append({"name":det.name})


obs, = sim.create_observations(detectors=dets)

pointings2det = lbs.scanning.get_pointings(obs,
                                       spin2ecliptic_quats=sim.spin2ecliptic_quats,
                                       detector_quats=detquats,
                                       bore2spin_quat=instr.bore2spin_quat,
                                      )

plt.figure()
plt.plot(pointings1det[0,0:31,0],lw=4)
plt.plot(pointings2det[0,0:31,0])

plt.figure()
plt.plot(pointings1det[0,0:31,0],lw=4)
plt.plot(pointings2det[0,0:31:2,0])
plt.plot(np.arange(14)+14,pointings2det[1,0:32:2,0])

[ziotom78]

Hmmm, it seems that the problem lies in the IMO itself. Lines 156086ff show that the quaternion associated with M03_030_QA_140T is (0, 0, 0, 1):

"metadata": {
  "name": "M03_030_QA_140T",
  "wafer": "M03",
  "pixel": 30,
  "pixtype": "MP1",
  "channel": "M1-140",
  "fwhm_arcmin": 30.8,
  "ellipticity": 0.0,
  "bandcenter_ghz": 140.0,
  "bandwidth_ghz": 42.0,
  "sampling_rate_hz": 19.0,
  "net_ukrts": 54.0,
  "pol_sensitivity_ukarcmin": 122.13653079313542,
  "fknee_mhz": 20.0,
  "fmin_hz": 0.01,
  "alpha": 1.0,
  "pol": "T",
  "orient": "Q",
  "quat": [
    0.0,
    0.0,
    0.0,
    1.0
  ]
}

However, lines 156053ff show that the quaternion for M03_030_QA_140B is (0, 0, 1/√2, 1/√2):

  "metadata": {
    "name": "M03_030_QA_140B",
    "wafer": "M03",
    "pixel": 30,
    "pixtype": "MP1",
    "channel": "M1-140",
    "fwhm_arcmin": 30.8,
    "ellipticity": 0.0,
    "bandcenter_ghz": 140.0,
    "bandwidth_ghz": 42.0,
    "sampling_rate_hz": 19.0,
    "net_ukrts": 54.0,
    "pol_sensitivity_ukarcmin": 122.13653079313542,
    "fknee_mhz": 20.0,
    "fmin_hz": 0.01,
    "alpha": 1.0,
    "pol": "B",
    "orient": "Q",
    "quat": [
      0.0,
      0.0,
      0.7071067811865475,
      0.7071067811865476
    ]
  },
  "quantity": "3a72c257-d447-4f10-a2e1-56c69ee1f90e",
  "spec_version": "1.0",
  "dependencies": []
}

[ziotom78]

No, sorry, I was wrong, litebird_sim encodes quaternions putting the scalar at the end, so it's perfectly fine. Now I see that your script points out to the possible culprit: a wrong mix of some matrix when using more than one detector. I'll look into this.

[ziotom78]

I believe I found the problem. The original purpose of Observation.get_pointings was to compute pointings for just one detector, as under the hood it uses all_compute_pointing_and_polangle, which assumes that there is only one detector.

However, in commit fa3cd39 the interface was changed to support multiple detectors, yet all_compute_pointing_and_polangle was left unchanged, apart from a call to numpy.reshape at the beginning of the function body. This change however had the only purpose to make the code still work with just one detector: its purpose was simply to verify that keyword parameters like detector_quat=… had become detector_quats=[…], but the … was always meant to be just one quaternion (as the tests in test/test_scanning.py show: there are no tests that verify what happens with more than one detector).

I'll try to figure out what is the best way to fix this.