Implementation of reduced shear in source selection for clusters

examples/cosmodc2/Cluster_pipelines/20deg2-in2p3-reduced.sub

@@ -0,0 +1,17 @@
+#!/usr/bin/bash
+#SBATCH --time=02:00:00
+#SBATCH --partition=hpc,lsst
+#SBATCH --ntasks=20
+#SBATCH --cpus-per-task=1
+#SBATCH --mem=128000
+
+module load conda
+#source /pbs/throng/lsst/software/desc/common/miniconda/setup_current_python.sh
+conda activate /sps/lsst/groups/clusters/cl_pipeline_project/conda_envs/txpipe_clp
+export HDF5_DO_MPI_FILE_SYNC=0
+export PYTHONPATH=/sps/lsst/users/ebarroso/TXPipe:$PYTHONPATH
+ceci CLClusterReducedShear-20deg2-CL.yml 
+
+#source ./conda/bin/activate
+#export HDF5_DO_MPI_FILE_SYNC=0
+#ceci examples/cosmodc2/Cluster_pipelines/pipeline-20deg2-CL-in2p3.yml

examples/cosmodc2/Cluster_pipelines/CLClusterReducedShear-20deg2-CL.yml

@@ -0,0 +1,66 @@
+#this step depends on where you run
+#for CCin2p3
+site:
+    name: cc-parallel
+    mpi_command: mpirun -n
+    max_threads: 20
+#for NERSC
+#site:
+#    name: cori-batch
+#    image: ghcr.io/lsstdesc/txpipe-dev
+
+
+#all the following steps should not depend on where you run 
+launcher:
+    name: mini
+    interval: 3.0
+modules: >
+    txpipe
+    rail.estimation.algos.bpz_lite
+
+python_paths: []
+
+stages:
+    - name: TXSourceSelectorMetadetect
+      nprocess: 20
+    - name: BPZliteInformer
+      nprocess: 1
+      aliases:
+        input: spectroscopic_catalog
+        model: photoz_model
+    - name: BPZliteEstimator
+      nprocess: 20
+      aliases:
+        model: photoz_model
+        input: shear_catalog
+        output: source_photoz_pdfs
+    - name: CLClusterBinningRedshiftRichness
+      nprocess: 1
+    - name: CLClusterShearCatalogs
+      nprocess: 1 #>1 does not work with mpi
+    - name: CLClusterEnsembleProfiles
+      nprocess: 10
+    - name: CLClusterSACC
+      nprocess: 1
+      aliases:
+          cluster_profiles: cluster_profiles
+
+output_dir: /sps/lsst/users/ebarroso/TXPipe/data/cosmodc2/outputs-20deg2-CL
+config: /sps/lsst/users/ebarroso/TXPipe/examples/cosmodc2/Cluster_pipelines/config-20deg2-CL.yml
+
+inputs:
+inputs:
+    # See README for paths to download these files
+    shear_catalog: /sps/lsst/groups/clusters/cl_pipeline_project/TXPipe_data/cosmodc2/20deg2/shear_catalog.hdf5
+    #photometry_catalog:  /sps/lsst/users/mricci/TXPipe_data/data_link/cosmodc2/20deg2/photometry_catalog.hdf5
+    fiducial_cosmology: /sps/lsst/groups/clusters/cl_pipeline_project/TXPipe_data/cosmodc2/fiducial_cosmology.yml
+    calibration_table: /sps/lsst/groups/clusters/cl_pipeline_project/TXPipe_data/cosmodc2/20deg2/sample_cosmodc2_w10year_errors.dat
+    spectroscopic_catalog: /sps/lsst/groups/clusters/cl_pipeline_project/TXPipe_data/cosmodc2/20deg2/spectroscopic_catalog.hdf5
+    cluster_catalog: /sps/lsst/groups/clusters/cl_pipeline_project/TXPipe_data/cosmodc2/20deg2/cluster_catalog.hdf5
+    #shear_tomography_catalog: ./data/example/outputs_metadetect/shear_tomography_catalog.hdf5
+    #source_photoz_pdfs: ./data/example/inputs/photoz_pdfs.hdf5
+
+resume: true
+log_dir: /sps/lsst/users/ebarroso/TXPipe/data/cosmodc2/logs
+pipeline_log: /sps/lsst/users/ebarroso/TXPipe/data/cosmodc2/log_20deg2.txt
+

examples/cosmodc2/Cluster_pipelines/config-20deg2-CL.yml

@@ -14,8 +14,8 @@ BPZliteInformer:
     zmin: 0.0
     zmax: 3.0
     nzbins: 301
-    columns_file: ./data/bpz_riz.columns
-    data_path: ./data/example/rail-bpz-inputs
+    columns_file: /sps/lsst/groups/clusters/cl_pipeline_project/TXPipe_data/bpz_riz.columns
+    data_path: /sps/lsst/groups/clusters/cl_pipeline_project/TXPipe_data/example/rail-bpz-inputs
     spectra_file: CWWSB4.list
     prior_band: i
     ref_band: i
