SAG: Stars and Galaxies Vol. 5, id. 4, 2022 December 28
c 2022. Center for Astronomy, University of Hyogo.
⃝

Data Reduction Process and Pipeline for the NIC Polarimetry Mode in
Python, NICpolpy
Yoonsoo P. Bach,1,2,∗ Masateru Ishiguro,1,2,∗∗ Jun Takahashi,3,∗∗∗ and Jooyeon Geem1,21Astronomy Program, Department of Physics and Astronomy, Seoul National University,
Gwanak-ro 1, Gwanak-gu, Seoul 08826, Republic of Korea
2
SNU Astronomy Research Center, Department of Physics and Astronomy, Seoul National University,
Gwanak-ro 1, Gwanak-gu, Seoul 08826, Republic of Korea
3
Nishi-Harima Astronomical Observatory, Center for Astronomy, University of Hyogo,
407–2 Nishigaichi, Sayo-cho, Hyogo 679–5313, Japan∗ysbach93@gmail.com, ∗∗ ishiguro@snu.ac.kr, ∗∗∗ takahashi@nhao.jp(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 asof 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 andin-depth descriptions are provided, as well as the descriptions of algorithms. The dome ﬂat frames (taken
on UT 2020-06-03) were analyzed, and the pixel positions vulnerable to ﬂat error were found. Using thedark and ﬂat frames, the detector parameters, gain factor (the conversion factor), and readout noise arealso updated. We found gain factor and readout noise are likely constants over pixel or “quadrant”.
Key words: methods: data analysis — methods: observational — techniques: image processing —techniques: polarimetric1.

Introduction

In this work, we describe a way to reduce Nishi-Harima
Astronomical Observatory (NHAO) Nishi-harima Infrared
Camera (NIC) polarimetric mode data. The process is implemented to NICpolpy1 , 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 ﬂat corrections, but excluding
photometry/polarimetry and error analysis. We analyzed
the dome ﬂat frames (UT 2020-06-03) and examined if any
pixels have large uncertainty due to imperfect ﬂat correction or hot pixels. Also, the eﬀect 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), wererecalculated using dark frames (UT 2019-11-21) and thedome ﬂat frames. These had already been determined in
Ishiguro et al. (2011) more than a decade ago and should
1

The stable version is registered to the Python Package
Index(PyPI):
https://pypi.org/project/NICpolpy/.
The development version is available via GitHub athttps://github.com/ysBach/NICpolpy.
See section 3 fordetails.

be updated due to a possible secular change. Throughout
this work, the three detectors associated with three ﬁlters,
J, H, and Ks , are referred to as J-, H-, and K-band, respectively.
In section 2, the raw image and reduction processes
implemented to NICpolpy package are described. Thenthe practical usage of the package is shown in section 3,
with additional descriptions of algorithms and implementations. Section 4 describes an analysis of polarimetricdome ﬂat images to ﬁnd any locations vulnerable to ﬂat
ﬁelding uncertainties. Finally, in section 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 ﬁle size
(a ﬁle 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 anarbitrary constant bias (e.g., few thousands) and save as2

Y. P. Bach et al.

Vol. 5,half is shifted by +512 pixel to the right.

part

FITS section

J-band, o-ray[534:683, 306:735]

2.1.

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]

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 theupper-right, upper-left, lower-right, and lower-left quadrants, respectively. For each column within each of theseboxes, a 2-sigma clipped median is calculated. This value
is subtracted from all the pixels in the same column in thesame 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 ﬁrst reduced
to level 1. Then after sigma-clipped median combining
these dark frames for each detector, any pixel with values
larger than max{5, Dmed + 5 σD } ADU or smaller thanmin{−5, Dmed − 5 σD } ADU are ﬂagged as bad-pixel foreach detector. Here, Dmed and σD are the 3-sigma clipped
median and standard deviation of the combined darkframe, respectively. Secondly, the sigma-clipped sample
standard deviation was calculated for each pixel position,
and any pixel with high variance was ﬂagged.
On UT 2020-06-03, 640 dome ﬂat 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 ﬂat countis diﬀerent in the two regions by up to ∼ 10%. Then the
pixels with too low value or too high variance are ﬂagged.
Combining all the ﬂagged pixels from these dark andﬂat frames, initial bad-pixel masks were generated foreach detector (total 3 mask ﬁles). These are called theimask in NICpolpy, and shown in ﬁgure 2. As a result,∼ 500 pixels are masked for each light-green RoI in ﬁgure
1. Since the total size after cropping the 10-pixel bezel is
130 × 410 ∼ 50, 000 pixel, ∼ 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 aremasked.

Table 1. The FITS sections of each lit part. The FITS section is in the format of [xmin : xmax ,ymin : ymax ] in 1-indexing(ﬁrst and last pixels inclusive).

an unsigned 16 bit integer (0−32, 767), because anyway
the maximum signal for NIC is < 2 × 104 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 ﬁrst 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 thispattern is removed, the image is called level 2 or lv2.
Then, usual dark subtraction, ﬂat ﬁelding, and bad pixel
interpolation results in level 3 or lv3. Finally, after an
optional fringe subtraction and cosmic-ray removal, theresult is called level 4 or lv4.
Figure 1 shows a sample FITS frame during the reduction steps. The raw frame consist of four parts of size
512 × 512 pixels, and we call them quadrants. As is visible, the light from the exposure falls onto a small portionin the right half of the full array, indicated as green boxes
in ﬁgure 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 sizeof 150 × 430 pixels. Hence, any region well outside theseboxes, including B1, B1’, B2, B2’, FFTo, and FFTe, aredark overscan areas. After processing until level 4, it isrecommended 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 aﬀectedby the imperfect ﬂat ﬁeld and can have arbitrarily large
or small pixel values, and can aﬀect the sky estimationin photometric/polarimetric measurements. We normally
crop out 10 pixel around the edge (so the ﬁnal analyzed
region is 130 × 410 pixels).
In this section, the detailed process to obtain the reduced image (levels 1−4) is described. B1, B1’, B2, andB2’ regions are used to make level 1 data, while FFToand FFTe regions are used for level 2 data. The FFToand FFTe regions have the same size (150 × 1024 pixels),
and are set such they overlap with the lit part when left

2.2.

Level 1 and Initial Bad-Pixel Mask

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, ﬁgure 28 of Ishiguro
et al. 2011 and ﬁgure 5 of Takahashi et al. 2018). Although
the speciﬁc reason is still unknown, this is likely an additive noise due to an imperfectly shielded external electric
signal. After thorough investigations, we could not ﬁnd a
way to remove the noise by tiling columns or rows of thefull 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 foundthat this pattern is quite similar in the left/right quadrants. That is, a direct subtraction of the left part from

