Lf calibration

txpipe/calibrate.py

@@ -127,15 +127,19 @@ def run(self):
                 # Cut down the data to just this selection for output
                 d = {name: data[name][w] for name in output_cols}
 
-                # Calibrate the shear columns
+                # Calibrate the shear columns 
                 if cat_type=='hsc':
                     d["g1"], d["g2"] = cal.apply(
                         d["g1"], d["g2"], d["c1"], d["c2"]
                     )
-                else:
+                elif cat_type=='lensfit':
                     d["g1"], d["g2"] = cal.apply(
-                        d["g1"], d["g2"], subtract_mean=subtract_mean_shear
+                        d["dec"],d["g1"], d["g2"], subtract_mean=subtract_mean_shear
                     )
+                else:
+                    d["g1"], d["g2"] = cal.apply(
+                    d["g1"], d["g2"]
+                )
 
                 # Write output, keeping track of sizes
                 splitter.write_bin(d, b)

txpipe/diagnostics.py

@@ -5,8 +5,6 @@
 from .utils.calibration_tools import (
     calculate_selection_response,
     calculate_shear_response,
-    apply_metacal_response,
-    apply_lensfit_calibration,
     MeanShearInBins,
     read_shear_catalog_type,
     metadetect_variants,
@@ -62,6 +60,7 @@ class TXSourceDiagnosticPlots(PipelineStage):
         "T_max": 4.0,
         "s2n_min": 10,
         "s2n_max": 300,
+        "psf_conv": False,
         "bands": "riz",
     }
 
@@ -140,6 +139,7 @@ def run(self):
         else:
 
             shear_cols = [
+                "dec",
                 "psf_g1",
                 "psf_g2",
                 "g1",
@@ -149,9 +149,6 @@ def run(self):
                 "T",
                 "weight",
                 "m",
-                "sigma_e",
-                "c1",
-                "c2",
             ] + [f"{shear_prefix}mag_{b}" for b in self.config["bands"]]
 
         shear_tomo_cols = ["bin"]
@@ -316,11 +313,14 @@ def plot_psf_size_shear(self):
         psf_prefix = self.config["psf_prefix"]
         delta_gamma = self.config["delta_gamma"]
         nbins = self.config["nbins"]
+
         psfT_min = self.config["psfT_min"]
         psfT_max = self.config["psfT_max"]
-        
+
         with self.open_input("shear_catalog") as c:
             col = c[f"shear/{psf_prefix}T_mean"][:]
+            if ((self.config["shear_catalog_type"] == 'lensfit') & (self.config['psf_conv']==True)):
+                col = col * 0.214**2
             psfT = col[(col > psfT_min) & (col < psfT_max)]
             psf_T_edges = self.BinEdges(psfT,nbins)
             del psfT
@@ -331,6 +331,7 @@ def plot_psf_size_shear(self):
             delta_gamma,
             cut_source_bin=True,
             shear_catalog_type=self.config["shear_catalog_type"],
+            psf_conv = self.config['psf_conv']
         )
 
         while True:
@@ -457,11 +458,14 @@ def plot_size_shear(self):
         shear_prefix = self.config["shear_prefix"]
         delta_gamma = self.config["delta_gamma"]
         nbins = self.config["nbins"]
+
         T_min = self.config["T_min"]
         T_max = self.config["T_max"]
 
         with self.open_input("shear_catalog") as c:
             col = c[f"shear/{shear_prefix}T"][:]
+            if ((self.config["shear_catalog_type"] == 'lensfit') & (self.config['psf_conv']==True)):
+                col = col * 0.214**2
             T = col[(col > T_min) & (col < T_max)]
             T_edges = self.BinEdges(T,nbins)
             del T
@@ -472,6 +476,7 @@ def plot_size_shear(self):
             delta_gamma,
             cut_source_bin=True,
             shear_catalog_type=self.config["shear_catalog_type"],
+            psf_conv = self.config['psf_conv']
         )
 
         while True:
@@ -635,6 +640,7 @@ def plot_g_histogram(self):
                 g2 = data["00/g2"]
                 w = data["00/weight"]
             elif cat_type == "lensfit":
