# Data Reduction Process and Pipeline for the NIC Polarimetry Mode in Python,
NICpolpy

Yoonsoo P. Bach Astronomy Program, Department of Physics and Astronomy, Seoul
National University,
Gwanak-ro 1, Gwanak-gu, Seoul 08826, Republic of Korea  SNU Astronomy Research
Center, Department of Physics and Astronomy, Seoul National University,
Gwanak-ro 1, Gwanak-gu, Seoul 08826, Republic of Korea<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Masateru Ishiguro Astronomy Program,
Department of Physics and Astronomy, Seoul National University,
Gwanak-ro 1, Gwanak-gu, Seoul 08826, Republic of Korea  SNU Astronomy Research
Center, Department of Physics and Astronomy, Seoul National University,
Gwanak-ro 1, Gwanak-gu, Seoul 08826, Republic of Korea  Jun Takahashi Nishi-
Harima Astronomical Observatory, Center for Astronomy, University of Hyogo,
407–2 Nishigaichi, Sayo-cho, Hyogo 679–5313, Japan  Jooyeon Geem Astronomy
Program, Department of Physics and Astronomy, Seoul National University,
Gwanak-ro 1, Gwanak-gu, Seoul 08826, Republic of Korea  SNU Astronomy Research
Center, Department of Physics and Astronomy, Seoul National University,
Gwanak-ro 1, Gwanak-gu, Seoul 08826, Republic of Korea

(Received 2022 October 26; Accepted 2022 December 13)

###### Abstract

A systematic way of data reduction for the NIC polarimetry mode has been
devised and implemented to an open software called NICpolpy in the programming
language python (tested on version 3.8–3.10 as of writing). On top of the
classical methods, including vertical pattern removal, a new way of diagonal
pattern (Fourier pattern) removal has been implemented. Each image undergoes
four reduction steps, resulting in “level 1” to “level 4” products, as well as
nightly calibration frames. A simple tutorial and in-depth descriptions are
provided, as well as the descriptions of algorithms. The dome flat frames
(taken on UT 2020-06-03) were analyzed, and the pixel positions vulnerable to
flat error were found. Using the dark and flat frames, the detector
parameters, gain factor (the conversion factor), and readout noise are also
updated. We found gain factor and readout noise are likely constants over
pixel or “quadrant”.

methods: data analysis — methods: observational — techniques: image processing
— techniques: polarimetric

##  1 Introduction

In this work, we describe a way to reduce Nishi-Harima Astronomical
Observatory (NHAO) Nishiharima Infrared Camera (NIC) polarimetric mode data.
The process is implemented to NICpolpy111 The stable version is registered to
the Python Package Index (PyPI): https://pypi.org/project/NICpolpy/. The
development version is available via GitHub at
https://github.com/ysBach/NICpolpy. See Sect. 3 for details. , and the
implementation details are given. The term reduction in this work is
restricted to preprocessing of image data, including artifact removal and the
standard dark and flat corrections, but excluding photometry/polarimetry and
error analysis. We analyzed the dome flat frames (UT 2020-06-03) and examined
if any pixels have large uncertainty due to imperfect flat correction or hot
pixels. Also, the effect of the instrumental rotator and half-wave plate (HWP)
angles are shown.

During the development of NICpolpy, the detector parameters, viz., gain factor
(conversion factor, unit of electrons per ADU) and readout noise (unit of
electrons), were recalculated using dark frames (UT 2019-11-21) and the dome
flat frames. These had already been determined in Ishiguro et al. (2011) more
than a decade ago and should be updated due to a possible secular change.
Throughout this work, the three detectors associated with three filters,
$\mathrm{J}$, $\mathrm{H}$, and $\mathrm{K_{s}}$, are referred to as J-, H-,
and K-band, respectively.

In Sect. 2, the raw image and reduction processes implemented to NICpolpy
package are described. Then the practical usage of the package is shown in
Sect. 3, with additional descriptions of algorithms and implementations. Sect.
4 describes an analysis of polarimetric dome flat images to find any locations
vulnerable to flat fielding uncertainties. Finally, in Sect. 5, a new
estimation of instrument parameters (gain and readout noise), as well as the
corresponding statistical analyses, are described.

##  2 Description of Image and Data Reduction

In NICpolpy, the original image is called level 0 or lv0. Unfortunately, level
0 images are saved in 32-bit integer format (BITPIX=32), unnecessarily
doubling the file size (a file is 4.2 MB instead of 2.1 MB). This is because
NIC does not produce bias frames, so the pixel values can be negative, while
the standard FITS format does not support a signed 16-bit integer. One
possibility is to add an arbitrary constant bias (e.g., few thousands) and
save as an unsigned 16-bit integer ($0\mathrm{-}32,767$), because anyway the
maximum signal for NIC is $<2\times 10^{4}\,\mathrm{ADU}$ (Ishiguro et al.,
2011). Another possibility is to use BZERO key in the FITS header. At this
moment, neither is implemented to NICpolpy because it is something to be
implemented when the data are saved in the first place.

After the column-pattern subtraction, it is now called level 1 or lv1. The
columnar pattern is called vertical pattern in NICpolpy. Then an additional
semi-periodic noise, the so-called Fourier pattern, across the detector frame
remains. This is completely random for each frame (although the exact cause is
unsolved) and has an amplitude of a few analog-to-digital units (ADU). After
this pattern is removed, the image is called level 2 or lv2. Then, usual dark
subtraction, flat fielding, and bad pixel interpolation results in level 3 or
lv3. Finally, after an optional fringe subtraction and cosmic-ray removal, the
result is called level 4 or lv4.

Fig. 1 shows a sample FITS frame during the reduction steps. The raw frame
consist of four parts of size $512\times 512\,\mathrm{pix}$, and we call them
quadrants. As is visible, the light from the exposure falls onto a small
portion in the right half of the full array, indicated as green boxes in Fig.
1b. The dark green part is called the lit part. The exact locations of these
dark green boxes are summarized in Table 1. In this format, each lit part has
size of $150\times 430\,\mathrm{pix}$. Hence, any region well outside these
boxes, including B1, B1’, B2, B2’, FFTo, and FFTe, are dark overscan areas.
After processing until level 4, it is recommended to crop only the central
region of the lit part, i.e., the light green boxes, called the region of
interest (RoI). This is because the edge pixels are affected by the imperfect
flat field and can have arbitrarily large or small pixel values, and can
affect the sky estimation in photometric/polarimetric measurements. We
normally crop out 10 pixels around the edge (so the final analyzed region is
$130\times 410\,\mathrm{pix}$).

part | FITS section
---|---
J-band, o-ray | [534:683, 306:735]
J-band, e-ray | [719:868, 306:735]
H-band, o-ray | [564:713, 331:760]
H-band, e-ray | [741:890, 331:760]
K-band, o-ray | [560:709, 341:770]
K-band, e-ray | [729:878, 341:770]
Table 1: The FITS sections of each lit part. The FITS section is in the format
of $[x_{\mathrm{min}}:x_{\mathrm{max}},y_{\mathrm{min}}:y_{\mathrm{max}}]$ in
1-indexing (first and last pixels inclusive). Figure 1: An example exposure
(2019-10-22 UT, Vesta). (a) The raw data, (b) schematic diagram describing the
full NIC frame in polarimetry mode, (c) level 1, (d) level 2, (e) level 3, and
(f) level 4. All images are in the same limits and scale. Numbers in (a) and
(b) are in the unit of a pixel, showing the dimensions of indicated areas.
Light comes only to some part on the right half of the full array, as in (a)
called the lit part. These are the green-boxed regions in (b). The regions
indicated by letter B in (a) and (b) are used to infer the vertical pattern,
resulting in (c). The left side of the array, indicated by FFTo and FFTe in
(a) and (b), is used to infer the Fourier pattern, resulting in (d). Figures
(c) and (d) are scaled identically such that the Fourier pattern along the
diagonal direction is exaggerated. Note that the pattern is reduced in the
FFTo/FFTe areas and lit parts in (d), compared to (c). After splitting
o-/e-ray regions, standard flat fielding, and bad-pixel interpolation, level 3
is obtained (e). Finally, a conventional cosmic-ray rejection technique is
applied to obtain level 4 (f). In this example, fringe removal was turned off
and no cosmic ray was detected, so the images in (e) and (f) are identical.
The dark green areas in (b) indicate the 10-pixel bezels that are ignored
(cropped) in any data analyses after level 4. See text for details.

In this section, the detailed process to obtain the reduced image (levels 1-4)
is described. B1, B1’, B2, and B2’ regions are used to make level 1 data,
while FFTo and FFTe regions are used for level 2 data. The FFTo and FFTe
regions have the same size ($150\times 1024\,\mathrm{pix}$), and are set such
they overlap with the lit part when left half is shifted by
$+512\,\mathrm{pix}$ to the right.

###  2.1 Level 1 and Initial Bad-Pixel Mask

Since NIC data contains random and additive vertical pattern (Ishiguro et al.,
2011), locations B1, B1’, B2, and B2’ were used to remove this vertical
pattern for the upper-right, upper-left, lower-right, and lower-left
quadrants, respectively. For each column within each of these boxes, a 2-sigma
clipped median is calculated. This value is subtracted from all the pixels in
the same column in the same quadrant. Then level 1 image is obtained.