id. 4

NICpolpy: NIC Polarimetry Reducer(a) raw(b) schematic
512

512

B1'

B1

3
(c) level 1

512

512

B1'

B1

50

50
124

150

FFTo

FFTe

FFTe

512

512

FFTo

B2'

B2

B2

150
100

100

(d) level 2

oray
FoV
e-ray
RoI

10 px bezels

B2'

150

oray
FoV
o-ray
RoI

430

512

512

150

124

(e) level 3(f) level 4

Fig. 1. An example exposure (2019-10-22 UT, Vesta). (a) The raw data, (b) schematic diagram describing the full NIC frame inpolarimetry 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 ofthe 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), isused 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 ﬂat ﬁelding, and bad-pixel interpolation, level 3 is obtained (e). Finally, a conventionalcosmic-ray rejection technique is applied to obtain level 4 (f). In this example, fringe removal was turned oﬀ 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.

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 diﬃcult to be applied directly, because(i) it takes a long time and (ii) the ﬁne-tuning of the parameters for dark (left) area was diﬃcult. Therefore, weapplied the following strategy.
As described in subsection 2.1, we already made theinitial pixel mask (imask) from dark frames of UT 201911-21. Whenever dark frames are available, those of thelongest 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 thenight, mmask = imask. This process can remove newly de-veloped bad pixels after UT 2019-11-212 . We use dmask
rather than just subtracting dark, because the absolutevalue of hot pixels often varies from frame to frame.
Afterwards, similar to, but slightly diﬀerent from, theprevious work Takahashi et al. (2021), the left half ofthe 1024 × 1024 pixels frame is used to estimate andremove 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 ﬂagged in
mmask are interpolated. This interpolation is similar
to IRAF.PROTO.FIXPIX, prioritizing the vertical direction
2

It is also possible to use the dark regions (outside the lit parts)of science exposure frames to ﬁnd additional bad pixels. Thismight 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.

4