@@ -37,12 +37,12 @@ BPZliteEstimator:
     zmax: 3.0
     dz: 0.01
     nzbins: 301
-    data_path: ./data/example/rail-bpz-inputs
+    data_path: /sps/lsst/groups/clusters/cl_pipeline_project/TXPipe_data/example/rail-bpz-inputs
     bands: [mag_r, mag_i, mag_z]
     err_bands: [mag_err_r, mag_err_i, mag_err_z]
     hdf5_groupname: shear/00
     nondetect_val: .inf
-    columns_file: ./data/bpz_riz.columns
+    columns_file: /sps/lsst/groups/clusters/cl_pipeline_project/TXPipe_data/bpz_riz.columns
     spectra_file: CWWSB4.list
     ref_band: mag_i
     prior_file: hdfn_gen
@@ -66,28 +66,27 @@ CLClusterShearCatalogs:
     chunk_rows: 100_000  # rows to read at once from source cat
     max_radius: 10.  # Mpc
     delta_z: 0.2   # redshift buffer
-    redshift_cut_criterion: ztrue
-    redshift_weight_criterion: ztrue
+    redshift_cut_criterion: zmode
+    redshift_weight_criterion: zmode
     redshift_cut_criterion_pdf_fraction: 0.9
     subtract_mean_shear: false
     coordinate_system: celestial #euclidean or celestial
     use_true_shear: false
-
+    delta_sigma: false
 
 CLClusterEnsembleProfiles:
     #radial bin definition
     r_min: 0.5  #in Mpc
     r_max: 10.  #in Mpc
-    nbins: 10  # number of bins
+    nbins: 3  # number of bins
     #type of profile
     delta_sigma_profile: true
     shear_profile: false
     magnification_profile: false
     coordinate_system: celestial #euclidean or celestial
-
+    units: 'arcmin'
+    angle_arcmin_min: 5
+    angle_arcmin_max: 30
 CLClusterSACC:
     survey_name: 'cosmodc2-20deg2-CL'
-    area: 440.0
-
-
-
+    area: 20.0

examples/cosmodc2/Cluster_pipelines/pipeline-20deg2-CL-in2p3.yml

@@ -49,20 +49,20 @@ stages:
 
 
 
-output_dir: ./data/cosmodc2/outputs-20deg2-CL
-config: ./examples/cosmodc2/Cluster_pipelines/config-20deg2-CL.yml
+output_dir: /sps/lsst/users/ebarroso/TXPipe/data/cosmodc2/outputs-20deg2-CL
+config: /sps/lsst/users/ebarroso/TXPipe/examples/cosmodc2/Cluster_pipelines/config-20deg2-CL.yml
 
 inputs:
     # See README for paths to download these files
-    shear_catalog: /sps/lsst/users/mricci/TXPipe_data/data_link/cosmodc2/20deg2/shear_catalog.hdf5
+    shear_catalog: /sps/lsst/groups/clusters/cl_pipeline_project/TXPipe_data/cosmodc2/20deg2/shear_catalog.hdf5
     #photometry_catalog:  /sps/lsst/users/mricci/TXPipe_data/data_link/cosmodc2/20deg2/photometry_catalog.hdf5
-    fiducial_cosmology: ./data/fiducial_cosmology.yml
-    calibration_table: /sps/lsst/users/mricci/TXPipe_data/data_link/cosmodc2/20deg2/sample_cosmodc2_w10year_errors.dat
-    spectroscopic_catalog: /sps/lsst/users/mricci/TXPipe_data/data_link/cosmodc2/20deg2/spectroscopic_catalog.hdf5
-    cluster_catalog: /sps/lsst/users/mricci/TXPipe_data/data_link/cosmodc2/20deg2/cluster_catalog.hdf5
+    fiducial_cosmology: /sps/lsst/users/ebarroso/TXPipe/data/fiducial_cosmology.yml
+    calibration_table: /sps/lsst/groups/clusters/cl_pipeline_project/TXPipe_data/cosmodc2/20deg2/sample_cosmodc2_w10year_errors.dat
+    spectroscopic_catalog: /sps/lsst/groups/clusters/cl_pipeline_project/TXPipe_data/cosmodc2/20deg2/spectroscopic_catalog.hdf5
+    cluster_catalog: /sps/lsst/groups/clusters/cl_pipeline_project/TXPipe_data/cosmodc2/20deg2/cluster_catalog.hdf5
     #shear_tomography_catalog: ./data/example/outputs_metadetect/shear_tomography_catalog.hdf5
     #source_photoz_pdfs: ./data/example/inputs/photoz_pdfs.hdf5
 resume: true
-log_dir: ./data/cosmodc2/logs
-pipeline_log: ./data/cosmodc2/log_20deg2.txt
+log_dir: /sps/lsst/users/ebarroso/TXPipe/data/cosmodc2/logs
+pipeline_log: /sps/lsst/users/ebarroso/TXPipe/data/cosmodc2/log_20deg2.txt
 