On UT 2019-11-21, 20 dark frames with an exposure time of 180 s were obtained
for all three detectors (J, H, and K bands). All these images were first
reduced to level 1. Then after sigma-clipped median combining these dark
frames for each detector, any pixel with values larger than
$\mathrm{max}\\{5,D_{\mathrm{med}}+5\sigma_{D}\\}\,\mathrm{ADU}$ or smaller
than $\mathrm{min}\\{-5,D_{\mathrm{med}}-5\sigma_{D}\\}\,\mathrm{ADU}$ are
flagged as bad-pixel for each detector. Here, $D_{\mathrm{med}}$ and
$\sigma_{D}$ are the 3-sigma clipped median and standard deviation of the
combined dark frame, respectively. Secondly, the sigma-clipped sample standard
deviation was calculated for each pixel position, and any pixel with high
variance was flagged.

On UT 2020-06-03, 640 dome flat frames with 2 s exposure time were obtained
(with varying instrumental rotation angle INSROT and half-wave plate rotator
angle POL-AGL1). All frames are reduced to level 1. Then similar to the dark
frames, the frames were median combined, and the sample standard deviation
maps were generated. Each of the o-/e-ray regions was normalized separately
before the combination because the average raw flat count is different in the
two regions by up to $\sim 10\%$. Then the pixels with too low value or too
high variance are flagged.

Combining all the flagged pixels from these dark and flat frames, initial bad-
pixel masks were generated for each detector (total 3 mask files). These are
called the imask in NICpolpy, and shown in Fig. 2. As a result, $\sim 500$
pixels are masked for each light-green RoI in Fig. 1. Since the total size
after cropping the 10-pixel bezel is $130\times 410\sim 50,000\,\mathrm{pix}$,
$\sim 1\%$ of the region is masked. It is found that the position we
traditionally put our science target is located where the least number of
pixels are masked.

Figure 2: The distribution of bad pixels in imask frames for J, H, and K
bands. The right panels show zoomed-in versions of the red box areas in the
left panels, i.e., around the lit parts (green boxes) in Fig. 1.

###  2.2 Nightly Master Mask and Level 2

Even after removing the vertical pattern, it is known that NIC has an
additional quasi-periodic artifact along the diagonal direction (for example,
Fig. 28 of Ishiguro et al. (2011) and Fig 5 of Takahashi (2018)). Although the
specific reason is still unknown, this is likely an additive noise due to an
imperfectly shielded external electric signal. After thorough investigations,
we could not find a way to remove the noise by tiling columns or rows of the
full frame or each quadrant into the 1-dimensional array, and the Fourier
transform. Accordingly, this noise is unlikely a perfect periodic wave.
Fortunately, it was found that this pattern is quite similar in the left/right
quadrants. That is, a direct subtraction of the left part from the right part
reduces the pattern. However, this simple subtraction is vulnerable to any
cosmic ray hit, amp-glow, newly developed hot-pixels, or even just standard
noise (e.g., read-out noise) in the left half of the array. In addition,
standard cosmic-ray rejection algorithms (e.g., van Dokkum (2001)) is
difficult to be applied directly, because (i) it takes a long time and (ii)
the fine-tuning of the parameters for dark (left) area was difficult.
Therefore, we applied the following strategy.

As described in Sect. 2.1, we already made the initial pixel mask (imask) from
dark frames of UT 2019-11-21. Whenever dark frames are available, those of the
longest exposure are used to make a similar pixel mask (called dmask). Then
the nightly master mask, mmask, is made by combining those: mmask =
imask|dmask (bitwise OR operation). If there are no dark frames at the nith,
mmask = imask. This process can remove newly developed bad pixels after UT
2019-11-21222 It is also possible to use the dark regions (outside the lit
parts) of science exposure frames to find additional bad pixels. This might be
useful when there are no nightly dark frames. However, it is less important
because the bad pixels in the RoI cannot be retrieved. Hence, it is not
implemented yet. . We use dmask rather than just subtracting dark, because the
absolute value of hot pixels often varies from frame to frame.

Afterwards, similar to, but slightly different from, the previous work
Takahashi et al. (2021), the left half of the $1024\times 1024\,\mathrm{pix}$
frame is used to estimate and remove the periodic pattern (also called the
Fourier pattern in NICpolpy). To minimize the computation time, this pattern
inference is made on FFTo/FFTe regions only rather than the full left-half
frame. The pixels within FFTo/FFTe regions that are flagged in mmask are
interpolated. This interpolation is similar to IRAF.PROTO.FIXPIX, prioritizing
the vertical direction over the horizontal one (the priority can be modified
by the user), implemented for $N$-dimensional array by using scipy.

In addition, the pixels that are likely cosmic-ray, newly developed bad-pixel,
or high-noise artifacts are replaced by the median-filtered value. Three
parameters are calculated from the image in this algorithm, called the Median-
filter Bad-Pixel Masking (MBPM):

$C_{1}=\mathcal{I}-\mathcal{M}~{},\\\ C_{2}=\mathcal{I}/\mathcal{M}~{},\\\
C_{3}=(\mathcal{I}-\mathcal{M})/\sigma(\mathcal{I})~{}.$ (1)

Here, $\mathcal{I}$ is the original image, $\mathcal{M}$ is the median
filtered image, and $\sigma(\mathcal{I})$ is the sigma-clipped standard
deviation at certain user-tunable region(s) of $\mathcal{I}$. In this work,
the pixels with $|C_{1}|>5\,[\mathrm{ADU}]$ or $C_{2}\notin[0.5,2]$ or
$|C_{3}|>5$ are flagged and replaced with the median-filtered value. The
regions for $\sigma(\mathcal{I})$ is [50:450, 100:900] by default in the FITS
section format. It is not iterated multiple times. The median-filter size was
$5\times 5\,\mathrm{pix}$, and a very small real number is added to
$\mathcal{M}$ to avoid zero-division. This MBPM algorithm is orders of
magnitude quicker than traditional cosmic-ray rejection algorithms (e.g., van
Dokkum (2001)) while removing artifacts/hot-pixels/cosmic-rays reasonably well
and is easily applicable to the dark region without fine-tuning of the
parameters 333 MBPM is not intended to detect cosmic rays in the presence of
celestial objects. Unlike algorithms that use edge detection and convolutions
to find cosmic rays, MBPM only uses the deviation from the median filtered
image to find outliers without any additional convolution or edge detection.
Thus, MBPM may easily flag stellar objects. MBPM is intended to quickly find
bad pixels or cosmic rays in dark frames or sky regions where no point source
is expected. MBPM is separately available as a function at
NICpolpy.ysfitsutilpy4nicpolpy.medfilt_bpm. .

After smoothing out these artifacts, a fast Fourier transform (FFT) is applied
to each column in the FFTo and FFTe regions. Then the components with
wavelength $>100\,\mathrm{pix}$ are extracted, so that two pattern maps with
sizes $150\times 1024\,\mathrm{pix}$ are made. These pattern maps are
subtracted from the right part, overlapping with o-ray and e-ray RoI,
respectively. This image undergoes an additional vertical pattern subtraction,
as in level 1, to remove any remaining vertical pattern (likely smaller than
$0.5\,\mathrm{ADU}$).

The resulting image is denoted as level 2. So far, it is found that this
process reduces the periodic pattern, especially for the frames having larger
contamination ($\lesssim\pm 5$ ADU, see, Fig. 1c). We confirmed that this
process does not change the final Stokes’ parameter values while greatly
reducing background fluctuation. These reduction processes would be beneficial
for extended source studies, e.g., the Moon, earth shine (Takahashi et al.,
2021) or comet polarimetry.

###  2.3 Master Flat and Level 3 and 4

We found no clear and severe difference in the combined flats for different
POL-AGL1. We hence did not split the master flat based on the rotator angles
(Sect. 4 and supporting material in Appendix B). The master flat was generated
by the following steps:

  1. 1.

sum combines each polarimetric set (four frames at four POL-AGL1),

  2. 2.

normalize each of the sum images based on the sigma-clipped median in the
o-ray RoI,

  3. 3.

median combines these by sigma-clipping.

Because the flat counts are not normalized separately for o-/e-ray regions,
the final flat has slightly different values for the o-/e-ray regions ($\sim
1$ and $\sim 1.1$, respectively). This difference is not important because NIC
samples the flux twice per each Stokes’ parameter (e.g., POL-AGL1 of 0 and
$45^{\circ}$ for $q$) to cancel out this offset effect.

After removing the instrumental noises (vertical and Fourier patterns), a
standard preprocessing, including dark subtraction and flat fielding, was
done. Using mmask, some pixels are interpolated. Then each frame is split into
two frames, o-ray and e-ray (each $150\times 430\,\mathrm{pix}$). This is
level 3 data. If the user demands, blank sky frames are extracted and combined
after scaling to make a master fringe frame. Otherwise, blank sky frames will
be discarded at this stage.

Finally, fringe subtraction and cosmic-ray rejection are done on each image,
saving the final result as level 4. Fringe is approximated as an additive
feature in imaging observation. In NICpolpy, fringe correction is implemented
either by scaling the normalized sky fringe or direct subtraction. The latter
is sometimes erroneous because the sky brightness (and even the fringe
pattern) can change rapidly. As it will be justified later in future work, the
fringe subtraction is implemented but omitted in real data reduction. Each
fringe frame is divided by the sigma-clipped average pixel value to eliminate
the overall sky brightness change. After this normalization, frames are
combined to make a master fringe frame (for the fine tuning, see Sect.
3.4.10). Since the fringe pattern is well known to be variable over time,
NICpolpy automatically combines fringe based on time stamps, too (Sect.
3.4.9). The HWP angle, POL-AGL1, is ignored when combining master fringe
because we found that separating fringe by HWP angle does not improve the
final results. For the cosmic-ray rejection, astroscrappy (McCully et al.,
2018), which implements L.A.Cosmic (van Dokkum, 2001), is used (version
1.1.1). After the investigation, the parameters similar to L.A.Cosmic are
found to be able to remove cosmic rays efficiently.