Y. P. Bach et al.
(a) J-bandout ﬁne-tuning of the parameters 3 .
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 pixel are extracted, so that two pattern maps withsizes 150 × 1024 pixels are made. These pattern maps
are subtracted from the right part, overlapping with oray and e-ray RoI, respectively. This image undergoes anadditional vertical pattern subtraction, as in level 1, to remove any remaining vertical pattern (likely smaller than
0.5 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 (<∼ ± 5 ADU, see, ﬁgure 1c). We conﬁrmed that thisprocess does not change the ﬁnal Stokes’ parameter values
while greatly reducing background ﬂuctuation. These reduction processes would be beneﬁcial for extended source
studies, e.g., the Moon, earth shine (Takahashi et al. 2021)
or comet polarimetry.

(b) J-band FoV

1000
700

800
600

600

400

500

200

400

0

300
(c) H-band(d) H-band FoV

1000
700

800
600

600

400

500

200

400

0

300
(e) K-band(f) K-band FoV

1000

2.3.

Master Flat and Level 3 and 4

700

800
600

600

400

500

200

400

0

Vol. 5,

0

250

500

750

1000

300

We found no clear and severe diﬀerence in the combined
ﬂats for diﬀerent POL-AGL1. We hence did not split themaster ﬂat based on the rotator angles (section 4 andsupporting material in Appendix 2). The master ﬂat was
generated by the following steps:
600

700

800

1. sum combines each polarimetric set (four frames atfour POL-AGL1),
2. normalize each of the sum images based on thesigma-clipped median in the o-ray RoI,
3. median combines these by sigma-clipping.

900

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

over the horizontal one (the priority can be modiﬁed 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 replacedby the median-ﬁltered value. Three parameters are calculated from the image in this algorithm, called the Medianﬁlter Bad-Pixel Masking (MBPM):
C1 = I − M , C2 = I/M , C3 = (I − M)/σ(I) . (1)Here, I is the original image, M is the median ﬁlteredimage, and σ(I) is the sigma-clipped standard deviation
at certain user-tunable region(s) of I. In this work, thepixels with |C1 | > 5 [ADU] or C2 ∈/ [0.5, 2] or |C3 | > 5 are
ﬂagged and replaced with the median-ﬁltered value. The
regions for σ(I) is [50:450, 100:900] by default in the FITS
section format. It is not iterated multiple times. The
median-ﬁlter size was 5 × 5 pixels, and a very small real
number is added to M to avoid zero-division. This MBPM
algorithm is orders of magnitude quicker than traditionalcosmic-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 with-

Because the ﬂat counts are not normalized separately foro-/e-ray regions, the ﬁnal ﬂat has slightly diﬀerent values
for the o-/e-ray regions (∼ 1 and ∼ 1.1, respectively). Thisdiﬀerence is not important because NIC samples the ﬂux
twice per each Stokes’ parameter (e.g., POL-AGL1 of 0◦
and 45◦ for q) to cancel out this oﬀset eﬀect.
After removing the instrumental noises (vertical andFourier patterns), a standard preprocessing, including
dark subtraction and ﬂat ﬁelding, was done. Using mmask,some pixels are interpolated. Then each frame is split into
two frames, o-ray and e-ray (each 150 × 430 pixel). Thisis 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 aredone on each image, saving the ﬁnal result as level 4.
3

MBPM is not intended to detect cosmic rays in the presence ofcelestial objects. Unlike algorithms that use edge detection andconvolutions to ﬁnd cosmic rays, MBPM only uses the deviation from the median ﬁltered image to ﬁnd outliers without any
additional convolution or edge detection. Thus, MBPM may
easily ﬂag stellar objects. MBPM is intended to quickly ﬁnd
bad pixels or cosmic rays in dark frames or sky regions where
no point source is expected. MBPM is separately available as afunction at NICpolpy.ysfitsutilpy4nicpolpy.medfilt bpm.

id. 4

NICpolpy: NIC Polarimetry Reducer

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 justiﬁed later in future
work, the fringe subtraction is implemented but omitted
in real data reduction. Each fringe frame is divided by thesigma-clipped average pixel value to eliminate the overallsky brightness change. After this normalization, frames
are combined to make a master fringe frame (for the ﬁnetuning, see sub-subsection 3.4.10). Since the fringe pattern is well known to be variable over time, NICpolpy
automatically combines fringe based on time stamps, too(sub-subsection 3.4.9). The HWP angle, POL-AGL1, is ignored when combining master fringe because we foundthat separating fringe by HWP angle does not improve theﬁnal results. For the cosmic-ray rejection, astroscrappy(McCully et al. 2018), which implements L.A.Cosmic (vanDokkum 2001), is used (version 1.1.1). After the investigation, the parameters similar to L.A.Cosmic are foundto be able to remove cosmic rays eﬃciently.
The blank sky (fringe) frame can approximately be regarded as a dark frame, especially for short exposures andwhen the sky is dark. This, however, is vulnerable to thechange in sky level, spatial polarization patterns in thesky, and typical pixel noise (readout noise and Poissonnoise) 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 ditheringobservation to construct a fringe map. This technique isalso diﬃcult for NIC because of the small ﬁeld-of-view, especially when sky level and/or fringe pattern vary quickly
relative to the exposure time. Finally, it can be shownthat the fringe pattern has a negligible eﬀect on Stokes’
parameter determination for beam-splitter type aperturepolarimetry when the eﬀect of fringe is a few-percent levelof the aperture sum and the true polarization degree of thetarget is <∼ few % (Bach et al. in prep.). Thus, we ignore
fringe in this work.
3.

as well as intermediate output ﬁles. 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 suﬃce most of the users’ needs.
3.1.

# 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 preﬁx that will be
used for result ﬁles (see subsection 3.2 and 3.3). The argument inputs, mflats and imasks indicates the regex-style
for input FITS ﬁles, master ﬂat ﬁles, and imask ﬁles, 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 secThis section describes the practical usage of the package. Since NICpolpy is registered to the Python Package
Index (PyPI), the installation is invoked once typing(2)

on the terminal4 . The development version is available via
GitHub5 . After the installation, calibration frames imask
and ﬂats are required. They are available in a separate
GitHub repository (Appendix 2).
First, the simplest way to use NICpolpy is brieﬂy discussed. Then the output FITS ﬁles and log ﬁles are explained. Finally, the core algorithmic implementations of
reduction processes (section 2) are described step by step,
4
5

The
terminal
must
have
https://pypi.org/project/pip/
https://github.com/ysBach/NICpolpy

pip

Basic Usage

The most recent development of NICpolpy as of writingis happening on Python 3.10 with macOS 12. However,since each part of the code was written OS-agnostic, itwill work without any problem on any recent major OS(Windows, macOS, Linux) and python versions (3.8 orlater).
First, import packages and initialize the reducer by giving the names of the directories where the original dataand calibration data are stored6 :Code Usage and Implementation$ pip install NICpolpy

5

installed:
6This example is also provided in Appendix 2.

6

Y. P. Bach et al.

• 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 (ﬂat-ﬁelding), T (trimmed), C(cosmic-ray rejection), R (fringe subtraction), P (badpixel interpolation), v (vertical pattern removal),
and f (Fourier pattern removal). The processingorder 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, andthe 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 diﬀerent
POL-AGL1 angles, and the counter increases after thecycle. 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).

_ = 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 diﬀerent Cell numbers arerecommended to run after checking the results from theprevious Cell. Note that, except for the input ﬁle directories, the user does not have to specify any speciﬁcparameters and quickly obtain the ﬁnal level 4 ﬁles. For
speciﬁc arguments for ﬁne-tuning, see subsection 3.4. The
comments with xx min or xx sec are the time spent on
Mac Book Pro with m1Pro chip (8 performance, 2 eﬃciency cores with 16 GB RAM).
From now on, the following simplifying notations areused:log directory = logs/<name> ,level <i> directory = lv<i>/<name> ,summ(X) = <name> summ X.csv ,plan(X) = <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 speciﬁc example means the ﬁle
SP 20190417 plan-MDARK.csv, inside the corresponding
log directory, logs/SP 20190417/. After running thecode cell’s above, multiple outputs are generated. Most
importantly, the FITS ﬁles are saved in the corresponding
level directory, lv<i>/<name>. Also, the log information
and all the calibration ﬁles are saved in the log directory,logs/<name>. The contents are described in the following subsections.
3.2.

Output: FITS Files and Level Directories

The resulting ﬁles will be saved in the corresponding level directories. At each level directory, a subdirectory for thumbnails of each FITS frame called
thumbs/, is generated when making the planner for thelevel <i+1> (see subsection 3.4). There can be other subdirectories, 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, andits tree structure is shown in ﬁgure 3. This is made toarchive virtually all necessary information to guarantee
the reproducibility of the user’s work. The calibration ﬁles(master ﬂat and mask ﬁles) are copied and saved, and all
the thumbnail images related to calibration will be saved
here. Inside here are multiple CSV ﬁles and directories7 :• 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 foreach level or master calibration frames, e.g.,
SP 20190417 summ lv0.csv.
• plan(X): The planner before proceeding to eachlevel or master calibration frame generation, e.g.,
SP 20190417 plan-MDARK.csv.

Any FITS ﬁle 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 ﬁlter (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 theheader.

Vol. 5,

3.4.

Implementation Details

Here, the details of implementation and the outputs
from each Cell number are described. As most of thetools do, NICpolpy also may undergo quick and breaking
changes when necessary. So for the most details, such as
tuning output thumbnail ﬁles, which do not aﬀect the scientiﬁc outcomes, it is recommended to look at the source
7

id. 4

NICpolpy: NIC Polarimetry Reducer__logs/SP_20190417/├── SP_20190417_plan-IFRIN.csv├── SP_20190417_plan-MDARK.csv├── SP_20190417_plan-MFRIN.csv├── SP_20190417_plan-MMASK.csv├── SP_20190417_plan-lv2.csv├── SP_20190417_plan-lv3.csv├── SP_20190417_plan-lv4.csv├── SP_20190417_summ_DMASK.csv├── SP_20190417_summ_IMASK.csv├── SP_20190417_summ_MDARK.csv├── SP_20190417_summ_MFLAT.csv├── SP_20190417_summ_MFRIN.csv├── SP_20190417_summ_MMASK.csv├── SP_20190417_summ_lv0.csv├── SP_20190417_summ_lv1.csv├── SP_20190417_summ_lv2.csv├── SP_20190417_summ_lv3.csv├── SP_20190417_summ_lv4.csv├── cal-dmask # 3.2 MB│├── darkmask_mdark_H_8.0s.fits│├── darkmask_mdark_J_8.0s.fits│└── darkmask_mdark_K_8.0s.fits├── cal-imask # 3.2 MB│├── mask_H.fits│├── mask_J.fits│└── mask_K.fits├── cal-mdark # 12.7 MB│├── mdark_H_8.0s.fits│├── mdark_J_8.0s.fits│└── mdark_K_8.0s.fits├── cal-mflat # 12.7 MB│├── mflat_H_20180507-lv1.fits│├── mflat_J_20180507-lv1.fits│└── mflat_K_20180507-lv1.fits├── cal-mfrin # 10 MB / 36 FITS│├── mfrin_H_Elias19_sky_001_e.fits
│
├── ⋮
│
└── mfrin_K_HDE283809_sky_001_o.fits7├── cal-mmask # 3.2 MB│├── mmask_H.fits│├── mmask_J.fits│└── mmask_K.fits├── thumbs_dmask # 68 kB│├── darkmask_mdark_H_8.0s.pdf│├── darkmask_mdark_J_8.0s.pdf│└── darkmask_mdark_K_8.0s.pdf├── thumbs_mdark # 81 kB│├── mdark_H_8.0s.pdf│├── mdark_J_8.0s.pdf│└── mdark_K_8.0s.pdf├── thumbs_mdark_plan # 2.3 MB / 60 PDF│├── h_20190417_0153_8.0_00.0_DARK-PROC-v_.pdf
│
├── ⋮
│
└── k_20190417_0172_8.0_00.0_DARK-PROC-v_.pdf├── thumbs_mfrin # 1 MB / 36 PDF│├── mfrin_H_Elias19_sky_001_e.pdf
│
├── ⋮
│
└── mfrin_K_HDE283809_sky_001_o.pdf├── thumbs_mfrin_plan # 4MB / 144 PDF│├── h_20190417_0029_8.0_00.0_HDE283809_sky-PROC-vDFTPC_q001_e.pdf
│
├── ⋮
│
└── k_20190417_0152_8.0_67.5_Elias19_sky-PROC-vDFTPC_u004_o.pdf├── thumbs_mmask # 75 kB│├── mmask_H.pdf│├── mmask_J.pdf│└── mmask_K.pdf└── thumbs_mmask_plan # 138 kB├── darkmask_mdark_H_8.0s.pdf├── darkmask_mdark_J_8.0s.pdf├── darkmask_mdark_K_8.0s.pdf├── mask_H.pdf├── mask_J.pdf└── mask_K.pdf
13 directories, 324 files
53.5 MB total

Fig. 3. The ﬁnal tree structure of the log directory, logs/SP 20190417, as an example. Total 13 directories, 324 ﬁles, and 53.5
MB. The size of each sub-directory is indicated as comments.

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 thelog directory, cal-imask/ and cal-mflat/, based on
imasks and mflats arguments, are made. These are notmade if arguments are not given (None). Also, the summary of level 0 ﬁles are made as summ(lv0) inside the log
directory. This contains some selected header informationfrom each FITS ﬁle. If imasks or mflats are given, corresponding summary ﬁles are also made: summ(MFLAT) andsumm(IMASK). These two contains little information, andmade 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 ﬁles, it ﬁnds those with OBJECTS ="DARK" in the header. If there is any nightly darkframes, they undergo vertical pattern removal (patternis estimated in the regions fitting sections, and using the statistic method, obtained by sigma-clipping withsigclip kw), and saved into dark/ inside level 1 directory. The planner plan(MDARK) is made. This plannercontains ﬁle name, timestamp, ﬁlter, and exposure time
information. There is a column named REMOVEIT, and if
the user ﬁlls 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 ﬁgure 4. The user may lookat the thumbnails and pinpoint strange dark frames toﬁll in REMOVEIT. For example, at one night, we could ﬁndmultiple dark frames have abnormally high median valuesof > 5 ADU only at K-band based on these thumbnails.
Such frames can be removed before making the masterdark in the next Cell.
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 combinedwith 3-sigma clipping (tunable by combine kw). The resulting master darks are saved in cal-mdark/, while thethumbnails are saved in thumbs mdark/ in the log directory. An example of thumbnails is shown in ﬁgure 5. The
summary of master darks is also saved as summ(MDARK).
For each master dark frame, any pixel that are outside[−10, 50] ADU range, or [0.01, 99.99] percentile range are
ﬂagged (tunable by dark thresh and dark percentile,
respectively).
The generated intermediate dmask ﬁles and theirthumbnails are saved in cal-dmask/ and thumbs dmask/,
and summary ﬁle 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8

Y. P. Bach et al.

Vol. 5,h_20190417_0156_8.0_00.0_DARK-PROC-v_.fits

1000
800

num
: 968256min: -27.0
max: 7955.0
avg: 0.5795739889144897
med: 0.0000
std: 28.219921112060547
madstd
: 4.4478
percentiles : ['0.0100', '99.9000', '99.9900']pct: ['-15.0000', '39.0000', '1368.0470']zmin: -10.0000
zmax: 10.0000
ext_lo: ['-27.0000', '-25.0000', '-23.0000']
ext_hi: ['5311.0000', '5790.0000', '7955.0000']

600
400
200
0

0

250

500

750

1000

Fig. 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 ofpixels, 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 speciﬁed 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.

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 and30, respectively.
Then the planner for nightly mmask is made asplan(MMASK) ﬁle, while the thumbnails for imask anddmask 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.
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 thata 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 tochange 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 patternis estimated in the regions fitting sections, and using the statistic method, obtained by sigma-clipping withsigclip kw), and save them to the level 1 directory withthe ﬁlename described in subsection 3.2. Then a summary ﬁle summ(lv1) will be saved in the log directory.
Since NIC loses linearity at around 8000 ADU (Ishiguro
et al. 2011), it ﬂags if a frame has many saturated pixels when making summ(lv1). For each level 1 FITS ﬁle,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) ﬁles.
When making plan(lv2), the REMOVEIT code value iscalculated by the following logic:
0 : Saturated pixels have minor eﬀects.
1 : Number of saturated pixels is larger thanmaxnsat oe (default 5) at one or more of o-ray or
e-ray lit parts in this FITS ﬁle.
2 : The total number of saturated pixels at the four lit
parts in POL-AGL1 = θ and θ + 45◦ frames is larger
than maxnsat qu (default 10). Both FITS ﬁles will
have REMOVEIT added by 2.
Those with code 1 will result in an unreliable aperture sumvalue of at least one of o-/e-ray. Those with 2 will result in
unreliable Stokes’ parameters. The ﬁnal code is the sumof these codes. Note that any frame with REMOVEIT!=0will be ignored from the next (level 2) procedure. Thatis, NICpolpy automatically ignores frames with many saturated pixels.
After running it, thumbnail ﬁles will also be saved inthumbs/ inside the level 1 directory. An example is shownin ﬁgure 6. Using the thumbnail, the user can modify theauto-generated REMOVEIT column in plan(lv2) to selectwhich frames to discard or use.
3.4.6. Cell 5: Make lv2 and lv3 Planner
The important default arguments are:npr.proc_lv2(

id. 4

NICpolpy: NIC Polarimetry Reducer

9

mdark_H_8.0s.fits

1000
800

num
: 968256min: -24.0
max: 7981.0
avg: 0.5669587850570679
med: 0.0000
std: 27.752891540527344
madstd
: 1.4826
percentiles : ['0.0100', '99.9000', '99.9900']pct: ['-13.0000', '39.8725', '1349.8625']zmin: -4.0000
zmax: 5.0000
ext_lo: ['-24.0000', '-23.5000', '-20.0000']
ext_hi: ['5295.5000', '5800.5000', '7981.0000']

600
400
200
0

0

250

500

750

1000

Fig. 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 ﬁgure 4.

0131 (Elias19 u019 22.5 )Each img ZSCALE'd. Order: Jo, Je, Ho, He, Ko, Ke# of saturated pix excluding MMASK: [0, 0, 29, 0, 0, 0]

400
300
200
100
0

0

200

400

600

800

Fig. 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 ﬁrst line is “<COUNTER> (<OBJECT> <QU><SETID> <POL-AGL1>◦ )” (see subsection 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.
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 haslarge 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 ﬁnding bad pixels can be tuned, too. However, parametersfor MBPM do not result in a large diﬀerence 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 ﬁlter information in theheader, corresponding master dark and ﬂat frames will befound from the log directory, and will be added as DARKFRM
and FLATFRM columns in plan(lv3) ﬁle. The thumbnailsof level 2 frames will be saved in thumbs/ in the level 2
directory (ﬁgure 7).
The generation of level 2 can take quite a long time. The
bad pixel ﬁnding and FFT are rather quick, while mostof 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 thecause for the Fourier pattern can be removed in the future.
Also, because the ﬁnal Stokes’ parameter values are notchanged much from this step8 , the user may completely
8

It is qualitatively understandable because NIC has a large

10

Y. P. Bach et al.

Vol. 5,0131 (Elias19 u019 22.5 )Each img ZSCALE'd. Order: Jo, Je, Ko, Ke# of saturated pix excluding MMASK: [0, 0, 0, 0]

400
300
200
100
0

0

200

400

600

800

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

ignore level 2 (remove proc lv2 line) and just proceedmaking 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 ﬂatﬁelding. In NIC, the ﬂat image near the edge of each lit
part has a very small ﬂat value. Thus, those edge region
will have very high pixel values in the processed image.
Because of this, by default, it ﬁxes ﬂat frame values < 0.3(flat mask) to a constant value 0.3 (flat fill) before
ﬂat-ﬁelding. From level 3, each FITS frame is split into
two FITS ﬁles, for o-ray and e-ray, speciﬁed by a suﬃx in
the ﬁle name (subsection 3.2). This greatly reduces thesize of each image (from 4.2 MB to two ﬁles of 0.28 MB).
Finally, using mmask for the corresponding ﬁlter (detector), bad pixel interpolation is conducted. It prioritizes
the vertical direction over the horizontal direction whenthe masked region has an identical length to both directions (so that the result is closer to IRAF.PROTO.FIXPIX).
After processing, the summary ﬁle 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 ﬁnding 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 (ﬁgure 8).
Also, the the planner plan(IFRIN) is made in thelog directory. In the planner, there are columns forparameters of astroscrappy (cosmic-ray rejection)
FWHM in pixel unit. Because of that, the aperture size should
be set large (∼ 40 pixel), and the positive and negative Fourier
pattern oﬀsets within the aperture will roughly cancel out.

h_20190417_0137_8.0_00.0_Elias19_sky-PROC-vfvDFTPC_q001_e.fits

400
num: 42900min: -131.1795654296875
max: 159.13392639160156
avg: 111.20191192626953
med: 111.2935
std: 10.510380744934082
madstd: 9.6655
percentiles : ['1.0000', '95.0000', '99.0000']pct: ['87.6504', '127.5034', '134.7212']zmin: 54.6311
zmax: 147.2469
ext_lo: ['-131.1796', '-120.4348', '-96.0463']
ext_hi: ['156.4678', '158.5584', '159.1339']

300
200
100
0

0

100

Fig. 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 theobject Elias19, at POL-AGL1 = 0.0. The statistics on theright part are identical to ﬁgure 4.

for each image. Currently, Gaussian or Moﬀat convolution options for ﬁne 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 andreadout noise values for each detector.
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 onplan(IFRIN), saves them into frin/ sub-directory inside the level 4 directory, and makes the next plannerﬁle plan(MFRIN). In plan(MFRIN), the most importantcolumn is FRINCID, which stands for fringe combination
identiﬁer. Any frames with the same FRINCID will be com-id. 4

NICpolpy: NIC Polarimetry Reducer
mfrin_H_Elias19_sky_001_o.fits

400

200
100
0

100

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

bined into one single master fringe frame. By default, it isset <FILTER> <OBJECT> <id> <OERAY>. Except for <id>,they are as in subsection 3.2. <id> is an integer that separates blank sky exposures that are separated by anotherexposure. For example, blank sky exposures before andafter 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 ﬁrst input image has no
scaling”, which is reproducible by scale to 0th=True9 .
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/) andsummary ﬁle summ(MFRIN). An example of a combined
master fringe thumbnail is shown in ﬁgure 9.
9

3.4.11. Cell 10: lv4 Planner
The important default arguments are:npr.plan_lv4(add_crrej=True, add_mfrin=True)num: 42900min: -1.9258685111999512
max: 1.3740648031234741
avg: 0.9812648296356201
med: 0.9819
std: 0.09484384208917618
madstd: 0.0853
percentiles : ['1.0000', '95.0000', '99.0000']pct: ['0.7748', '1.1247', '1.1885']
zmin: 0.4922
zmax: 1.3013
ext_lo: ['-1.9259', '-1.6863', '-1.4857']
ext_hi: ['1.3535', '1.3597', '1.3741']

300

0

11

Note that the results may be diﬀerent 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.

The ﬁnal level 4 processing is cosmic-ray removal
and fringe subtraction. Those can be turned on or
oﬀ by add crrej and add mfrin.
Here, we used
add mfrin=False. The planner plan(lv4) is made, andit contains multiple columns related to the cosmic-ray rejection algorithm as in plan(IFRIN) (see sub-subsection
3.4.8). If add mfrin=True, a column named FRINFRM isadded 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 ﬁnds any extended bright object
using SExtractor algorithm (Bertin & Arnouts 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
pixel (frin sep kw) is found. The fringe frame is thenscaled based on these pixel (sigma-clipped median is used;frin scale kw) and subtracted from the science exposure
image. Finally, cosmic-rays are removed based on the parameters speciﬁed in plan(lv4).
3.4.13. FITS Header
Finally, each FITS ﬁle 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 imagein o-/e-ray region.
• NSATPIX[O/E]: The number of saturated pixels in o/e-ray region, based on level 1 image (sub-subsection
3.4.5).
• FFTCUTWL: The minimum wavelength of Fourier pattern removal used for level 2 (sub-subsection 3.4.6)• LViFRM: The path to the level i frame related to thisﬁle.
• [DARK/FLAT/FRIN]FRM: The master dark, ﬂat, or
fringe frame used for this ﬁle.
• MASKFILE: The path to mmask ﬁle.
• MASKNPIX: The number of masked (ﬁxed) pixels inlevel 3 (sub-subsection 3.4.7).
• CRNFIX: The number of pixels ﬁxed by cosmic-ray rejection algorithm in level 4 (sub-subsection 3.4.12).
• others: YYYYMMDD, COUNTER, SETID, OERAY, PROCESS
as in subsection 3.2.

12

Y. P. Bach et al.

Also, as COMMENT and HISTORY, most of the minor logs
for this speciﬁc ﬁle are saved, with options used, timestamp, and the time taken for the step. Below is the lastpart 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_SSECHISTORY 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 ﬂat frames wereobtained on UT 2020-06-03 (total 640 frames per ﬁlter; 8
INSROT and 4 POL-AGL1 combinations and 20 FITS frames
per each combination). These frames are useful to checkif there is any region in the RoI the observer should avoid
putting the scientiﬁc targets. Also, these ﬂat frames wereused to re-calculate the gain factor (section 5).
The frames underwent these processes:
1. Reduce all frames to level 1.
2. Split o-/e-ray regions by cropping only the lit parts(ﬁgure 1).
3. Scale each frame by “median value/median value ofthe ﬁrst frame”.
The master ﬂat frame is obtained by considering all
640 frames per ﬁlter. The ﬂat frames were also groupedby diﬀerent 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 4groups per ﬁlter, respectively. The median and sample

Vol. 5,(a) Normalized flat(b) SNR

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

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 ﬂats (normalized) and SNR maps are shownin ﬁgure 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 Jband o-ray region in the mid-right part. These are visible
only in the SNR map but not in the master ﬂat map.
These doughnuts could also be found in the ﬂat frames
taken on UT 2018-05-07 using a similar technique.
An interesting feature appeared for ﬂats grouped basedon INSROT. Deﬁne the ratio map as
r1 (INSROT, POL-AGL1) =comb(INSROT, POL-AGL1)(3)
comb(any, POL-AGL1)

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 eﬀect of INSROT. An example for J-band oray region for POL-AGL1 = 0 are shown in ﬁgure 11. The
relative location of the brighter and darker part of thecommon 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 ﬂat may be diﬀerent fromthe true ﬂat at the observation by up to ∼ 2−3% at suchregions, depending on the INSROT at the exposure. For
safety, any observer is encouraged to put their scientiﬁc
object such that they are not aﬀected by these doughnuts.
Also, the vertical cut proﬁle has a gradient of ∼ ± 1%,
changing as a function of INSROT. The change in the lit
part is seen, too.

id. 4

NICpolpy: NIC Polarimetry Reducer

13

J-band, POL-AGL1 = 0.0°, o-ray.
-135°

-90°

-45°

0°

INSROT

+45°

+90°

+135°

+180°

Fig. 11. The r1 ratio map for J-band with HWP rotator angle 0◦ and o-ray region only. The numbers on top are the INSROT valuesof the group. The vertical and horizontal graphs show the cut proﬁles along red lines in the ﬂat images. The red lines in the graphs
are the cut proﬁle of ±10pixel around the red lines in the image, and black dotted lines are the minimum and maximum pixel values
within the ± 10 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 eﬀectof POL-AGL1:
comb(INSROT, POL-AGL1)
r2 (INSROT, POL-AGL1) =. (4)comb(INSROT, any)
In r2 maps, the dust doughnuts disappear. This is likely
because the dark-bright features of doughnuts cancel out
when combined over the INSROT. The combined ﬂats show
some large-scale pattern in the ﬁnal ﬂat, with ﬂuctuation< ± 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 ﬂat (ﬁgure 2), it is better toavoid the top-left part due to the over density of bad pixels
in H-band (e-ray) and K-band (e-ray). These regions maychange over time as new defects appear.
The results for diﬀerent cases (ﬁgures and videos forboth r1 and r2 ) are available in the supporting material
in Appendix 2.
5.

Gain and Readout Noise: Revisited

The gain (the conversion factor) and readout noise areimportant parameters in the error analysis and data reduction (e.g., used in fringe frame scaling, cosmic-ray rejection, and object detection algorithms). They were ﬁrst
reported in Ishiguro et al. (2011). According to them,the expected noise in the bias frame corresponds to thereadout noise divided by the gain factor, i.e., 5.4, 7.7, and
8.8 ADU for J-, H-, and K-bands, respectively. However,the sample standard deviations in the dark area in level1 images are smaller (roughly 4 ADU for all detectors).
This was the motivation for re-calculating the two detec-

tor parameters. We also brieﬂy tested if there is any hintthat these parameters are diﬀerent per pixel.
If the two ﬂat images have pixel values of F1 and F2
and two biases have B1 and B2 (all in ADU), the gain
and readout noise are(F̄1 + F̄2 ) − (B̄1 + B̄2 )g=[e/ADU](5)2
σF2 1 −F2 − σB
1 −B2
and
R=g

σB1 −B2
√
[e] .
2

(6)

Here X̄ means the average of all the pixels in the frame
X, and σX is the true standard deviation of the frameX,∑estimated from the sample standard deviation σX ≈[( i (Xi − X̄)2 ) / (N − 1)]1/2 . This method is called
Janesick’s method. The proof is given in Appendix 1.
In the case of NIC, the bias is eﬀectively 0, and thestandard deviation of the bias-diﬀerence map is obtainedby the diﬀerence map of short-exposure darks to minimize the Poisson noise eﬀect. We selected 20 level 1 darkframes with 2 s exposure time taken on UT 2019-10-22
for each of the three ﬁlters (detectors) to calculate R/g.
Next, the ﬂat 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 ﬂat frame, the bad pixels (imask,
subsection 2.1) are masked.
5.1.

R/g

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

14

Y. P. Bach et al.

Vol. 5,𝑓!.!$; & '(),& +()𝑓!.#$; & '(),& +()1%, 5%, 10% criteria lines𝑆), 𝜎,,
×𝑆,, 𝜎),

