Add cutoff for notch filter

notebooks/single_baseline_postprocessing_and_pspec.ipynb

@@ -152,6 +152,7 @@
     "TARGET_AVERAGING_TIME = float(os.environ.get(\"TARGET_AVERAGING_TIME\", 300.0)) # coherent integration time in seconds. Actual time might be less to so that all interleaves have the same number of samples averaged\n",
     "USE_CORR_MATRIX = os.environ.get(\"USE_CORR_MATRIX\", \"TRUE\").upper() == \"TRUE\" # Whether or not to use the correlation matrix to do the noise statistics correction factor\n",
     "CORR_MATRIX_FREQ_DECIMATION = int(os.environ.get(\"CORR_MATRIX_FREQ_DECIMATION\", 10)) # How much to decimate the frequency axis when getting the noise statistic correction factor\n",
+    "CORR_MATRIX_NOTCH_CUTOFF = float(os.environ.get(\"CORR_MATRIX_NOTCH_CUTOFF\", 30.)) # Maximum EW-projection in meters for including the notch filter in error covariance calculation\n",
     "\n",
     "# Power spectrum settings\n",
     "EFIELD_HEALPIX_BEAM_FILE = os.environ.get(\"EFIELD_HEALPIX_BEAM_FILE\", \"/lustre/aoc/projects/hera/h6c-analysis/IDR2/beams/NF_HERA_Vivaldi_efield_beam_healpix.fits\")\n",
@@ -1490,20 +1491,26 @@
     "                main_lobe_frop = get_frop_wrapper(pol=pol, stream_ind=stream_ind, band_ind=spw,\n",
     "                                                  dlst=dlst, bl_vec=bl_vec[cross_antpairpol],\n",
     "                                                  times=times)\n",
-    "                xtalk_frop = get_frop_wrapper(pol=pol, stream_ind=stream_ind, band_ind=spw,\n",
-    "                                              dlst=dlst, bl_vec=bl_vec[cross_antpairpol],\n",
-    "                                              times=times, rephase=False, mode=\"xtalk\", coherent_avg=True,\n",
-    "                                              t_avg=times[1] - times[0])\n",
     "                \n",
+    "                # Calculate whether the EW-baseline projection is short enough to care about the notch\n",
+    "                EW_proj = data.antpos[ANTPAIR[0]][0] - data.antpos[ANTPAIR[1]][0]\n",
+    "                if np.abs(EW_proj) < CORR_MATRIX_NOTCH_CUTOFF:\n",
+    "                    xtalk_frop = get_frop_wrapper(pol=pol, stream_ind=stream_ind, band_ind=spw,\n",
+    "                                                  dlst=dlst, bl_vec=bl_vec[cross_antpairpol],\n",
+    "                                                  times=times, rephase=False, mode=\"xtalk\", coherent_avg=True,\n",
+    "                                                  t_avg=times[1] - times[0])\n",
+    "\n",
     "\n",
-    "                frop = np.zeros_like(main_lobe_frop)\n",
-    "                Nfreqs_band = freq_slice.stop - freq_slice.start\n",
-    "                Nfreqs_calc = int(np.ceil(Nfreqs_band / CORR_MATRIX_FREQ_DECIMATION))\n",
+    "                    frop = np.zeros_like(main_lobe_frop)\n",
+    "                    Nfreqs_band = freq_slice.stop - freq_slice.start\n",
+    "                    Nfreqs_calc = int(np.ceil(Nfreqs_band / CORR_MATRIX_FREQ_DECIMATION))\n",
     "\n",
-    "                for freq_ind in range(Nfreqs_calc):\n",
-    "                    frop[:, :, freq_ind] = np.tensordot(main_lobe_frop[:, :, freq_ind],\n",
-    "                                                        xtalk_frop[:, :, freq_ind],\n",
-    "                                                        axes=1)\n",
+    "                    for freq_ind in range(Nfreqs_calc):\n",
+    "                        frop[:, :, freq_ind] = np.tensordot(main_lobe_frop[:, :, freq_ind],\n",
+    "                                                            xtalk_frop[:, :, freq_ind],\n",
+    "                                                            axes=1)\n",
+    "                else:\n",
+    "                    frop = main_lobe_frop\n",
     "                cov_here += frf.get_FRF_cov(frop, var)  # computes covariance for pI by adding ee and nn\n",
     "            if hd.pol_convention == \"avg\":\n",
     "                cov_here /= 4\n",