+                dec = data["dec"]
                 g1 = data["g1"]
                 g2 = data["g2"]
                 w = data["weight"]
@@ -645,8 +651,11 @@ def plot_g_histogram(self):
                 c2 = data['c2']
                 w = data["weight"]
 
-            if cat_type=='metacal' or cat_type=='metadetect' or cat_type=='lensfit':
+            if cat_type=='metacal' or cat_type=='metadetect':
                 g1, g2 = cal.apply(g1,g2)
+
+            elif cat_type=='lensfit':
+                g1, g2 = cal.apply(dec, g1,g2)
             else:
                 g1, g2 = cal.apply(g1,g2,c1,c2)
 

txpipe/maps.py

@@ -2,7 +2,7 @@
 from .data_types import TomographyCatalog, MapsFile, HDFFile, ShearCatalog
 import numpy as np
 from .utils import unique_list, choose_pixelization, rename_iterated
-from .utils.calibration_tools import read_shear_catalog_type, apply_metacal_response
+from .utils.calibration_tools import read_shear_catalog_type
 from .utils.calibrators import Calibrator
 from .mapping import ShearMapper, LensMapper, FlagMapper
 

txpipe/metadata.py

@@ -66,9 +66,12 @@ def copy_source_metadata(self, meta_file, metadata, area, area_sq_arcmin):
                 copy(shear_tomo_file, "response", "tracers", "R_2d", meta_file, metadata)
             elif shear_catalog_type == "lensfit":
                 copy(shear_tomo_file, "response", "tracers", "K", meta_file, metadata)
-                copy(shear_tomo_file, "response", "tracers", "C", meta_file, metadata)
+                copy(shear_tomo_file, "response", "tracers", "C_N", meta_file, metadata)
+                copy(shear_tomo_file, "response", "tracers", "C_S", meta_file, metadata)
                 copy(shear_tomo_file, "response", "tracers", "K_2d", meta_file, metadata)
-                copy(shear_tomo_file, "response", "tracers", "C_2d", meta_file, metadata)
+                copy(shear_tomo_file, "response", "tracers", "C_2d_N", meta_file, metadata)
+                copy(shear_tomo_file, "response", "tracers", "C_2d_S", meta_file, metadata)
+
             elif shear_catalog_type == "hsc":
                 copy(shear_tomo_file, "response", "tracers", "R", meta_file, metadata)
                 copy(shear_tomo_file, "response", "tracers", "K", meta_file, metadata)

txpipe/psf_diagnostics.py

@@ -7,12 +7,15 @@
     TomographyCatalog,
     RandomsCatalog,
     YamlFile,
+    TextFile
+
 )
 from parallel_statistics import ParallelHistogram, ParallelMeanVariance
 import numpy as np
 from .utils.calibration_tools import read_shear_catalog_type
-from .utils.calibration_tools import apply_metacal_response, apply_lensfit_calibration
 from .plotting import manual_step_histogram
+from .utils.calibrators import Calibrator
+
 
 STAR_PSF_USED = 0
 STAR_PSF_RESERVED = 1
@@ -483,6 +486,7 @@ def load_stars(self):
     def load_galaxies(self):
         # Columns we need from the shear catalog
         cat_type = read_shear_catalog_type(self)
+        _, cal = Calibrator.load(self.get_input("shear_tomography_catalog"))
 
         # Load tomography data
         with self.open_input("shear_tomography_catalog") as f:
@@ -526,23 +530,19 @@ def load_galaxies(self):
         # Change shear convention 
         if self.config["flip_g2"]:
             g2 *= -1
-
         # Apply calibration factor
         if cat_type == "metacal" or cat_type == "metadetect":
             print("Applying metacal/metadetect response")
-            g1, g2 = apply_metacal_response(R_total_2d, 0.0, g1, g2)
+            g1, g2 = cal.apply(g1, g2)
 
         elif cat_type == "lensfit":
             print("Applying lensfit calibration")
-            g1, g2, weight, _ = apply_lensfit_calibration(g1      = g1,
-                                                          g2      = g2,
-                                                          weight  = weight,
-                                                          sigma_e = sigma_e,
-                                                          m       = m
-                                                          )
+            g1, g2 = cal.apply(dec, g1,g2)
+
         else:
             print("Shear calibration type not recognized.")