Fig. 12. AnexampleofF-distribution
probability
distribution function (pdf) with degrees of freedom(n1 − 1, n2 − 1) = (19, 19). The rejection regions of signiﬁcance level α = 0.10 for a two-tailed test are shown by two
hatched regions.

5.1.1. Pixel-to-Pixel Variability
Are the R/g values dependent on each pixel, i.e., eachpixel has diﬀerent 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 σ 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′ ′H0 : σ ij = σ i j;′ ′H1 : σ ij ̸= σ i j ,(7)

we can test the null hypothesis (H0 ) using the F-statistic(e.g., Walpole et al. 2011, subsection 4.5). Deﬁning
S 2 σ2
F := 12 × 22 ∼ Fn1 −1,n2 −1 ,
S2 σ1

(8)

where ∼ means the statistic in the left-hand side follows
the distribution in the right-hand side, and Fn1 −1,n2 −1means the F-distribution of degrees of freedom (n1 −1, n2 − 1). n1 and n2 are the number of sample used toobtain S1 and S2 , respectively. Under H0 , the value ofthe statistic f = S12 / S22 . ﬁgure 12 visually depicts how
the rejection process works: If our calculated f value isinside one of those hatched regions, our null hypothesis isrejected with signiﬁcance level α = 0.10, i.e., σ1 ̸= σ2 . If
our statistic is not in the hatched regions, we can make noconclusion under the signiﬁcance level α = 0.10. A numberthat is more frequently used is the P-value: It measuresthe degree of rarity of the observation (f value) under H0 .
Consider two cases where one obtains f = f0.20;n1 −1,n2 −1
and f = f0.02;n1 −1,n2 −1 . the P-value will be 0.4, and 0.04,
respectively, which means the latter case is rarer than theformer one. Thus, the null hypothesis is rejected in thelatter case (H0 is unlikely to be true because such a rareobservation was made), while it is not in the former case.
Because we have ∼ 106 pixels at each ﬁlter, we have
∼ 1012 possible F-statistic values (pairs). To save time,
we randomly selected ∼ 105 pairs, calculated the P-value
for each pair, and saw the distribution of the P-values.