txpipe/extensions/cluster_counts/convert_to_sacc.py

@@ -14,9 +14,10 @@ class CLClusterSACC(PipelineStage):
     ]
 
     config_options = {
-        #radial bin definition
-        "r_min" : 0.2, #in Mpc
-        "r_max" : 5.0, #in Mpc
+        #distance bin definition
+        "d_min" : 0.5, #in Mpc or acrmin
+        "d_max" : 5.0, #in Mpc or arcmin
+        "cov_component": "tan_jk",
     }
 
     def run(self):
@@ -28,11 +29,12 @@ def run(self):
         survey_name = my_configs['survey_name']
         area = my_configs['area']
         output_filename = self.get_output("cluster_sacc_catalog", final_name=True)
-
+        first_bin = next(iter(data.values()))
+        profile_type = first_bin["profile_type"]        
         sacc_obj = sacc.Sacc()
         self.add_tracers(sacc_obj, data, survey_name, area)
         self.add_counts_data(sacc_obj, data, survey_name)
-        self.add_deltasigma_data(sacc_obj, data, survey_name)
+        self.add_shear_data(sacc_obj, data, survey_name, profile_type)
         self.add_covariance_data(sacc_obj, data)
 
         sacc_obj.to_canonical_order()
@@ -66,12 +68,12 @@ def get_bins(self, data: dict):
             bin_z_dict[bin_z] = z_edges
             bin_rich_dict[bin_rich] = rich_edges
 
-        radius_centers = np.array(data['bin_zbin_0_richbin_0']['clmm_cluster_ensemble'].stacked_data['radius'])
-        rmin = self.config_options['r_min']
-        rmax = self.config_options['r_max'] 
-        radius_edges = np.logspace(np.log10(rmin), np.log10(rmax), len(radius_centers) + 1)
-        for i in range(len(radius_edges) - 1):
-            bin_radius_dict[f'radius_{i}'] = (radius_edges[i], radius_edges[i+1], radius_centers[i])
+        distance_centers = np.array(data['bin_zbin_0_richbin_0']['clmm_cluster_ensemble'].stacked_data['radius'])
+        dmin = self.config_options['d_min']
+        dmax = self.config_options['d_max'] 
+        distance_edges = np.logspace(np.log10(dmin), np.log10(dmax), len(distance_centers) + 1)
+        for i in range(len(distance_edges) - 1):
+            bin_radius_dict[f'radius_{i}'] = (distance_edges[i], distance_edges[i+1], distance_centers[i])
 
         return bin_z_dict, bin_rich_dict, bin_radius_dict
 
@@ -103,12 +105,16 @@ def add_counts_data(self, sacc_obj, data: dict, survey_name: str):
             bin_z, bin_rich = self.transform_bin_string(bin_comb)
             sacc_obj.add_data_point(cluster_count, (survey_name, bin_rich, bin_z), int(bin_data['n_cl']))
 
-    def add_deltasigma_data(self, sacc_obj, data: dict, survey_name: str):
+    def add_shear_data(self, sacc_obj, data: dict, survey_name: str, profile_type: str):
         """
         Adds cluster shear (delta sigma) data to the SACC object.
         """
         import sacc
-        cluster_shear = sacc.standard_types.cluster_shear
+        print(profile_type)
+        if profile_type == 'reduced_shear':
+            cluster_shear = sacc.standard_types.cluster_shear
+        elif profile_type == 'delta_sigma':
+            cluster_shear = sacc.standard_types.cluster_delta_sigma    
         _, _, radius_bins = self.get_bins(data)
 
         for bin_comb, bin_data in data.items():
@@ -119,17 +125,20 @@ def add_deltasigma_data(self, sacc_obj, data: dict, survey_name: str):
 
     def add_covariance_data(self, sacc_obj, data: dict):
         """
-        Adds covariance data to the SACC object.
+        Adds full block-diagonal covariance data to the SACC object.
         """
         import sacc
+        from scipy.linalg import block_diag
+        cov_str = self.config["cov_component"]
         cluster_count = sacc.standard_types.cluster_counts
-        counts_points = np.array(sacc_obj.get_data_points(cluster_count))
-        counts_cov = np.array([point.value for point in counts_points])
+        counts_points = sacc_obj.get_data_points(cluster_count)
+
+        counts_values = np.array([point.value for point in counts_points])
+        counts_cov_block = np.diag(counts_values) 
 
-        deltasigma_cov = [
-            bin_data['clmm_cluster_ensemble'].cov['tan_sc'].diagonal()
+        ds_blocks = [
+            bin_data['clmm_cluster_ensemble'].cov[cov_str]
             for bin_data in data.values()
         ]
-
-        diag_cov_vector = np.concatenate([counts_cov.flatten(), np.array(deltasigma_cov).flatten()])
-        sacc_obj.add_covariance(np.diag(diag_cov_vector))
+        full_cov = block_diag(counts_cov_block, *ds_blocks)
+        sacc_obj.add_covariance(full_cov)