-
+        g1 = g1 
+        g2 = g2 
         return ra, dec, g1, g2, weight
 
     def tau_plots(self, tau_stats):
@@ -603,7 +603,8 @@ class TXRoweStatistics(PipelineStage):
 
     name = "TXRoweStatistics"
     parallel = False
-    inputs = [("star_catalog", HDFFile)]
+    inputs = [("star_catalog", HDFFile),
+             ("patch_centers", TextFile)]
     outputs = [
         ("rowe134", PNGFile),
         ("rowe25", PNGFile),
@@ -618,6 +619,8 @@ class TXRoweStatistics(PipelineStage):
         "bin_slop": 0.01,
         "sep_units": "arcmin",
         "psf_size_units": "sigma",
+        "star_type": 'PSF-reserved',
+        "var_method": 'bootstrap'
     }
 
     def run(self):
@@ -632,6 +635,8 @@ def run(self):
         rowe_stats = {}
         for t in STAR_TYPES:
             s = star_type == t
+            if STAR_TYPE_NAMES[t] != self.config.star_type:
+                continue
             rowe_stats[0, t] = self.compute_rowe(0, s, ra, dec, e_mod, e_mod)
             rowe_stats[1, t] = self.compute_rowe(1, s, ra, dec, de_psf, de_psf)
             rowe_stats[2, t] = self.compute_rowe(2, s, ra, dec, de_psf, e_mod)
@@ -660,8 +665,8 @@ def load_stars(self):
                     "measured_T"
                 ][:] ** 2
 
-            e_psf = np.array((e1psf, e2psf))
-            e_mod = np.array((e1mod,e2mod))
+            e_psf  = np.array((e1psf, e2psf))
+            e_mod  = np.array((e1mod,e2mod))
             de_psf = np.array((de1, de2))
 
             star_type = load_star_type(g)
@@ -680,10 +685,12 @@ def compute_rowe(self, i, s, ra, dec, q1, q2):
 
         corr = treecorr.GGCorrelation(self.config)
         cat1 = treecorr.Catalog(
-            ra=ra, dec=dec, g1=q1[0], g2=q1[1], ra_units="deg", dec_units="deg"
+            ra=ra, dec=dec, g1=q1[0], g2=q1[1], ra_units="deg", dec_units="deg",
+            patch_centers=self.get_input("patch_centers")
         )
         cat2 = treecorr.Catalog(
-            ra=ra, dec=dec, g1=q2[0], g2=q2[1], ra_units="deg", dec_units="deg"
+            ra=ra, dec=dec, g1=q2[0], g2=q2[1], ra_units="deg", dec_units="deg",
+            patch_centers=self.get_input("patch_centers")
         )
         corr.process(cat1, cat2)
         return corr.meanr, corr.xip, corr.varxip**0.5
@@ -695,6 +702,8 @@ def rowe_plots(self, rowe_stats):
 
         f = self.open_output("rowe0",wrapper=True,figsize=(10,6*len(STAR_TYPES)))
         for s in STAR_TYPES:
+            if STAR_TYPE_NAMES[s] != self.config.star_type:
+                continue
             ax = plt.subplot(len(STAR_TYPES), 1, s + 1)
 
             for j,i in enumerate([0]):
@@ -721,6 +730,8 @@ def rowe_plots(self, rowe_stats):
 
         f = self.open_output("rowe134", wrapper=True, figsize=(10, 6 * len(STAR_TYPES)))
         for s in STAR_TYPES:
+            if STAR_TYPE_NAMES[s] != self.config.star_type:
+                continue
             ax = plt.subplot(len(STAR_TYPES), 1, s + 1)
 
             for j, i in enumerate([1, 3, 4]):
@@ -747,6 +758,8 @@ def rowe_plots(self, rowe_stats):
 
         f = self.open_output("rowe25", wrapper=True, figsize=(10, 6 * len(STAR_TYPES)))
         for s in STAR_TYPES:
+            if STAR_TYPE_NAMES[s] != self.config.star_type:
+                continue
             ax = plt.subplot(len(STAR_TYPES), 1, s + 1)
             for j, i in enumerate([2, 5]): 
                 theta, xi, err = rowe_stats[i, s]