Fig. 13. Two examples of P-value calculations and their distributions for 124, 750 random pixel pairs. The settings for
obtaining these two ﬁgures are identical, except these twosimulations samples diﬀerent random pixel pairs. The threevertical lines from left to right mean signiﬁcance s of 0.01,
0.05, and 0.10, respectively.

ﬁgure 13 shows two trials of such test on K-band detector dark frames. We found no signiﬁcant population of
the histogram at small P-values, even at multiple trials
and diﬀerent detectors. This means we cannot reject thenull hypothesis in the vast majority of cases. Therefore,
we have no clear evidence to conclude that R/g is pixeldependent.
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 thathave values outside the 5-sigma for each frame are masked.
These two dark frames are subtracted, as in ﬁgure 14.
The histogram shows the distribution is very close to anormal distribution as √one would expect, and the standard
deviation divided by 2, i.e., R/g = 3.67 ADU for thisspeciﬁc case.
Since there are total 20 dark frames, one can maken(n−1)= 190 pairs of dark frames and see how the R/g is2
distributed. The result is shown in ﬁgure 15. Regardlessof the regions we chosen, the ﬁnal R/g is well constrainedto ∼ 3.7 ADU for all three detectors, with scatter muchless than 0.1 ADU.

id. 4

NICpolpy: NIC Polarimetry Reducer