The blank sky (fringe) frame can approximately be regarded as a dark frame,
especially for short exposures and when the sky is dark. This, however, is
vulnerable to the change in sky level, spatial polarization patterns in the
sky, and typical pixel noise (readout noise and Poisson noise) in sky frames
because we have few sky frames. Therefore, this function has not yet been
implemented to NICpolpy. Another approach is to conduct a dithering
observation to construct a fringe map. This technique is also difficult for
NIC because of the small field-of-view, especially when sky level and/or
fringe pattern vary quickly relative to the exposure time. Finally, it can be
shown that the fringe pattern has a negligible effect on Stokes’ parameter
determination for beam-splitter type aperture polarimetry when the effect of
fringe is a few-percent level of the aperture sum and the true polarization
degree of the target is $\lesssim\mathrm{few}\%$ (Bach et al. in prep.). Thus,
we ignore fringe in this work.

##  3 Code Usage and Implementation

This section describes the practical usage of the package. Since NICpolpy is
registered to the Python Package Index (PyPI), the installation is invoked
once typing

$ pip install NICpolpy (2)

on the terminal444 The terminal must have pip installed:
https://pypi.org/project/pip/ . The development version is available via
GitHub555 https://github.com/ysBach/NICpolpy . After the installation,
calibration frames imask and flats are required. They are available in a
separate GitHub repository (Appendix B).

First, the simplest way to use NICpolpy is briefly discussed. Then the output
FITS files and log files are explained. Finally, the core algorithmic
implementations of reduction processes (Sect. 2) are described step by step,
as well as intermediate output files. Although the detailed usage can include
many more arguments (and may change in future developments), here we describe
the core ideas and basic usage that will suffice most of the users’ needs.

###  3.1 Basic Usage

The most recent development of NICpolpy as of writing is happening on Python
3.10 with macOS 12. However, since each part of the code was written OS-
agnostic, it will work without any problem on any recent major OS (Windows,
macOS, Linux) and python versions (3.8 or later).

First, import packages and initialize the reducer by giving the names of the
directories where the original data and calibration data are stored666 This
example is also provided in Appendix B. :

    # import nicpolpy package
    import nicpolpy as nic

    # Do not print useless warnings
    import warnings
    from astropy.utils.exceptions import AstropyWarning
    warnings.filterwarnings(’ignore’,
        append=True, category=AstropyWarning)

    # Cell 0: Initialize the reducer
    npr = nic.NICPolReduc(
        name="SP_20190417",
        inputs="_original_32bit/190417/raw/*.fits",
        mflats="cal-flat_20180507-lv1/*.fits",
        imasks="masks/*.fits",
        verbose=1
    )

For the reducer, the variable name is the prefix that will be used for result
files (see Sect. 3.2 and 3.3). The argument inputs, mflats and imasks
indicates the regex-style for input FITS files, master flat files, and imask
files, relative to the current working directory, respectively. The next steps
are

      # Cell 1: Planner for master dark combination
      #   5 sec
      _ = npr.plan_mdark()
      # Cell 2: Combine dark & planner for master mask
      #   11 sec
      _ = npr.comb_mdark_dmask()
      _ = npr.plan_mmask()
      # Cell 3: Make master mask
      #   << 1 sec
      _ = npr.comb_mmask()
      # Cell 4: Make lv1 & planner for lv2
      #   35 sec/456 FITS
      _ = npr.proc_lv1()
      _ = npr.plan_lv2()
      # Cell 5: Make lv2 & planner for lv3
      #   4 min / 357 FITS
      _ = npr.proc_lv2()
      _ = npr.plan_lv3()
      # Cell 6: Make lv3
      #   2 min
      _ = npr.proc_lv3()
      # Cell 7: Planner for lv4 of fringe frames
      #   1 min
      _ = npr.plan_ifrin()
      # Cell 8: Planner for master fringe
      #   35 sec
      _ = npr.plan_mfrin()
      # Cell 9: Make master fringe
      #   few sec
      _ = npr.comb_mfrin()
      # Cell 10: Planner for lv4
      #   few sec
      _ = npr.plan_lv4(add_mfrin=False)
      # Cell 11: Make lv4
      #   2 min
      _ = npr.proc_lv4()

The code lines indicated by different Cell numbers are recommended to run
after checking the results from the previous Cell. Note that, except for the
input file directories, the user does not have to specify any specific
parameters and quickly obtain the final level 4 files. For specific arguments
for fine-tuning, see Sect. 3.4. The comments with xx min or xx sec are the
time spent on Mac Book Pro with m1Pro chip (8 performance, 2 efficiency cores
with 16 GB RAM).

From now on, the following simplifying notations are used:

$\displaystyle\textrm{log directory}=\texttt{\\_\\_logs/<name>}~{},$
$\displaystyle\textrm{level~{}}\texttt{<i>}\mathrm{~{}directory}=\texttt{\\_lv<i>/<name>}~{},$
$\displaystyle\texttt{summ(X)}=\texttt{<name>\\_summ\\_X.csv}~{},$
$\displaystyle\texttt{plan(X)}=\texttt{<name>\\_plan-X.csv}~{},$

where name is the argument given in Cell 0, <i> is an integer from 1 to 4, and
X is any name. For example, plan(MDARK) in this specific example means the
file SP_20190417_plan-MDARK.csv, inside the corresponding log directory,
__logs/SP_20190417/. After running the code cell’s above, multiple outputs are
generated. Most importantly, the FITS files are saved in the corresponding
level directory, _lv<i>/<name>. Also, the log information and all the
calibration files are saved in the log directory, __logs/<name>. The contents
are described in the following subsections.

###  3.2 Output: FITS Files and Level Directories

Any FITS file written by NICpolpy has the following naming convention:

      <FILTER>_<System YYYYMMDD>
      _<COUNTER:04d>_<OBJECT>_<EXPTIME:.1f>
      _<POL-AGL1:04.1f>_<INSROT:+04.0f>
      -PROC-{PROCESS}_{QU}{SETID:03d}[_{OERAY}].fits

Each part means

  * •

FILTER: Lower-cased filter (j, h, or k).

  * •

YYYYMMDD: The date when this image was obtained. Same as the original data
saved by NIC.

  * •

COUNTER: The image counter in 4 digits with leading 0’s. It starts from 1.

  * •

OBJECT: The object name as in FITS header key OBJECT.

  * •

EXPTIME: The exposure time in seconds as in FITS header key EXPTIME.

  * •

POL-AGL1: The HWP rotator angle as in FITS header key POL-AGL1. One of 00.0,
45.0, 22.5 or 67.5. It is xxxx if there is no POL-AGL1 in the header.

  * •

INSROT: The instrument rotator angle at exposure as in FITS header key INSROT.

  * •

PROCESS: The history of processing. It can contain D (dark subtraction), F
(flat-fielding), T (trimmed), C (cosmic-ray rejection), R (fringe
subtraction), P (bad pixel interpolation), v (vertical pattern removal), and f
(Fourier pattern removal). The processing order is from left to right.

  * •

QU: Either q or u, which means whether the frame will be used for the raw $q$
or $u$ parameter calculation. Former is when POL-AGL1 is either 00.0 or 45.0,
and the latter is when either 22.5 or 67.5.

  * •

SETID: The 3-digit integer counter with leading 0’s for the group (set) id for
polarimetric measurement. One set consists of four exposures with different
POL-AGL1 angles, and the counter increases after the cycle. It starts from 1.

  * •

OERAY: Either o or e, indicating whether this is o-ray or e-ray RoI (only for
level 3 and 4).

The resulting files will be saved in the corresponding level directories. At
each level directory, a sub-directory for thumbnails of each FITS frame called
thumbs/, is generated when making the planner for the level <i+1> (see Sect.
3.4). There can be other sub-directories, _lv1/<name>/dark, if there are
nightly darks, and _lv4/<name>/frin, if there are fringe (blank sky) frames.

###  3.3 Output: Log Directory (__logs/)

After initializing the reducer (Cell 0), __logs/<name> directory is made. This
is called the log directory, and its tree structure is shown in Fig. 3. This
is made to archive virtually all necessary information to guarantee the
reproducibility of the user’s work. The calibration files (master flat and
mask files) are copied and saved, and all the thumbnail images related to
calibration will be saved here. Inside here are multiple CSV files and
directories:

  * •

cal-xxxxx/: Contains calibration FITS frames.

  * •

thumbs_xxxx/: The thumbnails with certain statistics of each frame for quick-
look.

  * •

summ(X): The summary of FITS frames for each level or master calibration
frames, e.g., SP_20190417_summ_lv0.csv.

  * •

plan(X): The planner before proceeding to each level or master calibration
frame generation, e.g., SP_20190417_plan-MDARK.csv.

Figure 3: The final tree structure of the log directory, __logs/SP_20190417,
as an example. Total 13 directories, 324 files, and 53.5 MB. The size of each
sub-directory is indicated as comments.

###  3.4 Implementation Details