@@ -775,6 +788,8 @@ def save_stats(self, rowe_stats):
         g = f.create_group("rowe_statistics")
         for i in 0, 1, 2, 3, 4, 5:
             for s in STAR_TYPES:
+                if STAR_TYPE_NAMES[s] != self.config.star_type:
+                    continue
                 theta, xi, err = rowe_stats[i, s]
                 name = STAR_TYPE_NAMES[s]
                 h = g.create_group(f"rowe_{i}_{name}")
@@ -871,6 +886,7 @@ def load_galaxies(self):
         # Columns we need from the shear catalog
         # TODO: not sure of an application where we would want to use true shear but can be added
         cat_type = read_shear_catalog_type(self)
+        _, cal = Calibrator.load(self.get_input("shear_tomography_catalog"))
 
         # load tomography data
         with self.open_input("shear_tomography_catalog") as f:
@@ -913,12 +929,10 @@ def load_galaxies(self):
 
         if cat_type == "metacal" or cat_type == "metadetect":
             # We use S=0 here because we have already included it in R_total
-            g1, g2 = apply_metacal_response(R_total_2d, 0.0, g1, g2)
+            g1, g2 = cal.apply(g1,g2)
 
         elif cat_type == "lensfit":
-            g1, g2, weight, _ = apply_lensfit_calibration(
-                g1=g1, g2=g2, weight=weight, sigma_e=sigma_e, m=m
-            )
+            g1, g2  = cal.apply(dec, g1,g2)
         else:
             print("Shear calibration type not recognized.")
 

txpipe/source_selector.py

@@ -8,7 +8,7 @@
     TextFile,
 )
 from .utils import SourceNumberDensityStats, rename_iterated
-from .utils.calibration_tools import read_shear_catalog_type, apply_metacal_response
+from .utils.calibration_tools import read_shear_catalog_type
 from .utils.calibration_tools import (
     metacal_variants,
     metadetect_variants,
@@ -709,7 +709,6 @@ def compute_output_stats(self, calculator, mean, variance):
         # Like metacal, N_eff = N for metadetect
         return BinStats(N, Neff, mean_e, sigma_e, calibrator)
 
-
 class TXSourceSelectorLensfit(TXSourceSelectorBase):
     """
     Source selection and tomography for lensfit catalogs
@@ -733,6 +732,7 @@ def data_iterator(self):
         chunk_rows = self.config["chunk_rows"]
         bands = self.config["bands"]
         shear_cols = [
+            "dec",
             "psf_T_mean",
             "weight",
             "flags",
@@ -766,20 +766,24 @@ def setup_output(self):
         # This call to the super-class method defined above sets up most of the output
         # here, so the rest of this method only does things specific to this
         # calibration scheme
+        print("set up the mf lensfit output")
         outfile = super().setup_output()
         n = outfile["tomography/bin"].size
         nbin_source = outfile["tomography/counts"].size
         group = outfile.create_group("response")
         group.create_dataset("K", (nbin_source,), dtype="f")
-        group.create_dataset("C", (nbin_source, 2), dtype="f")
+        group.create_dataset("C_N", (nbin_source, 2), dtype="f")
+        group.create_dataset("C_S", (nbin_source, 2), dtype="f")
         group.create_dataset("K_2d", (1,), dtype="f")
-        group.create_dataset("C_2d", (2), dtype="f")
+        group.create_dataset("C_2d_N", (2), dtype="f")
+        group.create_dataset("C_2d_S", (2), dtype="f")
+
         return outfile
 
     def compute_output_stats(self, calculator, mean, variance):
-        K, C, N, Neff = calculator.collect(self.comm, allgather=True)
-        calibrator = LensfitCalibrator(K, C)
-        mean_e = C.copy()
+        K, C_N,C_S, N, Neff = calculator.collect(self.comm, allgather=True)
+        calibrator = LensfitCalibrator(K, C_N, C_S)
+        mean_e = (C_N+C_S)/2
         sigma_e = np.sqrt((0.5 * (variance[0] + variance[1]))) / (1 + K)
 
         return BinStats(N, Neff, mean_e, sigma_e, calibrator)