15(a) full frame with 50 pix bezels
350

3000

300

2500

250
200

2000

150

1500

100

1000

50

500

0

std/ 2 = 3.67

3500

0

100

0

20

10

0

10

20

(b) o-/e-ray regions only

Fig. 14. The diﬀerence between two random K-band darkframes (left; blue and red are −10 and +10 ADU, respectively), and the pixel value distribution histogram (right).
The ﬁtted Gaussian function is shown in the red dashed linewith its standard deviation information.

5.1.3. Quadrant-to-Quadrant Variability
Finally, we tested if the R/g depends on the quadrants. The technique, the same as ﬁgure 15, is applied
to each quadrant and checked if any clear diﬀerence could
be seen (ﬁgure 16). Although the distribution is diﬀerent
for quadrants in K-band, the diﬀerence is only ∼ 0.1ADU.
Therefore, we conclude there is no quadrant dependency
in R/g.
5.2.

Gain (The Conversion Factor) & Readout Noise

Hinted by the R/g calculation, we now assume both thegain and readout noise are constant values for all pixels
in each detector. As before, NIC has no bias, so B̄ = 0.
2
in equation (5) can be
Since R/g is a constant, σB
1 −B2
2replaced by 2 (R/g) , and(( )2 )
ijR
ij
ij f
F ∼N f ,
+2
(9)
,
g
g
for the ﬂat frame pixel value variable F ij and its measured
value f ij at the pixel (i,j). Here, N (µ,σ 2 ) is the normal
distribution with mean µ and variance σ 2 . The two terms
in the variance are the Poisson and readout noise terms.
Say the true pixel value (i.e., ﬂux) of F ij is constant.
Then one can calculate the average of the location F̄ ij
and sample standard deviation S ij to obtain
g=