Here, the details of implementation and the outputs from each Cell number are
described. As most of the tools do, NICpolpy also may undergo quick and
breaking changes when necessary. So for the most details, such as tuning
output thumbnail files, which do not affect the scientific outcomes, it is
recommended to look at the source code. Also, there are multiple arguments or
hard-coded parts that are useful merely for debugging purposes, which are not
described here.

####  3.4.1 Cell 0: Initialization

After the initialization, the log directory and level 1-4 directories are
automatically generated. Then inside the log directory, cal-imask/ and cal-
mflat/, based on imasks and mflats arguments, are made. These are not made if
arguments are not given (None). Also, the summary of level 0 files are made as
summ(lv0) inside the log directory. This contains some selected header
information from each FITS file. If imasks or mflats are given, corresponding
summary files are also made: summ(MFLAT) and summ(IMASK). These two contains
little information, and made only for the consistency.

####  3.4.2 Cell 1: Master Dark Planner

The important default arguments are:

    npr.plan_mdark(
        method="median",
        sigclip_kw=dict(sigma=2, maxiters=5),
        fitting_sections=None,
    )

Among level 0 files, it finds those with OBJECTS = "DARK" in the header. If
there is any nightly dark frames, they undergo vertical pattern removal
(pattern is estimated in the regions fitting_sections, and using the statistic
method, obtained by sigma-clipping with sigclip_kw), and saved into dark/
inside level 1 directory. The planner plan(MDARK) is made. This planner
contains file name, timestamp, filter, and exposure time information. There is
a column named REMOVEIT, and if the user fills here with any value other than
0, the frame will be ignored in the master dark combination. For each image, a
thumbnail is made with some statistical information and saved to
thumbs_mdark_plan/ in the log directory. An example is shown in Fig. 4. The
user may look at the thumbnails and pinpoint strange dark frames to fill in
REMOVEIT. For example, at one night, we could find multiple dark frames have
abnormally high median values of $>5\,\mathrm{ADU}$ only at K-band based on
these thumbnails. Such frames can be removed before making the master dark in
the next Cell.

Figure 4: An example of saved thumbnail images in
__logs/<name>/thumbs_mdark_plan. It is a DARK frame with 8 s exposure in
H-band, as shown in the title. The statistics are based on pixels ignoring
20-pixel edge bezels. The statistics are the number of pixels, minimum,
maximum, average, median, sample standard deviation, and median-absolute-
deviation (MAD) derived standard deviation of pixel values. pct are the
corresponding percentile pixel values specified by percentiles. zmin and zmax
are the famous zscale minimum/maximum values. ext_lo and ext_hi are the top
three pixel values in the lowest and highest regimes, respectively.

####  3.4.3 Cell 2: Master Dark and Master Mask Planner

The important default arguments are:

    npr.comb_mdark_dmask(
        combine_kw=dict(combine="med", reject="sc", sigma=3),
        dark_thresh=(-10, 50),
        dark_percentile=(0.01, 99.99),
    )
    npr.plan_mmask()

The master dark is then made by median combined with 3-sigma clipping (tunable
by combine_kw). The resulting master darks are saved in cal-mdark/, while the
thumbnails are saved in thumbs_mdark/ in the log directory. An example of
thumbnails is shown in Fig. 5. The summary of master darks is also saved as
summ(MDARK). For each master dark frame, any pixel that are outside
$[-10,50]\,\mathrm{ADU}$ range, or $[0.01,99.99]$ percentile range are flagged
(tunable by dark_thresh and dark_percentile, respectively).

The generated intermediate dmask files and their thumbnails are saved in cal-
dmask/ and thumbs_dmask/, and summary file as summ(DMASK) in the log
directory, respectively. Each dmask frame is in the format of uint8 (unsigned
8-bit integer). Pixel values are calculated by summing the following bits:

  * •

00000010 = 2: dark above upper threshold,

  * •

00000100 = 4: dark above upper percentile,

  * •

00001000 = 8: dark below lower threshold,

  * •

00010000 = 16: dark below lower percentile,

so that the minimum and maximum pixel values are 2 and 30, respectively.

Then the planner for nightly mmask is made as plan(MMASK) file, while the
thumbnails for imask and dmask are saved in thumbs_mmask_plan/ inside the log
directory, respectively. The user can remove certain imask or dmask for making
mmask, by changing REMOVEIT column in plan(MMASK) as before.

Figure 5: An example of saved thumbnail images in __logs/<name>/thumbs_mdark/.
It is the master dark frame with 8 s exposure in H-band, as shown in the
title. The statistics are as in Fig. 4.

####  3.4.4 Cell 3: Master Mask

The important default arguments are:

    npr.comb_mmask(combine_kw=dict(combine="sum", reject=None))

It is a simple step: Just summing the original imask and dmask to make mmask.
The user should check if mmask is reasonable at this step. If not, it is
likely that a few dark frames had an abnormality. Any pixels with non-zero
value in mmask will be regarded as bad-pixels. Final mmask, their thumbnails,
and summary are saved in cal-mmask/, thumbs_mmask/, and summ(MMASK) inside the
log directory, respectively. It is recommended not to change any argument
here.

####  3.4.5 Cell 4: Make lv1 and lv2 Planner

The important default arguments are:

    npr.proc_lv1(
        method="median",
        sigclip_kw=dict(sigma=2, maxiters=5),
        fitting_sections=None,
    )
    npr.plan_lv2(maxnsat_oe=5, maxnsat_qu=10)

Running proc_lv1 in this cell will process all frames (except for dark) to
remove the vertical pattern (arguments same as the master dark planner: the
pattern is estimated in the regions fitting_sections, and using the statistic
method, obtained by sigma-clipping with sigclip_kw), and save them to the
level 1 directory with the filename described in Sect. 3.2. Then a summary
file summ(lv1) will be saved in the log directory. Since NIC loses linearity
at around 8000 ADU (Ishiguro et al., 2011), it flags if a frame has many
saturated pixels when making summ(lv1). For each level 1 FITS file, the number
of saturated pixels in the o-/e-ray lit parts are counted. These numbers are
written in NSATPIXO and NSATPIXE columns, respectively, in both the summ(lv1)
and plan(lv2) files.

When making plan(lv2), the REMOVEIT code value is calculated by the following
logic:

  * 0

: Saturated pixels have minor effects.

  * 1

: Number of saturated pixels is larger than maxnsat_oe (default 5) at one or
more of o-ray or e-ray lit parts in this FITS file.

  * 2

: The total number of saturated pixels at the four lit parts in $\texttt{POL-
AGL1}=\theta\,\mathrm{and}\,\theta+45^{\circ}$ frames is larger than
maxnsat_qu (default 10). Both FITS files will have REMOVEIT added by 2.

Those with code 1 will result in an unreliable aperture sum value of at least
one of o-/e-ray. Those with 2 will result in unreliable Stokes’ parameters.
The final code is the sum of these codes. Note that any frame with REMOVEIT!=0
will be ignored from the next (level 2) procedure. That is, NICpolpy
automatically ignores frames with many saturated pixels.

After running it, thumbnail files will also be saved in thumbs/ inside the
level 1 directory. An example is shown in Fig. 6. Using the thumbnail, the
user can modify the auto-generated REMOVEIT column in plan(lv2) to select
which frames to discard or use.

Figure 6: An example of saved thumbnail images in _lv1/<name>/thumbs/. The J,
H, and K band o-/e-ray images are stacked into one image. In the title, the
first line is “<COUNTER> (<OBJECT> <QU><SETID> <POL-AGL1>∘)” (see Sect. 3.2
for description). The last line shows the number of saturated pixels, except
for those that have already masked in mmask. The red crosses in the image
indicate the saturated location. Many pixels in H-band o-ray from the target
are saturated, and the resulting Stokes’ parameter will be unreliable in the
H-band.

####  3.4.6 Cell 5: Make lv2 and lv3 Planner

The important default arguments are:

    npr.proc_lv2(
        cut_wavelength=100,
        bpm_kw=dict(
            size=5,
            med_sub_clip=[-5, 5],   # C_1 in Eq (1)
            med_rat_clip=[0.5, 2],  # C_2 in Eq (1)
            std_rat_clip=[-5, 5],   # C_3 in Eq (1)
            logical="and",
            std_model="std",
            std_section="[50:450, 100:900]",
            sigclip_kw=dict(sigma=3., maxiters=50, std_ddof=1)
        ),
    )
    npr.plan_lv3(use_lv1=False)

Using level 1 data, level 2 data is processed. The user has large degrees of
freedom, but the most important argument is the minimum wavelength of the
Fourier pattern (cut_wavelength) in pixel unit. By default, it is set as 100
pixels. The MBPM algorithm parameters for finding bad pixels can be tuned,
too. However, parameters for MBPM do not result in a large difference unless
they are tuned such that hot pixels or cosmic-rays are missed. proc_lv2 will
generate summ(lv2) in the log directory.

The planner for level 3 will be made after plan_lv3. Based on the exposure
time and filter information in the header, corresponding master dark and flat
frames will be found from the log directory, and will be added as DARKFRM and
FLATFRM columns in plan(lv3) file. The thumbnails of level 2 frames will be
saved in thumbs/ in the level 2 directory (Fig. 7).

Figure 7: An example of saved thumbnail images in _lv2/<name>/thumbs/,
identical to Fig. 6. H-band images are removed (non-zero REMOVEIT in
plan(lv2)). Note that the Fourier pattern is reduced compared to Fig. 6.