F̄ ij
.
S ij − (R/g)2(10)

Here, R/g = 3.7 ADU is already obtained previously andfound to be constant.
For the assumption (F ij is constant) to hold, one has toselect ﬂat frames taken at a similar time, the same HWP
angles (POL-AGL1), and the same INSROT because we onlyused the polarimetric ﬂat for this calculation. The ﬂat
ﬂux changes during the exposures (ﬁgure 17). The largechanges after every 80 frames are due to the change in

Fig. 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 ﬁrst (Q1 )and third (Q3 ) quantiles, and the black whiskers extend
until min{min{y}, Q1 − 1.5(Q3 − Q1 )} to the bottom andmax{max{y}, Q3 + 1.5(Q3 − Q1 )} to the top, where min{y}and max{y} are the minimum and maximum of the data, respectively. Black circles show the outliers (ﬂiers) outside thewhisker.

INSROT, and the scatter within the same INSROT is due tothe change in POL-AGL1 (semi-periodic for every 4 frames)
and/or the true change in the dome ﬂat.
There are 20 consecutive dome ﬂat exposures for thesame POL-AGL1 and same INSROT for each ﬁlter. For a, a few ﬂats with the same (POL-AGL1, INSROT) wereselected, and the change in the pixel statistics was calculated. ﬁgure 18 shows two such cases. The pixel valuechanges by <∼ 1% level. Hence, for those 20 frames, it canbe assumed that F ij is constant.
The median and standard deviation were obtained foreach 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 equation (10) to obtain a 2D 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 ﬁnal gain map and the gain value distributions are shown in ﬁgure 19.

16

Y. P. Bach et al.

Vol. 5,

Fig. 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 rightquadrants, respectively.

(a)

INSROT
0 -45-90 -135 +45 +90+135 +180

(b)

o-ray

INSROT
0 -45-90 -135 +45 +90+135 +180

e-ray

(c)

(d)

(e)

(f)

Fig. 17. The variability of ﬂat 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 thesigma-clipped standard deviation (right ordinate), for the RoIs (green regions in ﬁgure 1). The INSROT on top shows the instrument
rotator angle at the exposures.

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 constantwould be conﬁrmed separately in future work. However,
if we assume the two parameters are indeed constant, we

obtain gain and readout noise values for all detectors. Intable 2, the gain and readout noise parameters from thiswork are compared with the previous work (Ishiguro et al.
2011).
In the table, the last column shows the values adopted
in NICpolpy (callable as nic.GAIN and NIC.RDNOISE).

id. 4

NICpolpy: NIC Polarimetry Reducer

17

Fig. 18. The change in ﬂat 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ﬁgure 1).

Ishiguro et al. (2011)
Parameter

R/g [ADU]g [e/ADU]R [e]This work

Adopted value(mean ± std)

in NICpolpy

3.7 (± < 0.1)

(3.7)

Detector

Motor ON

Motor OFFJ-band

4.7 ± 0.5

5.4 ± 0.5

H-band8.8 ± 0.6

7.7 ± 0.4

K-bandN/A

8.8 ± 0.6

J-band9.8 ± 0.2

9.2 ± 0.2

9.9 ± 1.0

9.9

H-band

9.5 ± 0.2

9.8 ± 0.2

9.8 ± 0.9

9.8

K-bandN/A

9.4 ± 0.2

9.5 ± 0.9

9.5

J-band46 ± 550 ± 437 ± 437

H-band84 ± 575 ± 436 ± 336

K-band(> 250)83 ± 535 ± 335

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.

The oﬃcial 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 ﬂat
diﬀerence image is a linear function of signal (ﬁgure 2 ofIshiguro et al. 2011), the Poisson noise regime. The readout noise (or the variance without light) is obtained byextrapolating it to zero signal. This extrapolation, however, works only if the pixels in the ﬂat frame used to calculate variance have nearly identical pixel sensitivity (sothat they indeed have the same true ﬂat value). Moreover,in reality, there is always the R/g term in the variance, so
extrapolation using the equation ignoring R/g may havein an erroneous estimation of R. However, we note thatthe more important value, gain (the conversion factor),
was determined correctly.
Appendix 1.

Proof of Janesick’s Method

Consider that ﬂat and the bias images are F and B,
and any operation (e.g., addition or subtraction) meansa pixel-wise operation. Also, any artifacts, bad-pixels or
cosmic-rays are ignored. If the true bias level is b in ADU,

18

Y. P. Bach et al.
which means
F1 − F2 ∼ N

Vol. 5,(( )2 )
fR
0, 2 + 2
,
g
g

and

(A5)

( )2 )
Rf
2f, 2 + 4
. (A6)gg((F1 + F2 ) − (B1 + B2 ) ∼ N

The f /g terms are the Poisson noise term, and the R/gterms are the readout noise term, respectively. Fromthese,
2
σF2 1 −F2 − σB
≈2
1 −B2

f
,
g

(A7)

and
(F̄1 + F̄2 ) − (B̄1 + B̄2 ) ≈ 2f .

(A8)This proves equation (5).
Since the σ values
are approximated by the sample standard deviation,2
σF2 1 −F2 − σB
≈ np.std(F1-F2, ddof=1)**2 1 −B2
np.std(B1-B2, ddof=1)**2 and (F̄1 + F̄2 ) − (B̄1 + B̄2 ) =np.mean((F1+F2)-(B1+B2)) in python.
Appendix 2.

Supporting Materials

The basic calibration ﬁles (master ﬂats for UT 2018-0507 and 2020-06-04 and imask), ﬂat analyses results (section 4), and the data reduction example with sample data
(section 3) are available via the GitHub service10 . The
contents are:• cal-flat 20180507-lv1/: The master ﬂats taken
on UT 2018-05-07.
• cal-flat 20200603-lv1/: The master ﬂats taken
on UT 2020-06-04.
• masks/: The imask frames (subsection 2.1).
• flat analyses: The ﬁgures and videos of ﬂat
maps and their change over INSROT or POL-AGL1.
sm1.pptx ﬁle here describes the details of how theﬁgures are generated.
• example/: Sample data, reduction code and output
log directory contents (section 3).

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

any bias image will follow a normal distribution with readout noise:()
2
B ∼ N b, (R/g) [ADU] ,(A1)
so that()
2
B1 − B2 ∼ N 0, 2 (R/g) [ADU] .

(A2)Hence,2σB= 2(R/g)2 ,
1 −B2

(A3)

so equation (6) is proven. Since the left-hand side here isapproximated by the sample standard deviation, one may
2
write σB
≈ np.std(B1-B2, ddof=1)**2 in python.
1 −B2
The (raw) ﬂat image consists of photons with dark plus
bias level. Therefore, if f is the true ﬂat level plus darkin ADU, any ﬂat will follow a Poisson distribution, which
is roughly identical to a normal distribution when pixel
values are high:(( )2 )fR
F ∼ N f + b, +
[ADU]
(A4)
g
g

References
McCully, C., et al. 2018, astropy/astroscrappy: v1.0.5,
Zenodo, https://doi.org/10.5281/zenodo.1482019
Ishiguro, M., Takahashi, J., Zenno, T., Tokimasa, N. &Kuroda, T. 2011, Annu. Rep. Nishi-Harima Astron. Obs.,
21, 13
Takahashi, J., Zenno, T., Saito, T., & Itoh, Y. 2018, Stars And
Galaxies, 1, 17
van Dokkum, P. G. 2001, PASP, 113, 1420
Takahashi, J., Itoh, Y., Matsuo, T., Oasa, Y., Bach, Y., &Ishiguro, M. 2021, A&A, 653, A99
Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393
Walpole, R., Myers, R., Myers, S., & Ye, K. 2011, Essentials
of Probability & Statistics for Engineers & Scientists (New
Jersey: Prentice Hall)
10

https://github.com/ysBach/NICpolpy sag22sm