The generation of level 2 can take quite a long time. The bad pixel finding
and FFT are rather quick, while most of the time is spent for interpolating
the masked pixels by mmask. It is possible that this Fourier pattern removal
does not work for some cases we could not expect, or the cause for the Fourier
pattern can be removed in the future. Also, because the final Stokes’
parameter values are not changed much from this step777 It is qualitatively
understandable because NIC has a large FWHM in pixel unit. Because of that,
the aperture size should be set large ($\sim 40\,\mathrm{pix}$), and the
positive and negative Fourier pattern offsets within the aperture will roughly
cancel out. , the user may completely ignore level 2 (remove proc_lv2 line)
and just proceed making level 3 from level 1 data, by use_lv1=True.

####  3.4.7 Cell 6: Make lv3

The important default arguments are:

    npr.proc_lv3(flat_mask=0.3, flat_fill=0.3)

This is a standard process of dark subtraction and flat-fielding. In NIC, the
flat image near the edge of each lit part has a very small flat value. Thus,
those edge region will have very high pixel values in the processed image.
Because of this, by default, it fixes flat frame values $<0.3$ (flat_mask) to
a constant value $0.3$ (flat_fill) before flat-fielding. From level 3, each
FITS frame is split into two FITS files, for o-ray and e-ray, specified by a
suffix in the file name (Sect. 3.2). This greatly reduces the size of each
image (from 4.2 MB to two files of 0.28 MB). Finally, using mmask for the
corresponding filter (detector), bad pixel interpolation is conducted. It
prioritizes the vertical direction over the horizontal direction when the
masked region has an identical length to both directions (so that the result
is closer to IRAF.PROTO.FIXPIX). After processing, the summary file summ(lv3)
is made in the log directory.

####  3.4.8 Cell 7: lv4 Planner for Fringe Frames

The important default arguments are:

    npr.plan_ifrin(add_crrej=True)

The fringe frames are selected by finding any frame that has header OBJECT
value ends with _sky. For fringe frames, the thumbnails are saved into
thumbs_mfrin_plan/ in the log directory (Fig. 8). Also, the the planner
plan(IFRIN) is made in the log directory. In the planner, there are columns
for parameters of astroscrappy (cosmic-ray rejection) for each image.
Currently, Gaussian or Moffat convolution options for fine structure
calculation are not implemented. They are not important, as the image is
intended to be a blank sky frame without any source. Among the parameters,
sepmed=False makes astroscrappy runs similar to the original L.A.Cosmic, and
we found this is essential for properly removing cosmic-rays on NIC. Then
sigclip=4.5, sigfrac=5, objlim=1, satlevel=30000, cleantype="medmask",
fs="median" are used by default, with proper gain and readout noise values for
each detector.

Figure 8: An example of saved thumbnail images in
__logs/<name>/thumbs_mfrin_plan/. From the title, this is an H-band e-ray
blank sky frame, taken nearby the object Elias19, at POL-AGL1 = 0.0. The
statistics on the right part are identical to Fig. 4.

####  3.4.9 Cell 8: Make lv4 Fringe & Master Fringe Planner

There are no important arguments the user can tune. It makes level 4 image of
fringe frames based on plan(IFRIN), saves them into frin/ sub-directory inside
the level 4 directory, and makes the next planner file plan(MFRIN). In
plan(MFRIN), the most important column is FRINCID, which stands for fringe
combination identifier. Any frames with the same FRINCID will be combined into
one single master fringe frame. By default, it is set
<FILTER>_<OBJECT>_<id>_<OERAY>. Except for <id>, they are as in Sect. 3.2.
<id> is an integer that separates blank sky exposures that are separated by
another exposure. For example, blank sky exposures before and after the
exposures of the main target will be separated by <id> of 001 and 002
automatically. By default, POL-AGL1 The user can put non-zero value in
REMOVEIT to ignore it for master fringe combination.

####  3.4.10 Cell 9: Make Master Fringe

The important default arguments are:

    npr.comb_mfrin(
        combine_kw=dict(
          combine="med",
          reject="mm",
          n_minmax=(0, 1),
          scale="avg_sc",
          scale_section="[20:100, 20:400]",
          scale_to_0th=False,
        ),
    )

Based on the planner, master fringe is made by combining the same-FRINCID
frames. By default, each fringe frame is divided by sigma-clipped average
value (scale) within scale_section, so that the sky level variation is
removed. By default, in IRAF.IMMATCH.IMCOMBINE, “scale factors are normalized
so that the first input image has no scaling”, which is reproducible by
scale_to_0th=True888 Note that the results may be different from the IRAF
result because IRAF calculates statistics by selecting only a few pixels from
the region to reduce computation time and memory usage, and this behavior is
hard to modify. . This is useful when Poisson noise calculation is important.
By default, NICpolpy scales each image before combination. reject="mm",
n_minmax=(0,1) means it rejects 0 smallest and 1 largest pixel values, before
median (combine="med") combination.

The combined frames are saved in the log directory (cal-mfrin/) with
thumbnails (thumbs_mfrin/) and summary file summ(MFRIN). An example of a
combined master fringe thumbnail is shown in Fig. 9.

Figure 9: An example of saved thumbnail images in __logs/<name>/thumbs_mfrin/.
This is a result of combining four fringe frames, including that in Fig. 8.
The statistics on the right part are identical to Fig. 4.

####  3.4.11 Cell 10: lv4 Planner

The important default arguments are:

    npr.plan_lv4(add_crrej=True, add_mfrin=True)

The final level 4 processing is cosmic-ray removal and fringe subtraction.
Those can be turned on or off by add_crrej and add_mfrin. Here, we used
add_mfrin=False. The planner plan(lv4) is made, and it contains multiple
columns related to the cosmic-ray rejection algorithm as in plan(IFRIN) (see
Sect. 3.4.8). If add_mfrin=True, a column named FRINFRM is added to indicate
the corresponding fringe frame.

####  3.4.12 Cell 11: Make lv4

The important default arguments are:

    npr.proc_lv4(
        frin_bezels=20,
        frin_sep_kw=dict(minarea=np.pi*5**2),
        frin_scale_kw=dict(sigma=2.5),
    )

The important arguments are the parameters related to fringe subtraction.
First, NICpolpy loads the corresponding fringe frame based on FRINFRM in
Plan(lv4). Then it crops the 20-pixel edges (frin_bezels). Within the central
region, it finds any extended bright object using SExtractor algorithm (Bertin
et al., 1996). The variance map is calculated by proper gain and readout noise
parameters, and any object above the 5-sigma level and segment area of
equivalent circular radius of 5 pixels (frin_sep_kw) is found. The fringe
frame is then scaled based on these pixels (sigma-clipped median is used;
frin_scale_kw) and subtracted from the science exposure image. Finally,
cosmic-rays are removed based on the parameters specified in plan(lv4).

####  3.4.13 FITS Header

Finally, each FITS file contains much information on how it is processed
throughout NICpolpy. A few keywords are inserted:

  * •

[MIN/MAX]V<iii>[O/E] (e.g., MINV001E): <iii>-th minimum or maximum pixel
values in level 1 image in o-/e-ray region.

  * •

NSATPIX[O/E]: The number of saturated pixels in o-/e-ray region, based on
level 1 image (Sect. 3.4.5).

  * •

FFTCUTWL: The minimum wavelength of Fourier pattern removal used for level 2
(Sect. 3.4.6)

  * •

LViFRM: The path to the level i frame related to this file.

  * •

[DARK/FLAT/FRIN]FRM: The master dark, flat, or fringe frame used for this
file.

  * •

MASKFILE: The path to mmask file.

  * •

MASKNPIX: The number of masked (fixed) pixels in level 3 (Sect. 3.4.7).

  * •

CRNFIX: The number of pixels fixed by cosmic-ray rejection algorithm in level
4 (Sect. 3.4.12).

  * •

others: YYYYMMDD, COUNTER, SETID, OERAY, PROCESS as in Sect. 3.2.

Also, as COMMENT and HISTORY, most of the minor logs for this specific file
are saved, with options used, timestamp, and the time taken for the step.
Below is the last part of the header of a level 4 image:

    HISTORY ================== Level 1 (vertical pattern removal) ==================
    HISTORY Extrema pixel values found N(smallest, largest) = (5, 5) excluding mask
    HISTORY (__logs/SP_20190417/cal-mmask/mmask_H.fits) and bezel: ((20, 20), (20, 2
    HISTORY 0)) in xyz order. See MINViii[OE] and MAXViii[OE].
    HISTORY Saturated pixels found based on satlevel = 8000, excluding mask (__logs/
    HISTORY SP_20190417/cal-mmask/mmask_H.fits) and bezel: ((20, 20), (20, 20)) in x
    HISTORY yz order. See NSATPIX and SATLEVEL.
    HISTORY ..................................(dt = 0.009 s) 2022-09-25T09:32:30.900
    HISTORY Vertical pattern subtracted using [’[:, 100:250]’, ’[:, 850:974]’] by ta
    HISTORY king median with sigma-clipping in astropy (v 5.1), given {’sigma’: 2, ’
    HISTORY maxiters’: 5}. (using pixel mask) maskpath=’__logs/SP_20190417/cal-mmask
    HISTORY /mmask_H.fits’
    HISTORY ..................................(dt = 0.019 s) 2022-09-25T09:32:30.920
    COMMENT [yfu.update_process] Standard items for PROCESS includes B=bias, D=dark,
    COMMENT  F=flat, T=trim, W=WCS, O=Overscan, I=Illumination, C=CRrej, R=fringe, P
    COMMENT =fixpix, X=crosstalk.
    COMMENT User added items to PROCESS: v=vertical pattern.
    HISTORY ------------------------------------------------------------------------
    HISTORY ================== Level 2 (Fourier pattern removal) ===================
    HISTORY FIXPIX on the left part of the image
    HISTORY ..................................(dt = 0.512 s) 2022-09-25T09:34:56.037
    HISTORY Median filtered (convolved) frame calculated with {’size’: 5, ’footprint
    HISTORY ’: None, ’mode’: ’reflect’, ’cval’: 0.0, ’origin’: 0}
    HISTORY ..................................(dt = 0.168 s) 2022-09-25T09:34:56.207
    HISTORY Sky standard deviation (MB_SSKY) calculated by sigma clipping at MB_SSEC
    HISTORY T with {’sigma’: 3.0, ’maxiters’: 50, ’std_ddof’: 1}; used for std_ratio
    HISTORY  map calculation.
    HISTORY ..................................(dt = 0.020 s) 2022-09-25T09:34:56.228
    HISTORY [medfilt_bpm] Median-filter based Bad-Pixel Masking (MBPM) applied.
    HISTORY ..................................(dt = 0.004 s) 2022-09-25T09:34:56.232
    HISTORY FFT(left half) to get pattern map (see FFTCUTWL for the cut wavelength);
    HISTORY  subtracted from both left/right.
    HISTORY ..................................(dt = 0.002 s) 2022-09-25T09:34:56.241
    COMMENT User added items to PROCESS: f=fourier pattern.
    HISTORY Vertical pattern subtracted using [’[:, 100:250]’, ’[:, 850:974]’] by ta
    HISTORY king median with sigma-clipping in astropy (v 5.1), given {’sigma’: 2, ’
    HISTORY maxiters’: 5}. (using pixel mask) maskpath=’__logs/SP_20190417/cal-mmask
    HISTORY /mmask_H.fits’
    HISTORY ..................................(dt = 0.019 s) 2022-09-25T09:34:56.263
    HISTORY ------------------------------------------------------------------------
    HISTORY [yfu.darkcor] Dark subtracted (DARKFRM = __logs/SP_20190417/cal-mdark/md
    HISTORY ark_H_8.0s.fits)
    HISTORY ..................................(dt = 0.001 s) 2022-09-25T09:38:34.343
    HISTORY [yfu.flatcor] Flat pixels with ‘value < flat_mask = 0.3‘ are replaced by
    HISTORY  ‘flat_fill = 0.3‘.
    HISTORY .................................................2022-09-25T09:38:34.345
    HISTORY [yfu.flatcor] Flat corrected (FLATFRM = __logs/SP_20190417/cal-mflat/mfl
    HISTORY at_H_20180507-lv1.fits)
    HISTORY ..................................(dt = 0.002 s) 2022-09-25T09:38:34.346
    HISTORY [yfu.imslice] Sliced using ‘trimsec = ’[741:890, 331:760]’‘: converted t
    HISTORY o (slice(330, 760, None), slice(740, 890, None)).
    HISTORY ..................................(dt = 0.001 s) 2022-09-25T09:38:34.352
    HISTORY [fixpix] Pixel values interpolated.
    HISTORY ..................................(dt = 0.191 s) 2022-09-25T09:38:34.736
    HISTORY Cosmic-Ray rejection (CRNFIX=0 pixels fixed) by astroscrappy (v 1.1.1.de
    HISTORY v8+g783f217). Parameters: {’gain’: 9.8, ’readnoise’: 36, ’sigclip’: 4.5,
    HISTORY  ’sigfrac’: 5.0, ’objlim’: 1, ’satlevel’: 30000, ’niter’: 4, ’sepmed’: F
    HISTORY alse, ’cleantype’: ’medmask’, ’fsmode’: ’median’}
    HISTORY ..................................(dt = 0.041 s) 2022-09-25T09:42:40.204

##  4 Polarimetric Dome Flat Frames

A large number of polarimetric dome flat frames were obtained on UT 2020-06-03
(total 640 frames per filter; 8 INSROT and 4 POL-AGL1 combinations and 20 FITS
frames per each combination). These frames are useful to check if there is any
region in the RoI the observer should avoid putting the scientific targets.
Also, these flat frames were used to re-calculate the gain factor (Sect. 5).

The frames underwent these processes:

  1. 1.

Reduce all frames to level 1.

  2. 2.

Split o-/e-ray regions by cropping only the lit parts (Fig. 1).

  3. 3.

Scale each frame by “median value/median value of the first frame”.

The master flat frame is obtained by considering all 640 frames per filter.
The flat frames were also grouped by different combinations: 20 frames per
each (INSROT, POL-AGL1) pair, 80 frames per each INSROT, and 160 frames per
each POL-AGL1, which result in 32, 8, and 4 groups per filter, respectively.
The median and sample standard deviation of each pixel after 3-sigma clipping
are saved. The signal-to-noise (SNR) map is generated by dividing the median
map by the sample standard deviation map.

The master flats (normalized) and SNR maps are shown in Fig. 10. Two
interesting features are found: A common doughnut-shaped region at the top-
right corner in SNR maps and a large doughnut-shaped region in the J-band
o-ray region in the mid-right part. These are visible only in the SNR map but
not in the master flat map. These doughnuts could also be found in the flat
frames taken on UT 2018-05-07 using a similar technique.

Figure 10: The master flat frames for UT 2020-06-03 data. The top row (a) is
the normalized flat, and the bottom row (b) is the signal-to-noise ratio
distribution.

An interesting feature appeared for flats grouped based on INSROT. Define the
ratio map as

$r_{1}(\texttt{INSROT},\texttt{POL-
AGL1})=\frac{\mathrm{comb}(\texttt{INSROT},\texttt{POL-
AGL1})}{\mathrm{comb}(\mathrm{any},\texttt{POL-AGL1})}$ (3)

for the given (INSROT, POL-AGL1) pair, i.e., the combined image for the given
(INSROT, POL-AGL1) pair divided by the image combining all frames for the same
POL-AGL1, to see the effect of INSROT. An example for J-band o-ray region for
POL-AGL1 = 0 are shown in Fig. 11. The relative location of the brighter and
darker part of the common doughnut rotates with increasing INSROT. Hence, it
is likely that dust particles or defects in the common optical path are
responsible for the common doughnuts (top-right corner). The master flat may
be different from the true flat at the observation by up to $\sim
2\mathrm{-}3\%$ at such regions, depending on the INSROT at the exposure. For
safety, any observer is encouraged to put their scientific object such that
they are not affected by these doughnuts. Also, the vertical cut profile has a
gradient of $\sim\pm 1\%$, changing as a function of INSROT. The change in the
lit part is seen, too.

Figure 11: The $r_{1}$ ratio map for J-band with HWP rotator angle $0^{\circ}$
and o-ray region only. The numbers on top are the INSROT values of the group.
The vertical and horizontal graphs show the cut profiles along red lines in
the flat images. The red lines in the graphs are the cut profile of $\pm
10\,\mathrm{pixel}$ around the red lines in the image, and black dotted lines
are the minimum and maximum pixel values within the $\pm 10\,\mathrm{pixel}$
region. The ticks for the images indicate 100 pixels, while the ticks for the
cut prociles are $0.98$, $1.00$, and $1.02$.

A similar ratio map was made to investigate the effect of POL-AGL1:

$r_{2}(\texttt{INSROT},\texttt{POL-
AGL1})=\frac{\mathrm{comb}(\texttt{INSROT},\texttt{POL-
AGL1})}{\mathrm{comb}(\texttt{INSROT},\mathrm{any})}~{}.$ (4)

In $r_{2}$ maps, the dust doughnuts disappear. This is likely because the
dark-bright features of doughnuts cancel out when combined over the INSROT.
The combined flats show some large-scale pattern in the final flat, with
fluctuation $\lesssim\pm 0.5\%$ level for all three bands, except for the edge
regions.

To summarize, the users are recommended to avoid the top-right and middle-
right regions of the RoI to best use the three-band data. Also, considering
the bad pixel mask based on both the dark and flat (Fig. 2), it is better to
avoid the top-left part due to the over density of bad pixels in H-band
(e-ray) and K-band (e-ray). These regions may change over time as new defects
appear.

The results for different cases (figures and videos for both $r_{1}$ and
$r_{2}$) are available in the supporting material in Appendix B.

##  5 Gain and Readout Noise: Revisited

The gain (the conversion factor) and readout noise are important parameters in
the error analysis and data reduction (e.g., used in fringe frame scaling,
cosmic-ray rejection, and object detection algorithms). They were first
reported in Ishiguro et al. (2011). According to them, the expected noise in
the bias frame corresponds to the readout noise divided by the gain factor,
i.e., 5.4, 7.7, and 8.8 ADUs for J, H, and K-band, respectively. However, the
sample standard deviations in the dark area in level 1 images are smaller
(roughly 4 ADU for all detectors). This was the motivation for re-calculating
the two detector parameters. We also briefly tested if there is any hint that
these parameters are different per pixel.

If the two flat images have pixel values of $F_{1}$ and $F_{2}$ and two biases
have $B_{1}$ and $B_{2}$ (all in ADU), the gain and readout noise are

$g=\frac{(\bar{F}_{1}+\bar{F}_{2})-(\bar{B}_{1}+\bar{B}_{2})}{\sigma^{2}_{F_{1}-F_{2}}-\sigma^{2}_{B_{1}-B_{2}}}~{}\mathrm{[e/ADU]}$
(5)

and

$R=g\frac{\sigma_{B_{1}-B_{2}}}{\sqrt{2}}~{}\mathrm{[e]}~{}.$ (6)

Here $\bar{X}$ means the average of all the pixels in the frame $X$, and
$\sigma_{X}$ is the true standard deviation of the frame $X$, estimated from
the sample standard deviation
$\sigma_{X}\approx[(\sum_{i}(X_{i}-\bar{X})^{2})/(N-1)]^{1/2}$. This method is
called Janesick’s method. The proof is given in Appendix A.

In the case of NIC, the bias is effectively 0, and the standard deviation of
the bias-difference map is obtained by the difference map of short-exposure
darks to minimize the Poisson noise effect. We selected 20 level 1 dark frames
with 2 s exposure time taken on UT 2019-10-22 for each of the three filters
(detectors) to calculate $R/g$. Next, the flat frames taken on UT 2020-06-04
are used for $g$ calculation so that both $g$ and $R$ can be obtained. From
each dark and flat frame, the bad pixels (imask, Sect. 2.1) are masked.

###  5.1 $R/g$

Here, we discuss how we tested the pixel-to-pixel and quadrant-to-quadrant
variabilities of the $R/g$ value. Then the absolute values were obtained in
ADU for each detector.

####  5.1.1 Pixel-to-Pixel Variability

Are the $R/g$ values dependent on each pixel, i.e., each pixel has different
readout noise per gain? Say the sample standard deviation of the pixel at
(i,j) among the 20 dark frames is $S^{ij}$ and the true standard deviation is
$\sigma^{ij}$ (both in ADU). Since the dark frames have a very short exposure
time and most hot pixels are masked, the only dominant source of the scatter
is readout noise or $R/g$. If the hypotheses are set as

$H_{0}:\sigma^{ij}=\sigma^{i^{\prime}j^{\prime}}\quad;\quad
H_{1}:\sigma^{ij}\neq\sigma^{i^{\prime}j^{\prime}}~{},$ (7)

we can test the null hypothesis ($H_{0}$) using the F-statistic (e.g., Walpole
et al. (2011), Sect. 4.5). Defining

$F:=\frac{S_{1}^{2}}{S_{2}^{2}}\times\frac{\sigma_{2}^{2}}{\sigma_{1}^{2}}\sim
F_{n_{1}-1,n_{2}-1}~{},$ (8)

where $\sim$ means the statistic in the left-hand side follows the
distribution in the right-hand side, and $F_{n_{1}-1,n_{2}-1}$ means the
F-distribution of degrees of freedom $(n_{1}-1,n_{2}-1)$. $n_{1}$ and $n_{2}$
are the number of sample used to obtain $S_{1}$ and $S_{2}$, respectively.
Under $H_{0}$, the value of the statistic $f=S_{1}^{2}/S_{2}^{2}$. Fig. 12
visually depicts how the rejection process works: If our calculated $f$ value
is inside one of those hatched regions, our null hypothesis is rejected with
significance level $\alpha=0.10$, i.e., $\sigma_{1}\neq\sigma_{2}$. If our
statistic is not in the hatched regions, we can make no conclusion under the
significance level $\alpha=0.10$. A number that is more frequently used is the
P-value: It measures the degree of rarity of the observation ($f$ value) under
$H_{0}$. Consider two cases where one obtains $f=f_{0.20;n_{1}-1,n_{2}-1}$ and
$f=f_{0.02;n_{1}-1,n_{2}-1}$. the P-value will be 0.4, and 0.04, respectively,
which means the latter case is rarer than the former one. Thus, the null
hypothesis is rejected in the latter case ($H_{0}$ is unlikely to be true
because such a rare observation was made), while it is not in the former case.

Figure 12: An example of F-distribution probability distribution function
(pdf) with degrees of freedom $(n_{1}-1,n_{2}-1)=(19,19)$. The rejection
regions of significance level $\alpha=0.10$ for a two-tailed test are shown by
two hatched regions.

Because we have $\sim 10^{6}$ pixels at each filter, we have $\sim 10^{12}$
possible F-statistic values (pairs). To save time, we randomly selected $\sim
10^{5}$ pairs, calculated the P-value for each pair, and saw the distribution
of the P-values. Fig. 13 shows two trials of such test on K-band detector dark
frames. We found no significant population of the histogram at small P-values,
even at multiple trials and different detectors. This means we cannot reject
the null hypothesis in the vast majority of cases. Therefore, we have no clear
evidence to conclude that $R/g$ is pixel-dependent.

Figure 13: Two examples of P-value calculations and their distributions for
$124,750$ random pixel pairs. The settings for obtaining these two figures are
identical, except these two simulations samples different random pixel pairs.
The three vertical lines from left to right mean significance s of 0.01, 0.05,
and 0.10, respectively.

####  5.1.2 $R/g$ Value

As a quick test, we selected two random dark frames in K-band, cropped only
the o-ray lit part, and pixels that have values outside the 5-sigma for each
frame are masked. These two dark frames are subtracted, as in Fig. 14. The
histogram shows the distribution is very close to a normal distribution as one
would expect, and the standard deviation divided by $\sqrt{2}$, i.e.,
$R/g=3.67\,\mathrm{ADU}$ for this specific case.

Figure 14: The difference between two random K-band dark frames (left; blue
and red are $-10$ and $+10$ ADU, respectively), and the pixel value
distribution histogram (right). The fitted Gaussian function is shown in the
red dashed line with its standard deviation information.

Since there are total 20 dark frames, one can make $\frac{n(n-1)}{2}=190$
pairs of dark frames and see how the $R/g$ is distributed. The result is shown
in Fig. 15. Regardless of the regions we chosen, the final $R/g$ is well
constrained to $\sim 3.7\,\mathrm{ADU}$ for all three detectors, with scatter
much less than $0.1\,\mathrm{ADU}$.

Figure 15: The distribution of measured $R/g$ for three detectors. The
measurement used all pixels within the 50-pixel edge (a) and only the o-/e-ray
lit parts (b). The box-whisker plots with median are shown, and the violin
plot is overplotted. The boxes indicate the first ($Q_{1}$) and third
($Q_{3}$) quantiles, and the black whiskers extend until
$\min\\{\min\\{y\\},Q_{1}-1.5(Q_{3}-Q_{1})\\}$ to the bottom and
$\max\\{\max\\{y\\},Q_{3}+1.5(Q_{3}-Q_{1})\\}$ to the top, where $\min\\{y\\}$
and $\max\\{y\\}$ are the minimum and maximum of the data, respectively. Black
circles show the outliers (fliers) outside the whisker.

####  5.1.3 Quadrant-to-Quadrant Variability

Finally, we tested if the $R/g$ depends on the quadrants. The technique, the
same as Fig. 15, is applied to each quadrant and checked if any clear
difference could be seen (Fig. 16). Although the distribution is different for
quadrants in K-band, the difference is only $\sim 0.1\,\mathrm{ADU}$.
Therefore, we conclude there is no quadrant dependency in $R/g$.

Figure 16: The distribution of measured $R/g$ for three detectors at each
quadrant. The plots in the upper left, upper right, lower left, and lower
right show the $R/g$ measurement distribution of 190 pairs for the upper left,
upper right, lower left, and lower right quadrants, respectively.

###  5.2 Gain (The Conversion Factor) & Readout Noise

Hinted by the $R/g$ calculation, we now assume both the gain and readout noise
are constant values for all pixels in each detector. As before, NIC has no
bias, so $\bar{B}=0$. Since $R/g$ is a constant, $\sigma^{2}_{B_{1}-B_{2}}$ in
Eq. (5) can be replaced by $2(R/g)^{2}$, and

$F^{ij}\sim\mathcal{N}\left(f^{ij},\frac{f^{ij}}{g}+2\left(\frac{R}{g}\right)^{2}\right)~{},$
(9)

for the flat frame pixel value variable $F^{ij}$ and its measured value
$f^{ij}$ at the pixel (i,j). Here, $\mathcal{N}(\mu,\sigma^{2})$ is the normal
distribution with mean $\mu$ and variance $\sigma^{2}$. The two terms in the
variance are the Poisson and readout noise terms. Say the true pixel value
(i.e., flux) of $F^{ij}$ is constant. Then one can calculate the average of
the location $\bar{F}^{ij}$ and sample standard deviation $S^{ij}$ to obtain

$g=\frac{\bar{F}^{ij}}{S^{ij}-(R/g)^{2}}~{}.$ (10)

Here, $R/g=3.7\,\mathrm{ADU}$ is already obtained previously and found to be
constant.

For the assumption ($F^{ij}$ is constant) to hold, one has to select flat
frames taken at a similar time, the same HWP angles (POL-AGL1), and the same
INSROT because we only used the polarimetric flat for this calculation. The
flat flux changes during the exposures (Fig. 17). The large changes after
every 80 frames are due to the change in INSROT, and the scatter within the
same INSROT is due to the change in POL-AGL1 (semi-periodic for every 4
frames) and/or the true change in the dome flat.

Figure 17: The variability of flat counts in J-band (a, b), H-band (c, d), and
K-band (e, f), for o-ray (a, c, e) and e-ray (b, d, f), respectively. The
black and blue lines show the sigma-clipped average and median, respectively
(left ordinate) and the lines show the sigma-clipped standard deviation (right
ordinate), for the RoIs (green regions in Fig. 1). The INSROT on top shows the
instrument rotator angle at the exposures.

There are 20 consecutive dome flat exposures for the same POL-AGL1 and same
INSROT for each filter. For a , a few flats with the same (POL-AGL1, INSROT)
were selected, and the change in the pixel statistics was calculated. Fig. 18
shows two such cases. The pixel value changes by $\lesssim 1\%$ level. Hence,
for those 20 frames, it can be assumed that $F^{ij}$ is constant.

Figure 18: The change in flat frame pixel statistics after sigma-clipping. The
black and blue lines show the average and median, respectively (left ordinate)
and the red lines show the standard deviation (right ordinate) for the o-ray
RoI (the green region in Fig. 1).

The median and standard deviation were obtained for each pixel for 20 frames
having the same (POL-AGL1, INSROT). Then any pixels with too high standard
deviation or too low pixel values are removed. Next, the gain factor is
calculated using Eq. (10) to obtain a 2-D gain map. For all the 32
combinations of POL-AGL1 and INSROT, the 2-D map is generated for three
detectors. These 32 maps are then median combined with 3-sigma clipping. The
final gain map and the gain value distributions are shown in Fig. 19.

Figure 19: The 2-D gain maps for J-band (top), H-band (middle), and K-band
(bottom). The white pixels are the masked locations (NaN value). The right
column shows the corresponding gain factor distribution ([e/ADU]), best-fit
Gaussian curve, and the previous result (Ishiguro et al., 2011). In the legend
are the best-fit mean and standard deviations.

###  5.3 Summary and Comparison with Previous Work

We have shown that $R/g$ is independent of location on the detector and very
similar for all three detectors of NIC. Whether both $R$ and $g$, not only
$R/g$, are indeed constant would be confirmed separately in future work.
However, if we assume the two parameters are indeed constant, we obtain gain
and readout noise values for all detectors. In Table 2, the gain and readout
noise parameters from this work are compared with the previous work (Ishiguro
et al., 2011).

|  | Ishiguro et al. (2011) | This work | Adopted value
---|---|---|---|---
Parameter | Detector | Motor ON | Motor OFF | ($\mathrm{mean}\pm\mathrm{std}$) | in NICpolpy
| J-band | $4.7\pm 0.5$ | $5.4\pm 0.5$ |  |
$R/g$ [ADU] | H-band | $8.8\pm 0.6$ | $7.7\pm 0.4$ | $3.7(\pm<0.1)$ | $(3.7)$
| K-band | N/A | $8.8\pm 0.6$ |  |
| J-band | $9.8\pm 0.2$ | $9.2\pm 0.2$ | $9.9\pm 1.0$ | $9.9$
$g$ [e/ADU] | H-band | $9.5\pm 0.2$ | $9.8\pm 0.2$ | $9.8\pm 0.9$ | $9.8$
| K-band | N/A | $9.4\pm 0.2$ | $9.5\pm 0.9$ | $9.5$
| J-band | $46\pm 5$ | $50\pm 4$ | $37\pm 4$ | $37$
$R$ [e] | H-band | $84\pm 5$ | $75\pm 4$ | $36\pm 3$ | $36$
| K-band | $(>250)$ | $83\pm 5$ | $35\pm 3$ | $35$
Table 2: The parameter values. $R/g$ in Ishiguro et al. (2011) is calculated
from $R$ and $g$ with error-propagation, while the $R$ for this work is
calculated from $R/g$ and $g$. Error-propagation calculations were done
assuming null covariance.

In the table, the last column shows the values adopted in NICpolpy (callable
as nic.GAIN and NIC.RDNOISE). The official values had been the values in the
column “Motor OFF”, which is the measurement with the shutter motor power
turned down (data taken on UT 2011-07-28).

This discrepancy is likely due to the method used for determining $R$. They
assumed half the variance the flat difference image is a linear function of
signal (Fig. 2 of Ishiguro et al. (2011)), the Poisson noise regime. The
readout noise (or the variance without light) is obtained by extrapolating it
to zero signal. This extrapolation, however, works only if the pixels in the
flat frame used to calculate variance have nearly identical pixel sensitivity
(so that they indeed have the same true flat value). Moreover, in reality,
there is always the $R/g$ term in the variance, so extrapolation using the
equation ignoring $R/g$ may have in an erroneous estimation of $R$. However,
we note that the more important value, gain (the conversion factor), was
determined correctly.

##  Appendix A Proof of Janesick’s Method

Consider that flat and the bias images are $F$ and $B$, and any operation
(e.g., addition or subtraction) means a pixel-wise operation. Also, any
artifacts, bad-pixels or cosmic-rays are ignored. If the true bias level is
$b$ in ADU, any bias image will follow a normal distribution with readout
noise:

$B\sim\mathcal{N}\left(b,\left(R/g\right)^{2}\right)~{}\mathrm{[ADU]}~{},$
(A1)

so that

$B_{1}-B_{2}\sim\mathcal{N}\left(0,2\left(R/g\right)^{2}\right)~{}\mathrm{[ADU]}~{}.$
(A2)

Hence,

$\sigma^{2}_{B_{1}-B_{2}}=2(R/g)^{2}~{},$ (A3)

so Eq. (6) is proven. Since the left-hand side here is approximated by the
sample standard deviation, one may write $\sigma^{2}_{B_{1}-B_{2}}\approx$
np.std(B1-B2, ddof=1)**2 in python.

The (raw) flat image consists of photons with dark plus bias level. Therefore,
if $f$ is the true flat level _plus dark_ in ADU, any flat will follow a
Poisson distribution, which is roughly identical to a normal distribution when
pixel values are high:

$F\sim\mathcal{N}\left(f+b,\frac{f}{g}+\left(\frac{R}{g}\right)^{2}\right)~{}\mathrm{[ADU]}$
(A4)

which means

$F_{1}-F_{2}\sim\mathcal{N}\left(0,2\frac{f}{g}+2\left(\frac{R}{g}\right)^{2}\right)~{},$
(A5)

and

$(F_{1}+F_{2})-(B_{1}+B_{2})\sim\mathcal{N}\left(2f,2\frac{f}{g}+4\left(\frac{R}{g}\right)^{2}\right)~{}.$
(A6)

The $f/g$ terms are the Poisson noise term, and the $R/g$ terms are the
readout noise term, respectively. From these,

$\sigma^{2}_{F_{1}-F_{2}}-\sigma^{2}_{B_{1}-B_{2}}\approx 2\frac{f}{g}~{},$
(A7)

and

$(\bar{F}_{1}+\bar{F}_{2})-(\bar{B}_{1}+\bar{B}_{2})\approx 2f~{}.$ (A8)

This proves Eq (5). Since the $\sigma$ values are approximated by the sample
standard deviation, $\sigma^{2}_{F_{1}-F_{2}}-\sigma^{2}_{B_{1}-B_{2}}\approx$
np.std(F1-F2, ddof=1)**2 - np.std(B1-B2, ddof=1)**2 and
$(\bar{F}_{1}+\bar{F}_{2})-(\bar{B}_{1}+\bar{B}_{2})=$
np.mean((F1+F2)-(B1+B2)) in python.

##  Appendix B Supporting Materials

The basic calibration files (master flats for UT 2018-05-07 and 2020-06-04 and
imask), flat analyses results (Sect. 4), and the data reduction example with
sample data (Sect. 3) are available via the GitHub
service999https://github.com/ysBach/NICpolpy_sag22sm. The contents are:

  * •

cal-flat_20180507-lv1/: The master flats taken on UT 2018-05-07.

  * •

cal-flat_20200603-lv1/: The master flats taken on UT 2020-06-04.

  * •

masks/: The imask frames (Sect. 2.1).

  * •

flat_analyses: The figures and videos of flat maps and their change over
INSROT or POL-AGL1. sm1.pptx file here describes the details of how the
figures are generated.

  * •

example/: Sample data, reduction code and output log directory contents (Sect.
3).

## References

  * Bradley et al. (2022) Bradley, L., Sipőcz, B., Robitaille, T. et al. 2022, astropy/photutils: 1.5.0, Zenodo, https://doi.org/10.5281/zenodo.6825092
  * McCully et al. (2018) McCully, C., Crawford, S., Kovacs, G. et al. 2018, astropy/astroscrappy: v1.0.5, Zenodo, https://doi.org/10.5281/zenodo.1482019
  * Ishiguro et al. (2011) Ishiguro, M., Takahashi, J., Zenno, T., Tokimasa, N. & Kuroda, T. 2011, Annu. Rep. Nishi-Harima Astron. Obs., 21, 13
  * Takahashi (2018) Takahashi, J., Zenno, T., Saito & Itoh, Y. 2018, Stars And Galaxies, 1, 17
  * Takahashi et al. (2019) Takahashi, J. 2019, Stars And Galaxies, 2, 1
  * van Dokkum (2001) van Dokkum, P. 2001, PASP, 113, 1420
  * Takahashi et al. (2021) Takahashi, J., Itoh, Y., Matsuo, T., Oasa, Y., Bach, Y. & Ishiguro, M. 2021, A&A, 653, A99
  * Bertin et al. (1996) Bertin, E. & Arnouts, S. A&AS, 117, 393
  * Walpole et al. (2011) Walpole, R., Myers, R., Myers, S., & Ye, K. 2011, “Essentials of Probability & Statistics for Engineers & Scientists”