From 5f468854f7fce7a802eba4763bd6b49c59999bbc Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Mon, 6 May 2024 12:33:18 -0700
Subject: [PATCH 01/50] hyperbolic analysis implementation

---
 .../Code/RIFT/lalsimutils.py                  |  2 +
 ...te_event_parameter_pipeline_BasicIteration |  7 +++
 .../Code/bin/helper_LDG_Events.py             | 22 ++++++++
 .../integrate_likelihood_extrinsic_batchmode  |  8 ++-
 .../Code/bin/plot_posterior_corner.py         |  4 ++
 .../Code/bin/util_CleanILE.py                 |  8 ++-
 ...ctIntrinsicPosterior_GenericCoordinates.py | 31 ++++++++++-
 .../Code/bin/util_ILEdagPostprocess.sh        |  6 ++
 .../Code/bin/util_RIFT_pseudo_pipe.py         | 55 ++++++++++++++++++-
 .../Code/bin/util_SimInspiralToCoinc.py       |  2 +-
 10 files changed, 140 insertions(+), 5 deletions(-)

diff --git a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
index d69e44908..5e0e533c3 100644
--- a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
+++ b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
@@ -33,6 +33,8 @@
     import EOBRun_module
     has_external_teobresum=True
  except:
+    print('Failed to import EOBRun_module')
+    sys.exit(1)
     has_external_teobresum=False
     True; # print(" - no EOBRun (TEOBResumS) - ")
 info_use_resum_polarizations = False
diff --git a/MonteCarloMarginalizeCode/Code/bin/create_event_parameter_pipeline_BasicIteration b/MonteCarloMarginalizeCode/Code/bin/create_event_parameter_pipeline_BasicIteration
index 4d9af0b86..c55115128 100755
--- a/MonteCarloMarginalizeCode/Code/bin/create_event_parameter_pipeline_BasicIteration
+++ b/MonteCarloMarginalizeCode/Code/bin/create_event_parameter_pipeline_BasicIteration
@@ -211,6 +211,7 @@ parser.add_argument("--ile-runtime-max-minutes",default=None,type=int,help="If n
 parser.add_argument("--cip-exe",default=None,help="filename of CIP or equivalent executable. Will default to `which util_ConstructIntrinsicPosterior_GenericCoordinates` in low-level code")
 parser.add_argument("--cip-exe-G",default=None,help="filename of CIP or equivalent executable, as ALTERNATE CIP used when a 'G' iteration is requested ")
 parser.add_argument("--use-eccentricity",default=False,action='store_true')
+parser.add_argument("--use-hyperbolic",default=False,action='store_true')
 parser.add_argument("--use-tabular-eos-file",default=False,action='store_true')
 parser.add_argument("--test-exe",default=None,help="filename of test code or equivalent executable.  Must have a --test-output <fname> argument.  Used for convergence testing or other termination. NOT ACTIVE; see 'convergence_test_samples.py' for example")
 parser.add_argument("--plot-exe",default=None,help="filename of plot code or equivalent executable.  Will default to `which plot_posterior_corner.py`.  Default is to plot last set of samples")
@@ -769,6 +770,8 @@ if  (opts.last_iteration_extrinsic):
 con_arg_str = ''
 if opts.use_eccentricity:
     con_arg_str += " --eccentricity "
+if opts.use_hyperbolic:
+    con_arg_str += " --hyperbolic "
 if opts.use_tabular_eos_file:
     con_arg_str += " --tabular-eos-file "
 con_job, con_job_name = dag_utils.write_consolidate_sub_simple(tag='join',log_dir=None,arg_str=con_arg_str,base=opts.working_directory+"/iteration_$(macroiteration)_ile", target=opts.working_directory+'/consolidated_$(macroiteration)',universe=local_worker_universe,no_grid=no_worker_grid)
@@ -789,6 +792,8 @@ con_job.write_sub_file()
 unify_arg_str = ''
 if opts.use_eccentricity:
     unify_arg_str += " --eccentricity "
+if opts.use_hyperbolic:
+    unify_arg_str += " --hyperbolic "
 if opts.use_tabular_eos_file:
     unify_arg_str += " --tabular-eos-file "
 unify_job, unify_job_name = dag_utils.write_unify_sub_simple(tag='unify',log_dir='',arg_str=unify_arg_str, base=opts.working_directory, target=opts.working_directory+'/all.net',universe=local_worker_universe,no_grid=no_worker_grid)
@@ -1370,6 +1375,8 @@ for it in np.arange(it_start,opts.n_iterations):
         cmd += " --general-request-disk {} --ile-request-disk {} ".format(opts.general_request_disk,opts.ile_request_disk)
         if opts.use_eccentricity:
             cmd += " --use-eccentricity "
+        if opts.use_hyperbolic:
+            cmd += " --use-hyperbolic "
         if opts.cip_explode_jobs:
             cmd += " --cip-explode-jobs {} ".format(opts.cip_explode_jobs)
         if opts.cip_explode_jobs_dag:   # should be deafult
diff --git a/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py b/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
index 0e9d4d8f0..ac80f5994 100755
--- a/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
+++ b/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
@@ -198,6 +198,11 @@ def get_observing_run(t):
 parser.add_argument("--assume-matter-but-primary-bh",action='store_true',help="If present, the code will add options necessary to manage tidal arguments for the smaller body ONLY. (Usually pointless)")
 parser.add_argument("--internal-tabular-eos-file",type=str,default=None,help="Tabular file of EOS to use.  The default prior will be UNIFORM in this table!. NOT YET IMPLEMENTED (initial grids, etc)")
 parser.add_argument("--assume-eccentric",action='store_true',help="If present, the code will add options necessary to manage eccentric arguments. The proposed fit strategy and initial grid will allow for eccentricity")
+parser.add_argument("--assume-hyperbolic",action='store_true',help="If present, the code will add options necessary to manage hyperbolic arguments. The proposed fit strategy and initial grid will allow for hyperbolic orbits")
+parser.add_argument("--E0-max", default=1.060,type=float,help="Maximum range of 'E0' allowed.")
+parser.add_argument("--E0-min", default=1.0,type=float,help="Minimum range of 'E0' allowed.")
+parser.add_argument("--pphi0-max", default=5.4,type=float,help="Maximum range of 'p_phi0' allowed.")
+parser.add_argument("--pphi0-min", default=3.8,type=float,help="Minumum range of 'p_phi0' allowed.")
 parser.add_argument("--assume-nospin",action='store_true',help="If present, the code will not add options to manage precessing spins (the default is aligned spin)")
 parser.add_argument("--assume-precessing-spin",action='store_true',help="If present, the code will add options to manage precessing spins (the default is aligned spin)")
 parser.add_argument("--assume-volumetric-spin",action='store_true',help="If present, the code will assume a volumetric spin prior in its last iterations. If *not* present, the code will adopt a uniform magnitude spin prior in its last iterations. If not present, generally more iterations are taken.")
@@ -450,6 +455,9 @@ def get_observing_run(t):
     event_dict["s2z"] = P.s2z
     event_dict["P"] = P
     event_dict["epoch"]  = 0 # no estimate for now
+    if opts.assume_hyperbolic:
+        event_dict["E0"] = P.E0
+        event_dict["p_phi0"] = P.p_phi0
 elif opts.use_coinc: # If using a coinc through injections and not a GraceDB event.
     # Same code as used before for gracedb
     coinc_file = opts.use_coinc
@@ -1159,6 +1167,8 @@ def crit_m2(delta):
     helper_ile_args += " --data-start-time " + str(data_start_time) + " --data-end-time " + str(data_end_time)  + " --inv-spec-trunc-time 0  --window-shape " + str(window_shape)
 if opts.assume_eccentric:
     helper_ile_args += " --save-eccentricity "
+if opts.assume_hyperbolic:
+    helper_ile_args += " --save-hyperbolic "
 if opts.propose_initial_grid_fisher: # and (P.extract_param('mc')/lal.MSUN_SI < 10.):
     cmd  = "util_AnalyticFisherGrid.py  --inj-file-out  proposed-grid  "
     # Add standard downselects : do not have m1, m2 be less than 1
@@ -1291,6 +1301,12 @@ def lambda_m_estimate(m):
 
     if not (opts.force_initial_grid_size is None):
         grid_size = opts.force_initial_grid_size
+        
+    # hyperbolic grid:
+    if opts.assume_hyperbolic:
+        # note - currently uniform for testing purposes
+        cmd += f" --parameter E0 --parameter-range [{opts.E0_min},{opts.E0_max}] --parameter p_phi0 --parameter-range [{opts.pphi0_min},{opts.pphi0_max}] "
+        
     cmd += " --grid-cartesian-npts  " + str(int(grid_size))
     print(" Executing grid command ", cmd)
     os.system(cmd)
@@ -1387,6 +1403,8 @@ def lambda_m_estimate(m):
 helper_puff_args = " --parameter mc --parameter eta --fmin {} --fref {} ".format(opts.fmin_template,opts.fmin_template)
 if opts.assume_eccentric:
     helper_puff_args += " --parameter eccentricity "
+if opts.assume_hyperbolic:
+    helper_puff_args += " --parameter E0 --parameter p_phi0 "
 
 if event_dict["MChirp"] >25:
     # at high mass, mc/eta correlation weak, don't want to have eta coordinate degeneracy at q=1 to reduce puff proposals  near there
@@ -1707,6 +1725,10 @@ def lambda_m_estimate(m):
 if opts.assume_eccentric:
     with open("helper_convert_args.txt",'w+') as f:
         f.write(" --export-eccentricity ")
+        
+if opts.assume_hyperbolic:
+    with open("helper_convert_args.txt",'w+') as f:
+        f.write(" --export-hyperbolic ")
 
 if opts.propose_fit_strategy:
     with open("helper_puff_max_it.txt",'w') as f:
diff --git a/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode b/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode
index e62a7cce2..219d024c5 100755
--- a/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode
+++ b/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode
@@ -240,6 +240,7 @@ optp.add_option("--internal-waveform-fd-no-condition",action='store_true',help='
 optp.add_option("--internal-waveform-extra-lalsuite-args",type=str,default=None)
 optp.add_option("--verbose",action='store_true')
 optp.add_option("--save-eccentricity", action="store_true")
+optp.add_option("--save-hyperbolic", action="store_true")
 #
 # Add the integration options
 #
@@ -1633,6 +1634,9 @@ def analyze_event(P_list, indx_event, data_dict, psd_dict, fmax, opts,inv_spec_t
         if opts.save_eccentricity:
             # output format when eccentricity is being used
             numpy.savetxt(fname_output_txt, numpy.array([[event_id, m1, m2, P.s1x, P.s1y, P.s1z, P.s2x, P.s2y, P.s2z,  P.eccentricity, log_res+manual_avoid_overflow_logarithm, sqrt_var_over_res,sampler.ntotal, neff ]]))  #dict_return["convergence_test_results"]["normal_integral]"
+        if opts.save_hyperbolic:
+            # output format when hyperbolic is being used
+            numpy.savetxt(fname_output_txt, numpy.array([[event_id, m1, m2, P.s1x, P.s1y, P.s1z, P.s2x, P.s2y, P.s2z,  P.E0, P.p_phi0, log_res+manual_avoid_overflow_logarithm, sqrt_var_over_res,sampler.ntotal, neff ]]))
         elif not (P.lambda1>0 or P.lambda2>0):
           # output format when lambda is NOT used
           if not opts.pin_distance_to_sim:
@@ -1710,6 +1714,8 @@ def analyze_event(P_list, indx_event, data_dict, psd_dict, fmax, opts,inv_spec_t
             samples["alpha4"] =numpy.ones(samples["psi"].shape)*P.eccentricity
             samples["alpha5"] =numpy.ones(samples["psi"].shape)*P.lambda1
             samples["alpha6"] =numpy.ones(samples["psi"].shape)*P.lambda2
+            samples["psi3"] =numpy.ones(samples["psi"].shape)*P.E0
+            samples["beta"] =numpy.ones(samples["psi"].shape)*P.p_phi0
             # Below exist solely to placate XML export; new issue as of latest lalsuite say 7.15+ or so
             samples["alpha2"] = numpy.zeros(samples["psi"].shape)
             samples["alpha3"] = numpy.zeros(samples["psi"].shape)
@@ -1720,7 +1726,7 @@ def analyze_event(P_list, indx_event, data_dict, psd_dict, fmax, opts,inv_spec_t
 #        print(samples['t_ref'] - fiducial_epoch, len(samples['t_ref']))
       xmlutils.append_samples_to_xmldoc(xmldoc, samples)
         # Extra metadata
-      dict_out={"mass1": m1, "mass2": m2, "spin1z": P.s1z, "spin2z": P.s2z, "alpha4": P.eccentricity, "alpha5":P.lambda1, "alpha6":P.lambda2, "event_duration": sqrt_var_over_res, "ttotal": sampler.ntotal}
+      dict_out={"mass1": m1, "mass2": m2, "spin1z": P.s1z, "spin2z": P.s2z, "alpha4": P.eccentricity, "alpha5":P.lambda1, "alpha6":P.lambda2, "psi3":P.E0, "beta":P.p_phi0, "event_duration": sqrt_var_over_res, "ttotal": sampler.ntotal}
 #      if 'distance' in pinned_params:
 #            dict_out['distance'] = pinned_params["distance"]
       converged_result = False
diff --git a/MonteCarloMarginalizeCode/Code/bin/plot_posterior_corner.py b/MonteCarloMarginalizeCode/Code/bin/plot_posterior_corner.py
index a1fd60d31..91ed35b26 100755
--- a/MonteCarloMarginalizeCode/Code/bin/plot_posterior_corner.py
+++ b/MonteCarloMarginalizeCode/Code/bin/plot_posterior_corner.py
@@ -232,6 +232,7 @@ def render_coordinates(coord_names,logparams=[]):
 parser.add_argument("--lnL-cut",default=None,type=float)
 parser.add_argument("--sigma-cut",default=0.4,type=float)
 parser.add_argument("--eccentricity", action="store_true", help="Read sample files in format including eccentricity")
+parser.add_argument("--hyperbolic", action="store_true", help="Read sample files in format including hyperbolic params")
 parser.add_argument("--matplotlib-block-defaults",action="store_true",help="Relies entirely on user to set plot options for plot styles from matplotlibrc")
 parser.add_argument("--no-mod-psi",action="store_true",help="Default is to take psi mod pi. If present, does not do this")
 parser.add_argument("--verbose",action='store_true',help='print matplotlibrc data')
@@ -467,6 +468,9 @@ def render_coordinates(coord_names,logparams=[]):
 if opts.eccentricity:
     print(" Reading composite file, assuming eccentricity-based format ")
     field_names=("indx","m1", "m2",  "a1x", "a1y", "a1z", "a2x", "a2y", "a2z","eccentricity", "lnL", "sigmaOverL", "ntot", "neff")
+if opts.hyperbolic:
+    print(" Reading composite file, assuming hyperbolic-based format ")
+    field_names=("indx","m1", "m2",  "a1x", "a1y", "a1z", "a2x", "a2y", "a2z","E0", "p_phi0", "lnL", "sigmaOverL", "ntot", "neff")
 field_formats = [np.float32 for x in field_names]
 composite_dtype = [ (x,float) for x in field_names] #np.dtype(names=field_names ,formats=field_formats)
 # Load posterior files
diff --git a/MonteCarloMarginalizeCode/Code/bin/util_CleanILE.py b/MonteCarloMarginalizeCode/Code/bin/util_CleanILE.py
index cdccb43bb..8655bd947 100755
--- a/MonteCarloMarginalizeCode/Code/bin/util_CleanILE.py
+++ b/MonteCarloMarginalizeCode/Code/bin/util_CleanILE.py
@@ -29,6 +29,7 @@
 parser = argparse.ArgumentParser(usage="util_CleanILE.py fname1.dat fname2.dat ... ")
 parser.add_argument("fname",action='append',nargs='+')
 parser.add_argument("--eccentricity", action="store_true")
+parser.add_argument("--hyperbolic", action="store_true")
 #Askold: adding specification for tabular eos file
 parser.add_argument("--tabular-eos-file", action="store_true") 
 opts = parser.parse_args()
@@ -54,13 +55,16 @@
         if opts.eccentricity:
             indx, m1,m2, s1x,s1y,s1z,s2x,s2y,s2z,ecc, lnL, sigmaOverL, ntot, neff = line
             col_intrinsic = 10
+        if opts.hyperbolic:
+            indx, m1, m2, s1x, s1y, s1z, s2x, s2y, s2z, E0, p_phi0, lnL, sigmaOverL, ntot, neff = line
+            col_intrinsic = 11
         elif len(line) == 13 and (not tides_on) and (not distance_on):  # strip lines with the wrong length
             indx, m1,m2, s1x,s1y,s1z,s2x,s2y,s2z,lnL, sigmaOverL, ntot, neff = line
         elif  len(line) == 14:
             distance_on=True
             col_intrinsic=10
             indx, m1,m2, s1x,s1y,s1z,s2x,s2y,s2z,dist, lnL, sigmaOverL, ntot, neff = line
-        elif len(line)==15:
+        elif len(line)==15 and not opts.hyperbolic:
             tides_on  = True
             col_intrinsic =11
             indx, m1,m2, s1x,s1y,s1z,s2x,s2y,s2z, lambda1,lambda2,lnL, sigmaOverL, ntot, neff = line
@@ -94,6 +98,8 @@
 
     if opts.eccentricity:
         print(-1,  key[0],key[1], key[2], key[3],key[4], key[5],key[6], key[7], key[8], lnLmeanMinusLmax+lnLmax, sigmaNetOverL, np.sum(ntot), -1)
+    elif opts.hyperbolic:
+        print(-1,  key[0],key[1], key[2], key[3],key[4], key[5],key[6], key[7], key[8],key[9], lnLmeanMinusLmax+lnLmax, sigmaNetOverL, np.sum(ntot), -1)
     elif tides_on:
         print(-1,  key[0],key[1], key[2], key[3],key[4], key[5],key[6], key[7], key[8],key[9], lnLmeanMinusLmax+lnLmax, sigmaNetOverL, np.sum(ntot), -1)
     elif distance_on:
diff --git a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py
index bd9d3555b..eb96011ef 100755
--- a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py
+++ b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py
@@ -338,6 +338,11 @@ def extract_combination_from_LI(samples_LI, p):
 parser.add_argument("--internal-gmm-memory-chisquared-factor",default=None,type=float,help="Multiple of the number of degrees of freedom to save. 5 is a part in 10^6, 4 is 10^{-4}, and None keeps all up to lnL_offset.  Note that low-weight points can contribute notably to n_eff, and it can be dangerous to assume a simple chisquared likelihood!  Provided in case we need very long runs")
 parser.add_argument("--assume-eos-but-primary-bh",action='store_true',help="Special case of known EOS, but primary is a BH")
 parser.add_argument("--use-eccentricity", action="store_true")
+parser.add_argument("--use-hyperbolic", action="store_true")
+parser.add_argument("--E0-max", default=1.060,type=float,help="Maximum range of 'E0' allowed.")
+parser.add_argument("--E0-min", default=1.0,type=float,help="Minumum range of 'E0' allowed.")
+parser.add_argument("--pphi0-max", default=5.4,type=float,help="Maximum range of 'p_phi0' allowed.")
+parser.add_argument("--pphi0-min", default=3.8,type=float,help="Minimum range of 'p_phi0' allowed.")
 parser.add_argument("--tripwire-fraction",default=0.05,type=float,help="Fraction of nmax of iterations after which n_eff needs to be greater than 1+epsilon for a small number epsilon")
 
 # FIXME hacky options added by me (Liz) to try to get my capstone project to work.
@@ -357,6 +362,12 @@ def extract_combination_from_LI(samples_LI, p):
     opts.no_adapt_parameter =[] # needs to default to empty list
 ECC_MAX = opts.ecc_max
 ECC_MIN = opts.ecc_min
+
+E0_MAX = opts.E0_max
+E0_MIN = opts.E0_min
+PPHI0_MAX = opts.pphi0_max
+PPHI0_MIN = opts.pphi0_min
+
 no_plots = no_plots |  opts.no_plots
 lnL_shift = 0
 lnL_default_large_negative = -500
@@ -817,6 +828,12 @@ def eccentricity_prior(x):
 def precession_prior(x):
     return 0.5*np.ones(x.shape) # uniform over the interval [0.0, 2.0]
 
+def initial_energy_prior(x):
+    return np.ones(x.shape) / (E0_MAX-E0_MIN) # uniform over the interval [E0_MIN, E0_MAX]
+
+def initial_angmom_prior(x):
+    return np.ones(x.shape) / (PPHI0_MAX-PPHI0_MIN) # uniform over the interval [PPHI0_MIN, PPHI0_MAX]
+
 def unnormalized_uniform_prior(x):
     return np.ones(x.shape)
 def unnormalized_log_prior(x):
@@ -864,6 +881,8 @@ def normalized_zbar_prior(z):
     # Other priors
     'eccentricity':eccentricity_prior,
     'chi_pavg':precession_prior,
+    'E0':initial_energy_prior,
+    'p_phi0':initial_angmom_prior,
     'mu1': unnormalized_log_prior,
     'mu2': unnormalized_uniform_prior
 }
@@ -881,7 +900,9 @@ def normalized_zbar_prior(z):
   'lambda_plus':[0.01,lambda_plus_max],
   'lambda_minus':[-lambda_max,lambda_max],  # will include the true region always...lots of overcoverage for small lambda, but adaptation will save us.
   'eccentricity':[ECC_MIN, ECC_MAX],
-  'chi_pavg':[0.0,2.0],  
+  'chi_pavg':[0.0,2.0],
+  'E0':[E0_MIN,E0_MAX],
+  'p_phi0':[PPHI0_MIN,PPHI0_MAX],
   # strongly recommend you do NOT use these as parameters!  Only to insure backward compatibility with LI results
   'LambdaTilde':[0.01,5000],
   'DeltaLambdaTilde':[-500,500],
@@ -1576,6 +1597,11 @@ def fit_gp_sparse(x):
 elif opts.use_eccentricity:
     print(" Eccentricity input: [",ECC_MIN, ", ",ECC_MAX, "]")
     col_lnL += 1
+elif opts.use_hyperbolic:
+    # add two columns for hyperbolic params
+    print("E0 import: [",E0_MIN, ", ",E0_MAX, "]")
+    print("p_phi0 import: [",PPHI0_MIN, ", ",PPHI0_MAX, "]")
+    col_lnL += 2
 if opts.input_distance:
     print(" Distance input")
     col_lnL +=1
@@ -1669,6 +1695,9 @@ def fit_gp_sparse(x):
         P.eos_table_index = line[11]
     if opts.use_eccentricity:
         P.eccentricity = line[9]
+    if opts.use_hyperbolic:
+        P.E0 = line[9]
+        P.p_phi0 = line[10]
     if opts.input_distance:
         P.dist = lal.PC_SI*1e6*line[9]  # Incompatible with tides, note!
     
diff --git a/MonteCarloMarginalizeCode/Code/bin/util_ILEdagPostprocess.sh b/MonteCarloMarginalizeCode/Code/bin/util_ILEdagPostprocess.sh
index fabaa8351..90f76cd45 100755
--- a/MonteCarloMarginalizeCode/Code/bin/util_ILEdagPostprocess.sh
+++ b/MonteCarloMarginalizeCode/Code/bin/util_ILEdagPostprocess.sh
@@ -27,6 +27,12 @@ then
 else
     util_CleanILE.py ${RND}_tmp.dat $3 | sort -rg -k10 > $BASE_OUT.composite
 fi
+if [ "$3" == '--hyperbolic' ]
+then
+    util_CleanILE.py ${RND}_tmp.dat $3 | sort -rg -k12 > $BASE_OUT.composite
+else
+    util_CleanILE.py ${RND}_tmp.dat $3 | sort -rg -k10 > $BASE_OUT.composite
+fi
 
 # Manifest
 rm -f ${BASE_OUT}.manifest
diff --git a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py
index 59d18d36e..6416bec9e 100755
--- a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py
+++ b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py
@@ -157,6 +157,11 @@ def unsafe_parse_arg_string_dict(my_argstr):
 parser.add_argument("--assume-matter-but-primary-bh",action='store_true',help="If present, the code will add options necessary to manage tidal arguments for the smaller body ONLY. (Usually pointless)")
 parser.add_argument("--internal-tabular-eos-file",type=str,default=None,help="Tabular file of EOS to use.  The default prior will be UNIFORM in this table!")
 parser.add_argument("--assume-eccentric",action='store_true', help="Add eccentric options for each part of analysis")
+parser.add_argument("--assume-hyperbolic",action='store_true', help="Add hyperbolic options for each part of analysis")
+parser.add_argument("--force-E0-max",default=None,type=float,help="Provide this value to override the value of E0-max provided")
+parser.add_argument("--force-E0-min",default=None,type=float,help="Provide this value to override the value of E0-min provided")
+parser.add_argument("--force-pphi0-max",default=None,type=float,help="Provide this value to override the value of p_phi0-max provided")
+parser.add_argument("--force-pphi0-min",default=None,type=float,help="Provide this value to override the value of p_phi0-min provided")
 parser.add_argument("--assume-lowlatency-tradeoffs",action='store_true', help="Force analysis with various low-latency tradeoffs (e.g., drop spin 2, use aligned, etc)")
 parser.add_argument("--assume-highq",action='store_true', help="Force analysis with the high-q strategy, neglecting spin2. Passed to 'helper'")
 parser.add_argument("--assume-well-placed",action='store_true',help="If present, the code will adopt a strategy that assumes the initial grid is very well placed, and will minimize the number of early iterations performed. Not as extrme as --propose-flat-strategy")
@@ -469,6 +474,7 @@ def unsafe_parse_arg_string_dict(my_argstr):
 
 is_analysis_precessing =False
 is_analysis_eccentric =False
+is_analysis_hyperbolic =False
 if opts.approx == "SEOBNRv3" or opts.approx == "NRSur7dq2" or opts.approx == "NRSur7dq4" or (opts.approx == 'SEOBNv3_opt') or (opts.approx == 'IMRPhenomPv2') or (opts.approx =="SEOBNRv4P" ) or (opts.approx == "SEOBNRv4PHM") or (opts.approx == "SEOBNRv5PHM") or ('SpinTaylor' in opts.approx) or ('IMRPhenomTP' in opts.approx or ('IMRPhenomXP' in opts.approx)):
         is_analysis_precessing=True
 if opts.assume_precessing:
@@ -477,6 +483,8 @@ def unsafe_parse_arg_string_dict(my_argstr):
         is_analysis_precessing = False
 if opts.assume_eccentric:
         is_analysis_eccentric = True
+if opts.assume_hyperbolic:
+        is_analysis_hyperbolic = True
 
 
 dirname_run = gwid+ "_" + opts.calibration+ "_"+ opts.approx+"_fmin" + str(fmin) +"_fmin-template"+str(fmin_template) +"_lmax"+str(opts.l_max) + "_"+opts.spin_magnitude_prior
@@ -494,6 +502,8 @@ def unsafe_parse_arg_string_dict(my_argstr):
     dirname_run += "_with_matter"
 if opts.assume_eccentric:
     dirname_run += "_with_eccentricity"
+if opts.assume_hyperbolic:
+    dirname_run += "_with_hyperbolic"
 if opts.no_matter:
     dirname_run += "_no_matter"
 if opts.assume_highq:
@@ -617,6 +627,20 @@ def unsafe_parse_arg_string_dict(my_argstr):
         npts_it = 1500
 if is_analysis_eccentric:
     cmd += " --assume-eccentric "
+if is_analysis_hyperbolic:
+    cmd += " --assume-hyperbolic "
+    if not(opts.force_E0_max is None):
+        E0_max = opts.force_E0_max
+        cmd += " --E0-max {}  ".format(E0_max)
+    if not(opts.force_E0_min is None):
+        E0_min = opts.force_E0_min
+        cmd += " --E0-min {}  ".format(E0_min)
+    if not(opts.force_pphi0_max is None):
+        pphi0_max = opts.force_pphi0_max
+        cmd += " --pphi0-max {}  ".format(pphi0_max)
+    if not(opts.force_pphi0_min is None):
+        pphi0_min = opts.force_pphi0_min
+        cmd += " --pphi0-min {}  ".format(pphi0_min)
 if opts.assume_highq:
     cmd+= ' --assume-highq  --force-grid-stretch-mc-factor 2'  # the mc range, tuned to equal-mass binaries, is probably too narrow. Workaround until fixed in helper
     npts_it =1000
@@ -1040,6 +1064,23 @@ def unsafe_parse_arg_string_dict(my_argstr):
         if not(opts.force_ecc_min is None):
             ecc_min = opts.force_ecc_min
             line += " --ecc-min {}  ".format(ecc_min)
+    if opts.assume_hyperbolic:
+        if not(opts.internal_use_aligned_phase_coordinates):
+            line = line.replace('parameter mc', 'parameter mc --parameter E0 --parameter p_phi0 --use-hyperbolic')
+        else:
+            line = line.replace('parameter-nofit mc', 'parameter-nofit mc --parameter E0 --parameter p_phi0 --use-hyperbolic')
+        if not(opts.force_E0_max is None):
+            E0_max = opts.force_E0_max
+            line += " --E0-max {}  ".format(E0_max)
+        if not(opts.force_E0_min is None):
+            E0_max = opts.force_E0_min
+            line += " --E0-min {}  ".format(E0_min)
+        if not(opts.force_pphi0_max is None):
+            pphi0_max = opts.force_pphi0_max
+            line += " --pphi0-max {}  ".format(pphi0_max)
+        if not(opts.force_pphi0_min is None):
+            pphi0_max = opts.force_pphi0_min
+            line += " --pphi0-min {}  ".format(pphi0_min)
     if not(opts.manual_extra_cip_args is None):
         line += " {} ".format(opts.manual_extra_cip_args)  # embed with space on each side, avoid collisions
     line += "\n"
@@ -1102,6 +1143,16 @@ def unsafe_parse_arg_string_dict(my_argstr):
     puff_max_it +=5   # make sure we resolve the correlations
 if opts.assume_eccentric:
         puff_params += " --parameter eccentricity --downselect-parameter eccentricity --downselect-parameter-range '[0,0.9]' "
+if opts.assume_hyperbolic:
+        puff_params += " --parameter E0 --downselect-parameter E0 --downselect-parameter-range '[1.0,1.060]' --parameter E0 --downselect-parameter p_phi0 --downselect-parameter-range '[3.8,5.4]' "
+        if not(opts.force_E0_max is None and opts.force_E0_min is None):
+            E0_max = opts.force_E0_max
+            E0_min = opts.force_E0_min
+            puff_params.replace("--parameter E0 --downselect-parameter E0 --downselect-parameter-range '[1.0,1.060]'", f"--parameter E0 --downselect-parameter E0 --downselect-parameter-range '[{E0_min},{E0_max}]'")
+        if not(opts.force_pphi0_max is None and opts.force_pphi0_min is None):
+            pphi0_max = opts.force_pphi0_max
+            pphi0_min = opts.force_pphi0_min
+            puff_params.replace("--parameter pphi0 --downselect-parameter pphi0 --downselect-parameter-range '[3.8,5.4]'", f"--parameter pphi0 --downselect-parameter pphi0 --downselect-parameter-range '[{pphi0_min},{pphi0_max}]'")
 if opts.assume_highq:
     puff_params = puff_params.replace(' delta_mc ', ' eta ')  # use natural coordinates in the high q strategy. May want to do this always
     puff_max_it +=3
@@ -1190,7 +1241,7 @@ def unsafe_parse_arg_string_dict(my_argstr):
 if opts.use_subdags:
     cepp = "create_event_parameter_pipeline_AlternateIteration"
 cmd =cepp+ "  --ile-n-events-to-analyze {} --input-grid proposed-grid.xml.gz --ile-exe  `which integrate_likelihood_extrinsic_batchmode`   --ile-args `pwd`/args_ile.txt --cip-args-list args_cip_list.txt --test-args args_test.txt --request-memory-CIP {} --request-memory-ILE 4096 --n-samples-per-job ".format(n_jobs_per_worker,cip_mem) + str(npts_it) + " --working-directory `pwd` --n-iterations " + str(n_iterations) + " --n-iterations-subdag-max {} ".format(opts.internal_n_iterations_subdag_max) + "  --n-copies {} ".format(opts.ile_copies) + "   --ile-retries "+ str(opts.ile_retries) + " --general-retries " + str(opts.general_retries)
-if opts.assume_matter or opts.assume_eccentric:
+if opts.assume_matter or opts.assume_eccentric or opts.assume_hyperbolic:
     cmd +=  " --convert-args `pwd`/helper_convert_args.txt "
 if not(opts.ile_runtime_max_minutes is None):
     cmd += " --ile-runtime-max-minutes {} ".format(opts.ile_runtime_max_minutes)
@@ -1198,6 +1249,8 @@ def unsafe_parse_arg_string_dict(my_argstr):
     cmd+= " --puff-exe `which util_ParameterPuffball.py` --puff-cadence 1 --puff-max-it " + str(puff_max_it)+ " --puff-args `pwd`/args_puff.txt "
 if opts.assume_eccentric:
     cmd += " --use-eccentricity "
+if opts.assume_hyperbolic:
+    cmd += " --use-hyperbolic "
 if opts.calibration_reweighting and (not opts.bilby_pickle_file):
     cmd += " --calibration-reweighting --calibration-reweighting-exe `which calibration_reweighting.py` --bilby-ini-file {} --bilby-pickle-exe `which bilby_pipe_generation` ".format(str(opts.bilby_ini_file))
     if opts.calibration_reweighting_count:
diff --git a/MonteCarloMarginalizeCode/Code/bin/util_SimInspiralToCoinc.py b/MonteCarloMarginalizeCode/Code/bin/util_SimInspiralToCoinc.py
index 072a1feb7..8617d046e 100755
--- a/MonteCarloMarginalizeCode/Code/bin/util_SimInspiralToCoinc.py
+++ b/MonteCarloMarginalizeCode/Code/bin/util_SimInspiralToCoinc.py
@@ -6,7 +6,7 @@
 import json
 
 
-def ilwd_base():
+def ilwd_base(arg=None):
     return 0
 try:
     from  glue.ligolw import ilwd

From cc163fae841929d6d6628a770898ed6b450e26ce Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Mon, 6 May 2024 12:36:34 -0700
Subject: [PATCH 02/50] single_inject with hyperbolic options

---
 .../test/single_inject/single_inject_RIFT.py  | 503 ++++++++++++++++++
 1 file changed, 503 insertions(+)
 create mode 100644 MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py

diff --git a/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py b/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py
new file mode 100644
index 000000000..97127c45f
--- /dev/null
+++ b/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py
@@ -0,0 +1,503 @@
+#! /usr/bin/env python
+#
+# GOAL: single-button injection setup from ini file
+#
+#
+
+################################################################
+#### Code Description ####
+################################################################
+
+# This code serves to create a full RIFT injection workflow, pulling all relevant information (including waveform parameters) from a config.ini file.
+# Please see the example ini file (BBH_inject_example.ini)
+
+
+## Basic Code Process: ##
+
+# The code process is as follows, starting from reading in the ini file: (example: single_inject_RIFT.py -- /home/albert.einstein/path/to/config/BBH_inject_example.ini)
+# Create a MDC (mock data challenge) SimInspiral table xml from the waveform parameters
+# uses LALWriteFrame.py to create frame files for the specified detectors, stored in the signal_frames dir (also combined_frames if there is no noise)
+# if specified with --use-noise, adds noise frames and dumps the combined into the combined_frames dir
+# calculates SNR with util_FrameZeroNoiseSNR.py, uses network SNR as the initial guess
+# writes the coinc.xml table
+# runs util_RIFT_pseudo_pipe.py, which sets up the RIFT analysis workflow. In the resultant analysis directory the user needs to only submit the dag via condor_submit_dag
+
+
+## Options: ##
+
+# --use-osg will allow you to run on the IGWN pool of distributed computing resources, for which you need to have the environment variables:
+# SINGULARITY_RIFT_IMAGE: path to the RIFT Singularity image (container) that will run RIFT
+# SINGULARITY_BASE_EXE_DIR=/usr/local/bin/
+
+# --add-extrinsic will enable calculation of the extrinsic posteriors in the final iteration
+
+# --force-cpu will disable the GPU code pathway, this is useful for running smaller or custom jobs
+
+# --use-noise will add noise to the created signal frames
+
+# --bypass frames will skip to the util_RIFT_pseudo_pipe.py step, assuming that all the file structure is otherwise in place.
+# this is useful for recreating an analysis off of existing frames without waiting for them to generate.
+
+# --just-frames will exit the program after the SNR report is completed.
+# when combined with --bypass frames, allows for two-stage injection creation
+
+# --use-hyperbolic passes further options to util_RIFT_pseudo_pipe.py, enabling analysis of hyperbolic systems
+
+
+## Other notes: ##
+
+# if you're using a distance marginalization table, it will take a while to create the injection. You'll see:
+#  IntegrationWarning: The integral is probably divergent, or slowly convergent.
+# but it is working, give it a few minutes
+
+
+
+################################################################
+#### Preamble ####
+################################################################
+
+import numpy as np
+import argparse
+import os
+import sys
+sys.path.append(os.getcwd())
+import shutil
+import ast
+import glob
+import re
+import configparser
+
+from RIFT.misc.dag_utils import mkdir
+import RIFT.lalsimutils as lalsimutils
+import lal
+import lalsimulation as lalsim
+
+from gwpy.timeseries import TimeSeries
+from matplotlib import pyplot as plt
+
+# Backward compatibility
+from RIFT.misc.dag_utils import which
+lalapps_path2cache = which('lal_path2cache')
+if lalapps_path2cache == None:
+    lalapps_path2cache =  which('lalapps_path2cache')
+    
+################################################################
+#### Functions & Definitions ####
+################################################################
+
+
+def unsafe_config_get(config,args,verbose=False):
+    if verbose:
+        print(" Retrieving ", args, end=' ') 
+        print(" Found ",eval(config.get(*args)))
+    return eval( config.get(*args))
+
+def add_frames(channel,input_frame, noise_frame,combined_frame):
+    # note - this needs to be updated to work with the current install regardless of user
+    exe = os.path.abspath(os.path.join(os.path.dirname(os.path.abspath(__file__)), '../pp/add_frames.py'))
+    if not(os.path.dirname(exe) in os.environ['PATH'].split(os.pathsep)):
+        print('add_frames.py not found, adding to PATH...')
+        os.environ['PATH'] += os.pathsep + os.path.dirname(exe)
+    cmd = exe + " " + channel + " " + input_frame + " " + noise_frame + " " + combined_frame 
+    print(cmd)
+    os.system(cmd)
+
+lsu_MSUN=lal.MSUN_SI
+lsu_PC=lal.PC_SI        
+
+################################################################
+#### Options Parse ####
+################################################################
+        
+parser = argparse.ArgumentParser()
+parser.add_argument("--use-ini", default=None, type=str, required=True, help="REQUIRED. Pass ini file for parsing. Intended to reproduce lalinference_pipe functionality. Overrides most other arguments.")
+parser.add_argument("--use-osg",action='store_true',help="Attempt OSG operation. Command-line level option, not at ini level to make portable results")
+parser.add_argument("--add-extrinsic",action='store_true',help="Add extrinsic posterior.  Corresponds to --add-extrinsic --add-extrinsic-time-resampling --batch-extrinsic for pipeline")
+parser.add_argument("--force-cpu",action='store_true',help="Forces avoidance of GPUs")
+parser.add_argument("--use-noise",action='store_true',help="Combine clean signal injection with fiducial noise frames. Note that you must have pycbc installed.")
+parser.add_argument("--bypass-frames",action='store_true',help="Skip making mdc and frame files, use with caution")
+parser.add_argument("--just-frames",action='store_true',help="Stop after making frame files, use with caution")
+parser.add_argument("--use-hyperbolic",action='store_true',help="Adds hyperbolic options, requires TEOBResumS approx. NOTE: development still in progress.")
+opts =  parser.parse_args()
+
+config = configparser.ConfigParser(allow_no_value=True) #SafeConfigParser deprecated from py3.2
+
+ini_path = os.path.abspath(opts.use_ini)
+
+# read in config
+config.read(ini_path)
+
+bypass_frames=opts.bypass_frames 
+just_frames=opts.just_frames
+use_noise=opts.use_noise 
+
+dist_marg_arg = config.get('rift-pseudo-pipe','internal-marginalize-distance') 
+dist_marg = True if dist_marg_arg.lower() == 'true' else False
+
+################################################################
+#### Initialize ####
+################################################################
+
+# Create, go to working directory
+working_dir = config.get('init','working_directory')
+print(" Working dir ", working_dir)
+mkdir(working_dir)
+os.chdir(working_dir)
+working_dir_full = os.getcwd()
+
+# load IFO list from config
+ifos = unsafe_config_get(config,['analysis','ifos'])
+
+################################################################
+#### Load Injection Parameters from config into lalsimutils ####
+################################################################
+
+P_list = []
+P = lalsimutils.ChooseWaveformParams()
+
+## Waveform Model ##
+approx_str = config.get('engine','approx')
+P.approx=lalsim.GetApproximantFromString(approx_str)
+
+## Intrinsic Parameters ##
+P.m1 = float(config.get('injection-parameters','m1'))*lsu_MSUN
+P.m2 = float(config.get('injection-parameters','m2'))*lsu_MSUN
+P.s1x = float(config.get('injection-parameters','s1x'))
+P.s1y = float(config.get('injection-parameters','s1y'))
+P.s1z = float(config.get('injection-parameters','s1z'))
+P.s2x = float(config.get('injection-parameters','s2x'))
+P.s2y = float(config.get('injection-parameters','s2y'))
+P.s2z = float(config.get('injection-parameters','s2z'))
+if opts.use_hyperbolic:
+    P.E0 = float(config.get('injection-parameters','E0'))
+    P.p_phi0 = float(config.get('injection-parameters','p_phi0'))
+## Extrinsic Parameters ##
+P.dist = float(config.get('injection-parameters','dist'))*1e6*lsu_PC
+fiducial_event_time = float(config.get('injection-parameters','event_time')) # for frame addition
+P.tref = fiducial_event_time
+P.theta = float(config.get('injection-parameters','dec'))
+P.phi = float(config.get('injection-parameters','ra'))
+P.incl = float(config.get('injection-parameters','incl'))
+P.psi = float(config.get('injection-parameters','psi'))
+P.phiref = float(config.get('injection-parameters','phase'))
+# Other waveform settings
+P.fmin = float(config.get('injection-parameters','fmin'))
+
+P_list.append(P)
+
+# create mdc.xml.gz - contains all injection params
+if bypass_frames:
+    print("Skipping mdc creation")
+    pass
+else:
+    lalsimutils.ChooseWaveformParams_array_to_xml(P_list,"mdc")
+
+    
+    
+################################################################    
+#### Create frame files and cache file ####
+################################################################
+
+if bypass_frames:
+    print("Skipping frame creation")
+    indx = 0
+    pass
+else:
+    # create clean signal frames #
+    print(" === Writing signal files === ")
+
+    mkdir('signal_frames')
+    # must be compatible with duration of the noise frames
+    t_start = int(fiducial_event_time)-150
+    t_stop = int(fiducial_event_time)+150
+
+    indx = 0 # need an event number, meaningless otherwise
+    target_subdir = 'signal_frames/event_{}'.format(indx) # this is where the signal frames will go
+    mkdir(working_dir_full+"/"+target_subdir)
+    os.chdir(working_dir_full+"/"+target_subdir)
+    # here we're in the signal_frames/event_0 directory.
+    # Loop over instruments, write frame files and cache
+    for ifo in ifos:
+        cmd = "util_LALWriteFrame.py --inj " + working_dir_full+"/mdc.xml.gz --event {} --start {}  --stop {}  --instrument {} --approx {}".format(indx, t_start,t_stop,ifo, approx_str)
+        print(cmd)
+        os.system(cmd)
+        
+    # Here we plot the signal frames    
+    gwf_files = [os.path.abspath(file) for file in os.listdir(os.getcwd()) if file.endswith('.gwf')]
+    
+    for index, path in enumerate(gwf_files):
+        file = os.path.basename(path)
+        # assign correct channel name
+        if file.startswith('H'):
+            channel = 'H1:FAKE-STRAIN'
+        elif file.startswith('L'):
+            channel = 'L1:FAKE-STRAIN'
+        elif file.startswith('V'):
+            channel = 'V1:FAKE-STRAIN'
+            
+        data = TimeSeries.read(path, channel)
+        strain = np.array(data.value)
+        sample_times = np.array(data.times)
+        
+        # plot domain - don't do this for combined frames
+        max_amp = np.argmax(strain)        
+        plot_domain = [sample_times[max_amp] - 1.0, sample_times[max_amp] + 0.5]
+        
+        plt.plot(sample_times, data, label=f'{channel}')
+        plt.xlim(plot_domain)
+        plt.legend(loc='best')
+        plt.savefig(f'{channel}_plot.png', bbox_inches='tight')
+        plt.close()
+        
+    if not use_noise:
+        # copy the contents of signal_frames/event_0 to combined_frames/event_0
+        os.chdir(working_dir_full)
+        mkdir('combined_frames')
+        os.chdir(working_dir_full)
+        target_subdir='combined_frames/event_{}'.format(indx)
+        shutil.copytree(working_dir_full+'/signal_frames/event_{}'.format(indx), target_subdir)
+        
+    else:
+
+        print(" === Using fiducial noise frames === ")
+
+        noise_frame_dir = config.get("make_data","fiducial_noise_frames")
+
+        print(" === Joining synthetic signals to reference noise === ")
+
+        seglen_actual = t_stop - t_start
+        os.chdir(working_dir_full)
+        mkdir('combined_frames')
+        os.chdir(working_dir_full)
+        target_subdir='combined_frames/event_{}'.format(indx)
+        mkdir(target_subdir)
+        os.chdir(working_dir_full+"/"+target_subdir)
+        for ifo in ifos:
+            fname_input = working_dir_full+"/signal_frames/event_{}/{}-fake_strain-{}-{}.gwf".format(indx,ifo[0],int(t_start),seglen_actual)
+            fname_noise= noise_frame_dir + "/" + ifo + ("/%08d.gwf" % indx)
+            fname_output = working_dir_full+"/{}/{}-combined-{}-{}.gwf".format(target_subdir,ifo[0],int(t_start),seglen_actual)
+            print(ifo+":FAKE-STRAIN")
+            print(fname_input)
+            print(fname_noise)
+            print(fname_output)
+            add_frames(ifo+":FAKE-STRAIN",fname_input, fname_noise,fname_output)
+            
+        # Here we plot the combined frames    
+        gwf_files = [os.path.abspath(file) for file in os.listdir(os.getcwd()) if file.endswith('.gwf')]
+
+        for index, path in enumerate(gwf_files):
+            file = os.path.basename(path)
+            # assign correct channel name
+            if file.startswith('H'):
+                channel = 'H1:FAKE-STRAIN'
+            elif file.startswith('L'):
+                channel = 'L1:FAKE-STRAIN'
+            elif file.startswith('V'):
+                channel = 'V1:FAKE-STRAIN'
+            
+            data = TimeSeries.read(path, channel)
+
+            plt.plot(data, label=f'{channel}')
+            plt.legend(loc='best')
+            plt.savefig(f'{channel}_plot.png', bbox_inches='tight')
+            plt.close()
+            
+    # write frame paths to cache
+    os.chdir(working_dir_full)
+    target_subdir='combined_frames/event_{}'.format(indx)
+    os.chdir(working_dir_full+"/"+target_subdir)
+    cmd = "/bin/find .  -name '*.gwf' | {} > signals.cache".format(lalapps_path2cache)
+    os.system(cmd)
+    
+    
+if not opts.use_noise:
+    # Calculate SNRs for zero-noise frames
+    target_subdir='combined_frames/event_{}'.format(indx)
+    target_subdir_full = working_dir_full+"/"+target_subdir
+    os.chdir(target_subdir_full)
+    if os.path.exists(target_subdir_full + '/snr-report.txt'):
+        print('SNR report already exists, skipping')
+        pass
+    else:
+        cmd = "util_FrameZeroNoiseSNR.py --cache signals.cache --psd lalsim.SimNoisePSDaLIGOZeroDetHighPower"
+        os.system(cmd)
+    with open('snr-report.txt', 'r') as file:
+        snr_dict = ast.literal_eval(file.read())
+    file.close()
+    snr_guess = snr_dict['Network'] # passed to pseudo_pipe as --force-hint-snr
+else:
+    # need a way to estimate SNR for the noise frames, skipping for now
+    print('Skipping SNR guess with noise frames')
+    pass
+    
+if opts.just_frames:
+    print('Exiting after frame creation!')
+    sys.exit(0)
+
+################################################################
+#### RIFT Workflow Generation ####
+################################################################ 
+
+
+# check for distance marginalization lookup table - primarily for testing
+if bypass_frames:
+    os.chdir(working_dir_full)
+    if os.path.exists(working_dir_full+'/'+'distance_marginalization_lookup.npz'):
+        print("Using existing distance marginalization table")
+        dist_marg_exists = True
+        dist_marg_path = working_dir_full + "/distance_marginalization_lookup.npz"
+    else:
+        if dist_marg:
+            print("Distance marginalization table not found, it will be made during the workflow")
+        dist_marg_exists = False
+else:
+    dist_marg_exists = False
+    
+
+os.chdir(working_dir_full)
+
+
+# write coinc file
+cmd = "util_SimInspiralToCoinc.py --sim-xml mdc.xml.gz --event {}".format(indx)
+for ifo in ifos:
+    cmd += "  --ifo {} ".format(ifo)
+os.system(cmd)
+
+# Designate rundir 
+rundir = config.get('init','rundir')
+    
+#point to signals.cache
+target_file = 'combined_frames/event_{}/signals.cache'.format(indx)
+
+# RIFT analysis setup
+cmd = 'util_RIFT_pseudo_pipe.py --use-coinc `pwd`/coinc.xml --use-ini {} --use-rundir `pwd`/{} --fake-data-cache `pwd`/{}'.format(ini_path, rundir, target_file)
+if opts.add_extrinsic:
+    cmd += ' --add-extrinsic --add-extrinsic-time-resampling --batch-extrinsic '
+if not opts.use_noise:
+    cmd += ' --force-hint-snr {} '.format(snr_guess)
+if dist_marg_exists:
+    cmd += ' --internal-marginalize-distance-file {} '.format(dist_marg_path)
+if opts.force_cpu:
+    if opts.use_osg:
+        print("Don't force CPUs on the OSG, exiting.")
+        sys.exit(1)
+    else:
+        print('Avoiding GPUs')
+        cmd += ' --ile-no-gpu '
+if opts.use_osg:
+    print('Configuring for OSG operation')
+    cmd += ' --use-osg '
+    cmd += ' --use-osg-file-transfer '
+    cmd += ' --ile-retries 10 '
+    # note: --use-osg-cip not used by default - can set this in the ini file
+if opts.use_hyperbolic:
+    print('Adding hyperbolic options')
+    cmd += ' --assume-hyperbolic '
+
+    
+os.system(cmd)   
+          
+################################################################
+#### Post-Workflow Cleanup ####
+################################################################
+
+# copy over PSDs
+psd_dict = unsafe_config_get(config,["make_psd",'fiducial_psds'])
+for psd in psd_dict:
+    psd_path = psd_dict[psd]
+    shutil.copy2(psd_path, rundir)
+
+# copy config file into rundir/reproducibility
+repro_path = working_dir_full+"/"+rundir+"/"+"reproducibility"
+shutil.copy2(ini_path, repro_path)
+
+if opts.use_osg:
+    # make frames_dir and copy frame files for transfer
+    os.chdir(working_dir_full+"/"+rundir)
+    mkdir('frames_dir')
+    os.chdir(working_dir_full+"/combined_frames/event_{}".format(indx))
+    for filename in os.listdir(os.getcwd()):
+        if filename.endswith('.gwf'):
+            pathname = os.path.join(os.getcwd(), filename)
+            if os.path.isfile(pathname):
+                shutil.copy2(pathname, working_dir_full+"/"+rundir+"/frames_dir")
+                
+    
+# for ease of testing, copy the distance marginalization table to the top-level dir
+if dist_marg:
+    if os.path.exists(working_dir_full+'/'+'distance_marginalization_lookup.npz'):
+        pass
+    else:
+        print('Copying distance marginalization')
+        os.chdir(working_dir_full+"/"+rundir)
+        shutil.copy2('distance_marginalization_lookup.npz', working_dir_full)
+    
+
+    
+# change name of dag
+os.chdir(working_dir_full+"/"+rundir)
+os.rename('marginalize_intrinsic_parameters_BasicIterationWorkflow.dag', working_dir + '.dag')
+
+
+## Update ILE requirements to avoid bad hosts ##
+
+excluded_hosts = []
+
+host_file = '/home/richard.oshaughnessy/igwn_feedback/rift_avoid_hosts.txt'
+
+with open(host_file, 'r') as file:
+    for line in file:
+        hostname = line.strip()
+        excluded_hosts.append(f"(TARGET.Machine =!= '{hostname}')")
+
+avoid_string = "&&".join(excluded_hosts)
+
+# List of file paths to modify
+file_paths = ["ILE.sub","ILE_puff.sub","iteration_4_cip/ILE.sub","iteration_4_cip/ILE_puff.sub"]
+
+if opts.add_extrinsic:
+    file_paths.append("ILE_extr.sub")
+    
+if opts.use_hyperbolic:
+    file_paths = ["ILE.sub","ILE_puff.sub","iteration_2_cip/ILE.sub","iteration_2_cip/ILE_puff.sub"]
+
+for file_path in file_paths:
+    # Construct the full file path
+    full_path = os.path.join(os.getcwd(), file_path)
+
+    # Read the contents of the file
+    with open(full_path, "r") as file:
+        lines = file.readlines()
+        
+    # Find the index of the line containing the 'requirements' command
+    requirements_index = None
+    for i, line in enumerate(lines):
+        if line.strip().startswith("requirements"):
+            requirements_index = i
+            break
+            
+    # Find the index of the line containing the 'queue' command
+    queue_index = len(lines) - 1
+    for i, line in enumerate(reversed(lines)):
+        if line.strip().startswith("queue"):
+            queue_index = len(lines) - 1 - i
+            break    
+    
+    # Modify lines
+    if requirements_index is not None:
+        #requirements_line = lines[requirements_index]
+        if opts.use_osg:
+            requirements_line = "requirements = (HAS_SINGULARITY=?=TRUE)" + "&&" + avoid_string + "\n"
+            lines[requirements_index] = requirements_line
+            # Insert the 'require_gpus' line before the 'queue' command
+            lines.insert(queue_index, "require_gpus = Capability >= 3.5\n")
+        else:
+            requirements_line = "requirements = " + avoid_string + "\n"
+            lines[requirements_index] = requirements_line
+            
+    # Write the modified contents back to the file
+    with open(full_path, "w") as file:
+        file.writelines(lines)
+
+#
+#

From 849af53d7c7da138aa229bfd25bd42d4c458e8d5 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Tue, 7 May 2024 13:10:03 -0700
Subject: [PATCH 03/50] further changes for hyperbolic analysis

---
 .../Code/bin/convert_output_format_ile2inference | 16 ++++++++++++++++
 .../Code/bin/helper_LDG_Events.py                | 10 +++++-----
 .../test/single_inject/single_inject_RIFT.py     | 13 ++++++++++++-
 3 files changed, 33 insertions(+), 6 deletions(-)
 mode change 100644 => 100755 MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py

diff --git a/MonteCarloMarginalizeCode/Code/bin/convert_output_format_ile2inference b/MonteCarloMarginalizeCode/Code/bin/convert_output_format_ile2inference
index 6c306aee4..b0133dd46 100755
--- a/MonteCarloMarginalizeCode/Code/bin/convert_output_format_ile2inference
+++ b/MonteCarloMarginalizeCode/Code/bin/convert_output_format_ile2inference
@@ -39,6 +39,7 @@ optp.add_option("--export-eos",action='store_true',help="Include EOS index param
 optp.add_option("--export-cosmology",action='store_true',help="Include source frame masses and redshift")
 optp.add_option("--export-weights",action='store_true',help="Include a field 'weights' equal to L p/ps")
 optp.add_option("--export-eccentricity", action="store_true", help="Include eccentricity")
+optp.add_option("--export-hyperbolic", action="store_true", help="Include E0 and p_phi0")
 optp.add_option("--with-cosmology",default="Planck15",help="Specific cosmology to use")
 optp.add_option("--use-interpolated-cosmology",action='store_true',help="Specific cosmology to use")
 optp.add_option("--convention",default="RIFT",help="RIFT|LI")
@@ -84,6 +85,8 @@ if opts.convention == 'LI':
       print( " weights ", )
   if opts.export_eccentricity:
       print( "eccentricity ",)
+  if opts.export_hyperbolic:
+      print( "E0 p_phi0 ",)  
   else:
       print('')
   for fname in args:
@@ -96,6 +99,10 @@ if opts.convention == 'LI':
 
     if opts.export_eccentricity:
         ecc = [row.alpha4 for row in points]
+        
+    if opts.export_hyperbolic:
+        p_phi0 = [row.beta for row in points]
+        E0 = [row.psi3 for row in points]
     
     wt = np.exp(like)*p/ps
     neff_here = np.sum(wt)/np.max(wt)  # neff for this file.  Assumes NOT mixed samples: dangerous
@@ -171,6 +178,8 @@ if opts.convention == 'LI':
             print(wt[indx],end=' ')
         if opts.export_eccentricity:
             print(ecc[indx])
+        if opts.export_hyperbolic:
+            print(E0[indx], p_phi0[indx])
         else:
             print('')
 
@@ -192,6 +201,8 @@ if opts.export_weights:
     print( " weights ",end=' ')
 if opts.export_eccentricity:
     print( "eccentricity ",)
+if opts.export_hyperbolic:
+    print( "E0 p_phi0 ",)
 else:
     print('')
 for fname in args:
@@ -233,6 +244,9 @@ for fname in args:
 
     if opts.export_eccentricity:
         ecc = [row.alpha4 for row in points]
+    if opts.export_hyperbolic:
+        p_phi0 = [row.beta for row in points]
+        E0 = [row.psi3 for row in points]
     
     wt = np.exp(like)*p/ps
 
@@ -275,6 +289,8 @@ for fname in args:
             print(wt[indx],end=' ')
         if opts.export_eccentricity:
             print(ecc[indx],)
+        if opts.export_hyperbolic:
+            print(E0[indx], p_phi0[indx],)
         else:
             print('')
 #    print pt.geocent_end_time + 1e-9* pt.geocent_end_time_ns, pt.coa_phase,  pt.inclination, pt.polarization, pt.longitude,pt.latitude, pt.distance, ind like[indx]
diff --git a/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py b/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
index ac80f5994..6f5ba78b6 100755
--- a/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
+++ b/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
@@ -1299,13 +1299,13 @@ def lambda_m_estimate(m):
     if "SNR" in event_dict:
         grid_size *= np.max([1,event_dict["SNR"]/15])  # more grid points at higher amplitude. Yes, even though we also contract the paramete range
 
-    if not (opts.force_initial_grid_size is None):
-        grid_size = opts.force_initial_grid_size
-        
     # hyperbolic grid:
     if opts.assume_hyperbolic:
-        # note - currently uniform for testing purposes
-        cmd += f" --parameter E0 --parameter-range [{opts.E0_min},{opts.E0_max}] --parameter p_phi0 --parameter-range [{opts.pphi0_min},{opts.pphi0_max}] "
+        # note - currently testing
+        cmd += f" --random-parameter E0 --random-parameter-range [{opts.E0_min},{opts.E0_max}] --random-parameter p_phi0 --random-parameter-range [{opts.pphi0_min},{opts.pphi0_max}] "
+
+    if not (opts.force_initial_grid_size is None):
+        grid_size = opts.force_initial_grid_size
         
     cmd += " --grid-cartesian-npts  " + str(int(grid_size))
     print(" Executing grid command ", cmd)
diff --git a/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py b/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py
old mode 100644
new mode 100755
index 97127c45f..c11374895
--- a/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py
+++ b/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py
@@ -481,7 +481,14 @@ def add_frames(channel,input_frame, noise_frame,combined_frame):
     for i, line in enumerate(reversed(lines)):
         if line.strip().startswith("queue"):
             queue_index = len(lines) - 1 - i
-            break    
+            break
+
+    # Find the index of the line containing 'request_memory'
+    memory_index = None
+    for i, line in enumerate(lines):
+        if line.strip().startswith("request_memory"):
+            memory_index = i
+            break
     
     # Modify lines
     if requirements_index is not None:
@@ -494,6 +501,10 @@ def add_frames(channel,input_frame, noise_frame,combined_frame):
         else:
             requirements_line = "requirements = " + avoid_string + "\n"
             lines[requirements_index] = requirements_line
+
+    # Modify the 'request_memory' line if opts.use_hyperbolic is true
+    if opts.use_hyperbolic and memory_index is not None:
+        lines[memory_index] = "request_memory = 16384M\n"
             
     # Write the modified contents back to the file
     with open(full_path, "w") as file:

From c2d455adf41e7c68e1f3bf791164953bc648fc94 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Wed, 8 May 2024 08:31:00 -0700
Subject: [PATCH 04/50] prevent duplicate coords in puffball

---
 MonteCarloMarginalizeCode/Code/bin/util_ParameterPuffball.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/MonteCarloMarginalizeCode/Code/bin/util_ParameterPuffball.py b/MonteCarloMarginalizeCode/Code/bin/util_ParameterPuffball.py
index 30a239b67..ca979919c 100755
--- a/MonteCarloMarginalizeCode/Code/bin/util_ParameterPuffball.py
+++ b/MonteCarloMarginalizeCode/Code/bin/util_ParameterPuffball.py
@@ -60,7 +60,7 @@
     opts.random_parameter = []
 
 # Extract parameter names
-coord_names = opts.parameter # Used  in fit
+coord_names = list(set(opts.parameter)) # Used  in fit
 #if opts.parameter_nofit:
 #    coord_names = coord_names + opts.parameter_nofit
 if coord_names is None:

From 0e6e5c789cb2d08a8f1092ad333ce6f5d543beb3 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Fri, 10 May 2024 06:30:22 -0700
Subject: [PATCH 05/50] fixed memory for ILE_extr

---
 .../Code/test/single_inject/single_inject_RIFT.py           | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py b/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py
index c11374895..6bcf2980f 100755
--- a/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py
+++ b/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py
@@ -454,12 +454,12 @@ def add_frames(channel,input_frame, noise_frame,combined_frame):
 
 # List of file paths to modify
 file_paths = ["ILE.sub","ILE_puff.sub","iteration_4_cip/ILE.sub","iteration_4_cip/ILE_puff.sub"]
-
-if opts.add_extrinsic:
-    file_paths.append("ILE_extr.sub")
     
 if opts.use_hyperbolic:
     file_paths = ["ILE.sub","ILE_puff.sub","iteration_2_cip/ILE.sub","iteration_2_cip/ILE_puff.sub"]
+    
+if opts.add_extrinsic:
+    file_paths.append("ILE_extr.sub")
 
 for file_path in file_paths:
     # Construct the full file path

From 3215dd223379d42f08d2b436f14f4fc623bb3884 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Tue, 14 May 2024 11:26:00 -0700
Subject: [PATCH 06/50] xmlutils + pseudo_pipe event_dict fix

---
 MonteCarloMarginalizeCode/Code/RIFT/misc/xmlutils.py        | 2 ++
 MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py | 2 ++
 2 files changed, 4 insertions(+)

diff --git a/MonteCarloMarginalizeCode/Code/RIFT/misc/xmlutils.py b/MonteCarloMarginalizeCode/Code/RIFT/misc/xmlutils.py
index 140753255..e7f13f970 100644
--- a/MonteCarloMarginalizeCode/Code/RIFT/misc/xmlutils.py
+++ b/MonteCarloMarginalizeCode/Code/RIFT/misc/xmlutils.py
@@ -40,6 +40,8 @@ def assign_time(row, t):
     "joint_prior": "alpha2",
     "joint_s_prior": "alpha3",
     "eccentricity":"alpha4",
+    "E0":"psi3",
+    "p_phi0":"beta",    
     "lambda1":"alpha5",
     "lambda2":"alpha6",
     "spin1x":"spin1x",
diff --git a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py
index 6416bec9e..771f5b8ff 100755
--- a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py
+++ b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py
@@ -93,6 +93,8 @@ def retrieve_event_from_coinc(fname_coinc):
     event_dict["s1z"] = row.spin1z
     event_dict["s2z"] = row.spin2z
     event_dict["eccentricity"] = row.alpha4
+    event_dict["E0"] = row.psi3
+    event_dict["p_phi0"] = row.beta
     event_dict["IFOs"] = list(set(ifo_list))
     max_snr_idx = snr_list.index(max(snr_list))
     event_dict['SNR'] = snr_list[max_snr_idx]

From 79e952131d41a06a198d3916b2ec71ff21639e5d Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Tue, 14 May 2024 11:36:26 -0700
Subject: [PATCH 07/50] xmlutils fix 2

---
 MonteCarloMarginalizeCode/Code/RIFT/misc/xmlutils.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/MonteCarloMarginalizeCode/Code/RIFT/misc/xmlutils.py b/MonteCarloMarginalizeCode/Code/RIFT/misc/xmlutils.py
index e7f13f970..42a26a8c9 100644
--- a/MonteCarloMarginalizeCode/Code/RIFT/misc/xmlutils.py
+++ b/MonteCarloMarginalizeCode/Code/RIFT/misc/xmlutils.py
@@ -54,7 +54,7 @@ def assign_time(row, t):
 
 # FIXME: Find way to intersect given cols with valid cols when making table.
 # Otherwise, we'll have to add them manually and ensure they all exist
-sim_valid_cols = ["simulation_id", "inclination", "longitude", "latitude", "polarization", "geocent_end_time", "geocent_end_time_ns", "coa_phase", "distance", "mass1", "mass2", "alpha1", "alpha2", "alpha3", "alpha4", "alpha5", "alpha6", "spin1x", "spin1y", "spin1z", "spin2x", "spin2y", "spin2z"]
+sim_valid_cols = ["simulation_id", "inclination", "longitude", "latitude", "polarization", "geocent_end_time", "geocent_end_time_ns", "coa_phase", "distance", "mass1", "mass2", "alpha1", "alpha2", "alpha3", "alpha4", "alpha5", "alpha6", "spin1x", "spin1y", "spin1z", "spin2x", "spin2y", "spin2z", "psi3", "beta"]
 sngl_valid_cols = [ "event_id", "snr", "tau0", "tau3"]
 multi_valid_cols = ["process_id", "event_id", "snr"]
 

From 56600c1c041fa5d5d37ebc828d70e2648ebbdcd5 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Wed, 22 May 2024 06:28:51 -0700
Subject: [PATCH 08/50] extrinsic and submit fixes

---
 .../Code/RIFT/misc/xmlutils.py                |  2 +-
 .../Code/bin/helper_LDG_Events.py             | 15 +---
 .../Code/bin/util_RIFT_pseudo_pipe.py         |  4 --
 .../test/single_inject/single_inject_RIFT.py  | 71 +++++++++++++++++++
 4 files changed, 74 insertions(+), 18 deletions(-)

diff --git a/MonteCarloMarginalizeCode/Code/RIFT/misc/xmlutils.py b/MonteCarloMarginalizeCode/Code/RIFT/misc/xmlutils.py
index 1dca463b1..5cd052424 100644
--- a/MonteCarloMarginalizeCode/Code/RIFT/misc/xmlutils.py
+++ b/MonteCarloMarginalizeCode/Code/RIFT/misc/xmlutils.py
@@ -56,7 +56,7 @@ def assign_time(row, t):
 
 # FIXME: Find way to intersect given cols with valid cols when making table.
 # Otherwise, we'll have to add them manually and ensure they all exist
-sim_valid_cols = ["simulation_id", "inclination", "longitude", "latitude", "polarization", "geocent_end_time", "geocent_end_time_ns", "coa_phase", "distance", "mass1", "mass2", "alpha1", "alpha2", "alpha3", "alpha4", "alpha5", "alpha6", "psi3", "beta" "spin1x", "spin1y", "spin1z", "spin2x", "spin2y", "spin2z"]
+sim_valid_cols = ["simulation_id", "inclination", "longitude", "latitude", "polarization", "geocent_end_time", "geocent_end_time_ns", "coa_phase", "distance", "mass1", "mass2", "alpha1", "alpha2", "alpha3", "alpha4", "alpha5", "alpha6", "psi3", "beta", "spin1x", "spin1y", "spin1z", "spin2x", "spin2y", "spin2z"]
 
 sngl_valid_cols = [ "event_id", "snr", "tau0", "tau3"]
 multi_valid_cols = ["process_id", "event_id", "snr"]
diff --git a/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py b/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
index 6a83db1df..5f6ede05e 100755
--- a/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
+++ b/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
@@ -987,10 +987,6 @@ def crit_m2(delta):
         eta_min = q_min/(1.+q_min)**2
     if 'ecc_min' in engine_dict:
         ecc_range_str = "  ["+str(engine_dict['force_ecc_min'])+","+str(engine_dict['force_ecc_max'])+"]"
-    if 'E0_min' in engine_dict:
-        E0_range_str = "  ["+str(engine_dict['force_E0_min'])+","+str(engine_dict['force_E0_max'])+"]"
-    if 'pphi0_min' in engine_dict:
-        pphi0_range_str = "  ["+str(engine_dict['force_pphi0_min'])+","+str(engine_dict['force_pphi0_max'])+"]"
         
 mc_range_str = "  ["+str(mc_min_tight)+","+str(mc_max_tight)+"]"  # Use a tight placement grid for CIP
 if not(opts.manual_mc_min is None):
@@ -1203,8 +1199,7 @@ def crit_m2(delta):
     if opts.assume_eccentric:
         cmd += " --random-parameter eccentricity --random-parameter-range " + ecc_range_str
     if opts.assume_hyperbolic:
-        cmd += " --random-parameter E0 --random-parameter-range " + E0_range_str
-        cmd += " --random-parameter p_phi0 --random-parameter-range " + pphi0_range_str
+        cmd += f" --random-parameter E0 --random-parameter-range [{opts.E0_min},{opts.E0_max}] --random-parameter p_phi0 --random-parameter-range [{opts.pphi0_min},{opts.pphi0_max}] "
     if "SNR" in event_dict:
         grid_size *= np.max([1,event_dict["SNR"]/15])  # more grid points at higher amplitude. Yes, even though we also contract the paramete range
     if not (opts.force_initial_grid_size is None):
@@ -1260,8 +1255,7 @@ def crit_m2(delta):
     if opts.assume_eccentric:
         cmd += " --random-parameter eccentricity --random-parameter-range " + ecc_range_str
     if opts.assume_hyperbolic:
-        cmd += " --random-parameter E0 --random-parameter-range " + E0_range_str
-        cmd += " --random-parameter p_phi0 --random-parameter-range " + pphi0_range_str
+        cmd += f" --random-parameter E0 --random-parameter-range [{opts.E0_min},{opts.E0_max}] --random-parameter p_phi0 --random-parameter-range [{opts.pphi0_min},{opts.pphi0_max}] "
     if opts.internal_tabular_eos_file:
         cmd += " --tabular-eos-file {} ".format(opts.internal_tabular_eos_file)
         grid_size *=2  # larger grids needed for discrete realization scenarios
@@ -1309,11 +1303,6 @@ def lambda_m_estimate(m):
     if "SNR" in event_dict:
         grid_size *= np.max([1,event_dict["SNR"]/15])  # more grid points at higher amplitude. Yes, even though we also contract the paramete range
 
-    # hyperbolic grid:
-    if opts.assume_hyperbolic:
-        # note - currently testing
-        cmd += f" --random-parameter E0 --random-parameter-range [{opts.E0_min},{opts.E0_max}] --random-parameter p_phi0 --random-parameter-range [{opts.pphi0_min},{opts.pphi0_max}] "
-
     if not (opts.force_initial_grid_size is None):
         grid_size = opts.force_initial_grid_size
         
diff --git a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py
index 3fca4427c..ce09ba67e 100755
--- a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py
+++ b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py
@@ -169,10 +169,6 @@ def unsafe_parse_arg_string_dict(my_argstr):
 parser.add_argument("--internal-tabular-eos-file",type=str,default=None,help="Tabular file of EOS to use.  The default prior will be UNIFORM in this table!")
 parser.add_argument("--assume-eccentric",action='store_true', help="Add eccentric options for each part of analysis")
 parser.add_argument("--assume-hyperbolic",action='store_true', help="Add hyperbolic options for each part of analysis")
-parser.add_argument("--force-E0-max",default=None,type=float,help="Provide this value to override the value of E0-max provided")
-parser.add_argument("--force-E0-min",default=None,type=float,help="Provide this value to override the value of E0-min provided")
-parser.add_argument("--force-pphi0-max",default=None,type=float,help="Provide this value to override the value of p_phi0-max provided")
-parser.add_argument("--force-pphi0-min",default=None,type=float,help="Provide this value to override the value of p_phi0-min provided")
 parser.add_argument("--assume-lowlatency-tradeoffs",action='store_true', help="Force analysis with various low-latency tradeoffs (e.g., drop spin 2, use aligned, etc)")
 parser.add_argument("--assume-highq",action='store_true', help="Force analysis with the high-q strategy, neglecting spin2. Passed to 'helper'")
 parser.add_argument("--assume-well-placed",action='store_true',help="If present, the code will adopt a strategy that assumes the initial grid is very well placed, and will minimize the number of early iterations performed. Not as extrme as --propose-flat-strategy")
diff --git a/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py b/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py
index 6bcf2980f..8353716bc 100755
--- a/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py
+++ b/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py
@@ -512,3 +512,74 @@ def add_frames(channel,input_frame, noise_frame,combined_frame):
 
 #
 #
+## special configuration for hyperbolic analysis in the IGWN pool
+if opts.use_hyperbolic and opts.use_osg:
+    
+    print('Modifying submit files for hyperbolics in the IGWN pool')
+    
+    # need to edit all ILE and CIP submit files
+    ILE_file_paths = ["ILE.sub","ILE_puff.sub","iteration_2_cip/ILE.sub","iteration_2_cip/ILE_puff.sub"]
+    CIP_file_paths = ["CIP.sub", "CIP_worker.sub", "CIP_0.sub", "CIP_worker0.sub", "CIP_1.sub", "CIP_worker1.sub", "CIP_2.sub", "CIP_worker2.sub", "iteration_2_cip/CIP.sub", "iteration_2_cip/CIP_worker.sub"]
+    
+
+    if opts.add_extrinsic:
+        ILE_file_paths.append("ILE_extr.sub")
+        
+        
+    total_file_paths = ILE_file_paths + CIP_file_paths
+    
+    
+    def modify_submit_file(input_file):
+        
+        current_singularity_image = 'rift_o4b_jl-chadhenshaw-teobresums_eccentric-2024-05-14_12-02-56.sif'
+        sif_path = 'osdf:///igwn/cit/staging/james.clark/' + current_singularity_image
+        
+        with open(input_file, 'r') as file:
+            lines = file.readlines()
+    
+        lines_to_remove = ['getenv', '+SingularityBindCVMFS', '+SingularityImage']
+        lines_to_add = [f'MY.SingularityImage = "{current_singularity_image}"', 'use_oauth_services = scitokens']
+        modifications = {'request_disk': 'request_disk = 3G', 'requirements': 'requirements = (HAS_SINGULARITY=?=TRUE)&&(IS_GLIDEIN=?=TRUE)', 'transfer_input_files': None}
+
+        modified_lines = []
+        queue_line = None
+
+        for line in lines:
+            stripped_line = line.strip()
+
+            # remove lines
+            if any(stripped_line.startswith(keyword) for keyword in lines_to_remove):
+                continue
+
+            # modify lines
+            modified = False
+            for key, new_value in modifications.items():
+                if stripped_line.startswith(key):
+                    if key == 'transfer_input_files':
+                        line = stripped_line + ',' + sif_path + '\n'
+                    else:
+                        line = new_value + '\n'
+                    modified_lines.append(line)
+                    modified = True
+                    break
+                    
+            if stripped_line.startswith('queue'):
+                queue_line = line
+            elif not modified:
+                modified_lines.append(line)
+                
+
+        # add new lines
+        modified_lines.extend([line + '\n' for line in lines_to_add])
+        
+        if queue_line:
+            modified_lines.append(queue_line)
+
+        with open(input_file, 'w') as file:
+            file.writelines(modified_lines)
+            
+        return
+    
+    
+    for submit_file in total_file_paths:
+        modify_submit_file(submit_file)
\ No newline at end of file

From 7e28d3f530c03ea4ab1235f0c753a70a22edd675 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Tue, 28 May 2024 10:58:17 -0700
Subject: [PATCH 09/50] peak finding for scatter/capture/plunge tapering

---
 .../Code/RIFT/lalsimutils.py                  | 57 ++++++++++++++++---
 1 file changed, 48 insertions(+), 9 deletions(-)

diff --git a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
index f7ce8c1dd..34cf02e19 100644
--- a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
+++ b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
@@ -3020,10 +3020,30 @@ def hoft(P, Fp=None, Fc=None,**kwargs):
             ht.data.data[0:n_samp] *= vectaper
             # If Scattering waveform, tapers same amount as early taper
             # If Capture waveform, no end taper
-            if np.abs(ht.data.length-nmax) > 3e3:
-                n_samp2=n_samp
-                vectaper2= 0.5 - 0.5*np.cos(np.pi* (1-np.arange(n_samp2+1)/(1.*n_samp2)))
-                ht.data.data[-(n_samp2+1):] *= vectaper2
+            
+            # peak finding to determine system type
+            peaks, props = signal.find_peaks(ht.data.data, height = 0.25*np.abs(np.amax(ht.data.data)))
+            peak_heights = props['peak_heights']
+            indices_to_keep = set()
+            sorted_indices = np.argsort(peak_heights)[::-1]
+            tol = 300 # hardcoded peak spacing tolerance. if spacing is less than tol, discard the peak
+            for i in sorted_indices:
+                peak = peaks[i]                
+                keep = True
+                for kept_index in indices_to_keep:
+                    if abs(peaks[kept_index] - peak) <= tol:
+                        keep = False
+                        break                
+                if keep:
+                    indices_to_keep.add(i)                    
+            filtered_peaks = peaks[list(indices_to_keep)]
+            
+            if len(filtered_peaks) == 1:
+                # check if a scatter or a plunge            
+                if np.abs(ht.data.length-nmax) > 3e3:
+                    n_samp2=n_samp
+                    vectaper2= 0.5 - 0.5*np.cos(np.pi* (1-np.arange(n_samp2+1)/(1.*n_samp2)))
+                    ht.data.data[-(n_samp2+1):] *= vectaper2
                 
     if P.deltaF is not None:
         TDlen = int(1./P.deltaF * 1./P.deltaT)
@@ -3654,11 +3674,30 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil
                 hlm[mode].data.data[0:n_samp] *= vectaper
             # If Capture waveform, no end taper
             # If Scattering waveform, tapers end same amount as early taper
-            if np.abs(hlm[(2,2)].data.length-nmax) > 3e3:
-                n_samp2=n_samp
-                vectaper2= 0.5 - 0.5*np.cos(np.pi* (1-np.arange(n_samp2+1)/(1.*n_samp2)))
-                for mode in modes_used_new2:
-                    hlm[mode].data.data[-(n_samp2+1):] *= vectaper2
+            # peak finding to determine system type
+            peaks, props = signal.find_peaks(hlm[(2,2)].data.data, height = 0.25*np.amax(hlm[(2,2)].data.data))
+            peak_heights = props['peak_heights']
+            indices_to_keep = set()
+            sorted_indices = np.argsort(peak_heights)[::-1]
+            tol = 300 # hardcoded tolerance
+            for i in sorted_indices:
+                peak = peaks[i]                
+                keep = True
+                for kept_index in indices_to_keep:
+                    if abs(peaks[kept_index] - peak) <= tol:
+                        keep = False
+                        break                
+                if keep:
+                    indices_to_keep.add(i)                    
+            filtered_peaks = peaks[list(indices_to_keep)]
+            
+            if len(filtered_peaks) == 1:
+                #check for scatter waveform            
+                if np.abs(hlm[(2,2)].data.length-nmax) > 3e3:
+                    n_samp2=n_samp
+                    vectaper2= 0.5 - 0.5*np.cos(np.pi* (1-np.arange(n_samp2+1)/(1.*n_samp2)))
+                    for mode in modes_used_new2:
+                        hlm[mode].data.data[-(n_samp2+1):] *= vectaper2
         for mode in modes_used_new2:
             if not (P.deltaF is None):
                 TDlen = int(1./P.deltaF * 1./P.deltaT)

From bfe7531cdc89f9802e6a65406ddf129dff06dbd6 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Tue, 4 Jun 2024 13:21:07 -0700
Subject: [PATCH 10/50] hyperbolic params in coinc

---
 .../Code/bin/util_SimInspiralToCoinc.py                    | 2 ++
 .../Code/test/single_inject/single_inject_RIFT.py          | 7 +++++--
 2 files changed, 7 insertions(+), 2 deletions(-)

diff --git a/MonteCarloMarginalizeCode/Code/bin/util_SimInspiralToCoinc.py b/MonteCarloMarginalizeCode/Code/bin/util_SimInspiralToCoinc.py
index 8617d046e..be2ac6d56 100755
--- a/MonteCarloMarginalizeCode/Code/bin/util_SimInspiralToCoinc.py
+++ b/MonteCarloMarginalizeCode/Code/bin/util_SimInspiralToCoinc.py
@@ -126,6 +126,8 @@ def _empty_row(obj):
         sngl.snr = 20.  # made up, needed for some algorithms to work
 
     sngl.alpha4 = P.eccentricity
+    sngl.psi3 = P.E0
+    sngl.beta = P.p_phi0
     # add to table
     sngl_table.append(sngl)
 
diff --git a/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py b/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py
index 8353716bc..00436f2c5 100755
--- a/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py
+++ b/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py
@@ -531,8 +531,11 @@ def add_frames(channel,input_frame, noise_frame,combined_frame):
     
     def modify_submit_file(input_file):
         
-        current_singularity_image = 'rift_o4b_jl-chadhenshaw-teobresums_eccentric-2024-05-14_12-02-56.sif'
-        sif_path = 'osdf:///igwn/cit/staging/james.clark/' + current_singularity_image
+        #current_singularity_image = 'rift_o4b_jl-chadhenshaw-teobresums_eccentric-2024-05-14_12-02-56.sif'
+        #sif_path = 'osdf:///igwn/cit/staging/james.clark/' + current_singularity_image
+        
+        current_singularity_image = 'rift_o4b_jl-2024-05-29_09-54-05.sif'
+        sif_path = 'osdf:///igwn/cit/staging/chad.henshaw/' + current_singularity_image
         
         with open(input_file, 'r') as file:
             lines = file.readlines()

From 30406fab4b048d26e74912220cc01fe2c7bebb01 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Wed, 5 Jun 2024 05:14:08 -0700
Subject: [PATCH 11/50] removed sys exit when EOBRun_module not found

---
 MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
index 34cf02e19..4d8967b5d 100644
--- a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
+++ b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
@@ -34,7 +34,6 @@
     has_external_teobresum=True
  except:
     print('Failed to import EOBRun_module')
-    sys.exit(1)
     has_external_teobresum=False
     True; # print(" - no EOBRun (TEOBResumS) - ")
 info_use_resum_polarizations = False

From cdd680375bbc955d00ea67641f31992712668fd9 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Wed, 5 Jun 2024 12:13:11 -0700
Subject: [PATCH 12/50] removed problematic print statement

---
 MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
index 4d8967b5d..9bff20a00 100644
--- a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
+++ b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
@@ -33,7 +33,7 @@
     import EOBRun_module
     has_external_teobresum=True
  except:
-    print('Failed to import EOBRun_module')
+    #print('Failed to import EOBRun_module')
     has_external_teobresum=False
     True; # print(" - no EOBRun (TEOBResumS) - ")
 info_use_resum_polarizations = False

From b45c476b6777ee040f99bd6546ed01b03cc8d78b Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Tue, 16 Jul 2024 09:22:49 -0700
Subject: [PATCH 13/50] hyperbolic epoch fix

---
 .../Code/RIFT/lalsimutils.py                  | 76 ++++++++++++++++++-
 .../test/single_inject/single_inject_RIFT.py  |  2 +-
 2 files changed, 75 insertions(+), 3 deletions(-)

diff --git a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
index 9bff20a00..092737f89 100644
--- a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
+++ b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
@@ -2879,6 +2879,7 @@ def hoft(P, Fp=None, Fc=None,**kwargs):
             }
         elif (P.eccentricity == 0.0):
             print("Using hyperbolic call RIFT O4b branch")
+            hyp_wav = True # convenient way to know if the waveform is hyperbolic
             pars = {
                 'M'                  : M1+M2,
                 'q'                  : M1/M2,
@@ -2933,7 +2934,42 @@ def hoft(P, Fp=None, Fc=None,**kwargs):
         print("Starting EOBRun_module")
         t, hptmp, hctmp, hlmtmp, dyn = EOBRun_module.EOBRunPy(pars)
         print("EOBRun_module done")
-        hpepoch = -P.deltaT*np.argmax(np.abs(hptmp)**2+np.abs(hctmp)**2)
+        
+        if not(hyp_wav):
+            ## Set the epoch for non-hyperbolic cases ##
+            hpepoch = -P.deltaT*np.argmax(np.abs(hptmp)**2+np.abs(hctmp)**2)            
+        else:
+            ## custom epoch for the hyperbolic case ##            
+            # wf amplitude
+            amp = np.abs(hptmp)**2+np.abs(hctmp)**2
+            amp_times = P.deltaT*np.arange(len(amp))
+            amp_max_ind = np.argmax(amp)        
+            amp_norm = amp / amp[amp_max_ind] # normalize amplitude for peak finding                
+            # peak finding to determine system type
+            peaks, props = signal.find_peaks(amp_norm, height = 0.25)    
+            peak_heights = props['peak_heights']
+            # filtering out peaks so we only keep the local maxima
+            indices_to_keep = set()
+            sorted_indices = np.argsort(peak_heights)[::-1]
+            tol = int(pars['srate_interp'] / 13.65) # 300 samples at srate of 4096 - minimum distance between peaks.
+            for i in sorted_indices:
+                peak = peaks[i]                
+                keep = True
+                for kept_index in indices_to_keep:
+                    if abs(peaks[kept_index] - peak) <= tol:
+                        keep = False
+                        break                
+                if keep:
+                    indices_to_keep.add(i)                    
+            filtered_peaks = peaks[list(indices_to_keep)]
+            # parsing number of peaks after filtering against distance tolerance
+            if len(filtered_peaks) == 1:
+                # scatter case OR plunge case, we can set the epoch normally
+                hpepoch = -P.deltaT*np.argmax(np.abs(hptmp)**2+np.abs(hctmp)**2)
+            else:
+                # capture case, we need to force the epoch to be the last peak
+                hpepoch = -P.deltaT*filtered_peaks[-1]
+                
         hplen = len(hptmp)
         hp = {}
         hc = {}
@@ -3488,6 +3524,7 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil
             
         elif (P.eccentricity == 0.0):
             print("Using hyperbolic call RIFT O4b branch")
+            hyp_wav = True # convenient way to know if the waveform is hyperbolic
             pars = {
                 'M'                  : M1+M2,
                 'q'                  : M1/M2,
@@ -3548,7 +3585,42 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil
         t, hptmp, hctmp, hlmtmp, dym = EOBRun_module.EOBRunPy(pars)
         print("EOBRun_module done")
         k_list_orig = hlmtmp.keys()
-        hpepoch = -P.deltaT*np.argmax(np.abs(hptmp)**2+np.abs(hctmp)**2)
+        
+        if not(hyp_wav):
+            ## Set the epoch for non-hyperbolic cases ##
+            hpepoch = -P.deltaT*np.argmax(np.abs(hptmp)**2+np.abs(hctmp)**2)            
+        else:
+            ## custom epoch for the hyperbolic case ##            
+            # wf amplitude
+            amp = np.abs(hptmp)**2+np.abs(hctmp)**2
+            amp_times = P.deltaT*np.arange(len(amp))
+            amp_max_ind = np.argmax(amp)        
+            amp_norm = amp / amp[amp_max_ind] # normalize amplitude for peak finding                
+            # peak finding to determine system type
+            peaks, props = signal.find_peaks(amp_norm, height = 0.25)    
+            peak_heights = props['peak_heights']
+            # filtering out peaks so we only keep the local maxima
+            indices_to_keep = set()
+            sorted_indices = np.argsort(peak_heights)[::-1]
+            tol = int(pars['srate_interp'] / 13.65) # 300 samples at srate of 4096 - minimum distance between peaks.
+            for i in sorted_indices:
+                peak = peaks[i]                
+                keep = True
+                for kept_index in indices_to_keep:
+                    if abs(peaks[kept_index] - peak) <= tol:
+                        keep = False
+                        break                
+                if keep:
+                    indices_to_keep.add(i)                    
+            filtered_peaks = peaks[list(indices_to_keep)]
+            # parsing number of peaks after filtering against distance tolerance
+            if len(filtered_peaks) == 1:
+                # scatter case OR plunge case, we can set the epoch normally
+                hpepoch = -P.deltaT*np.argmax(np.abs(hptmp)**2+np.abs(hctmp)**2)
+            else:
+                # capture case, we need to force the epoch to be the last peak
+                hpepoch = -P.deltaT*filtered_peaks[-1]
+                
         hlmlen = len(hptmp)
         hlm = {}
         hlmtmp2 = {}
diff --git a/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py b/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py
index 00436f2c5..21687c298 100755
--- a/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py
+++ b/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py
@@ -15,7 +15,7 @@
 ## Basic Code Process: ##
 
 # The code process is as follows, starting from reading in the ini file: (example: single_inject_RIFT.py -- /home/albert.einstein/path/to/config/BBH_inject_example.ini)
-# Create a MDC (mock data challenge) SimInspiral table xml from the waveform parameters
+# Create a MDC (mock data challenge) SimInspiral table xml from the waveform parameterrs
 # uses LALWriteFrame.py to create frame files for the specified detectors, stored in the signal_frames dir (also combined_frames if there is no noise)
 # if specified with --use-noise, adds noise frames and dumps the combined into the combined_frames dir
 # calculates SNR with util_FrameZeroNoiseSNR.py, uses network SNR as the initial guess

From e4ac410867174a1752d41b509f3d2b6aa522b8c9 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Fri, 26 Jul 2024 07:23:57 -0700
Subject: [PATCH 14/50] mtot coords

---
 .../Code/RIFT/lalsimutils.py                  |  1 +
 .../Code/bin/helper_LDG_Events.py             | 34 ++++++++++++++++---
 .../Code/bin/util_RIFT_pseudo_pipe.py         | 16 ++++++++-
 3 files changed, 46 insertions(+), 5 deletions(-)

diff --git a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
index 092737f89..e4fe77889 100644
--- a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
+++ b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
@@ -3495,6 +3495,7 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil
         M1=P.m1/lal.MSUN_SI
         M2=P.m2/lal.MSUN_SI
         nu=M1*M2/((M1+M2)**2)
+        hyp_wav = False
         if (P.eccentricity == 0.0 and P.E0 ==0.0):
             print("Using ResumS master; not eccentric")
             pars = {
diff --git a/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py b/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
index 5f6ede05e..c7fee81af 100755
--- a/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
+++ b/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
@@ -157,6 +157,7 @@ def get_observing_run(t):
 parser.add_argument("--limit-mc-range",default=None,type=str,help="For PP plots, or other analyses requiring a specific mc range (eg ini file), bounding the limit *above*.  Allows the code to auto-select its mc range as usual, then takes the intersection with this limit")
 parser.add_argument("--scale-mc-range",type=float,default=None,help="If using the auto-selected mc, scale the ms range proposed by a constant factor. Recommend > 1. . ini file assignment will override this.")
 parser.add_argument("--force-eta-range",default=None,type=str,help="For PP plots. Enforces initial grid placement inside this region")
+parser.add_argument("--force-mtot-range",default=None,type=str,help="For PP plots, hyperbolic analysis, or other analyses requiring a specific mtot range (eg ini file). Enforces initial grid placement inside this region. Passed directly to MOG and CIP.")
 parser.add_argument("--allow-subsolar", action='store_true', help="Override limits which otherwise prevent subsolar mass PE")
 parser.add_argument("--use-legacy-gracedb",action='store_true')
 parser.add_argument("--event-time",type=float,default=None)
@@ -203,6 +204,7 @@ def get_observing_run(t):
 parser.add_argument("--E0-min", default=1.0,type=float,help="Minimum range of 'E0' allowed.")
 parser.add_argument("--pphi0-max", default=5.4,type=float,help="Maximum range of 'p_phi0' allowed.")
 parser.add_argument("--pphi0-min", default=3.8,type=float,help="Minumum range of 'p_phi0' allowed.")
+parser.add_argument("--use-mtot-coords",action='store_true',help="Configures CIP and PUFF for mtot instead of mc. REQUIRES --force-mtot-range.")
 parser.add_argument("--assume-nospin",action='store_true',help="If present, the code will not add options to manage precessing spins (the default is aligned spin)")
 parser.add_argument("--assume-precessing-spin",action='store_true',help="If present, the code will add options to manage precessing spins (the default is aligned spin)")
 parser.add_argument("--assume-volumetric-spin",action='store_true',help="If present, the code will assume a volumetric spin prior in its last iterations. If *not* present, the code will adopt a uniform magnitude spin prior in its last iterations. If not present, generally more iterations are taken.")
@@ -252,6 +254,12 @@ def get_observing_run(t):
 parser.add_argument("--verbose",action='store_true')
 opts=  parser.parse_args()
 
+if opts.use_mtot_coords:
+    if opts.force_mtot_range is None:
+        print('Using the mtot coords requires a specified range!')
+        print('Specify the mtot range with --force-mtot-range!')
+        sys.exit(1)
+
 if opts.assume_matter_but_primary_bh:
     opts.assume_matter=True
 
@@ -996,6 +1004,9 @@ def crit_m2(delta):
 mc_range_str_cip = " --mc-range ["+str(mc_min)+","+str(mc_max)+"]"
 if not(opts.force_mc_range is None):
     mc_range_str_cip = " --mc-range " + opts.force_mc_range
+if not(opts.force_mtot_range is None):
+    mtot_range_str = opts.force_mtot_range
+    mtot_range_str_cip = " --mtot-range " + opts.force_mtot_range
 elif opts.limit_mc_range:
     line_here = list(map(float,opts.limit_mc_range.replace('[','').replace(']','').split(',') ))
     mc_min_lim,mc_max_lim = line_here
@@ -1216,13 +1227,18 @@ def crit_m2(delta):
     # add basic mass parameters
     cmd  = "util_ManualOverlapGrid.py  --fname proposed-grid --skip-overlap "
     mass_string_init = " --random-parameter mc --random-parameter-range   " + mc_range_str + "  --random-parameter delta_mc --random-parameter-range '[" + str(delta_grid_min) +"," + str(delta_grid_max) + "]'  "
+    if not(opts.force_mtot_range is None):
+        mass_string_init = " --random-parameter mtot --random-parameter-range   " + mtot_range_str + "  --random-parameter delta_mc --random-parameter-range '[" + str(delta_grid_min) +"," + str(delta_grid_max) + "]'  "
     cmd+= mass_string_init
     # Add standard downselects : do not have m1, m2 be less than 1
     if not(opts.force_mc_range is None):
         # force downselect based on this range
         cmd += " --downselect-parameter mc --downselect-parameter-range " + opts.force_mc_range 
     if not(opts.force_eta_range is None):
-        cmd += " --downselect-parameter eta --downselect-parameter-range " + opts.force_eta_range 
+        cmd += " --downselect-parameter eta --downselect-parameter-range " + opts.force_eta_range
+    if not(opts.force_mtot_range is None):
+        # force downselect based on this range
+        cmd += " --downselect-parameter mtot --downselect-parameter-range " + opts.force_mtot_range
     cmd += " --fmin " + str(opts.fmin_template)
     if opts.data_LI_seglen and not (opts.no_enforce_duration_bound):  
         cmd += " --enforce-duration-bound " + str(opts.data_LI_seglen)
@@ -1400,6 +1416,8 @@ def lambda_m_estimate(m):
 
 puff_max_it=0
 helper_puff_args = " --parameter mc --parameter eta --fmin {} --fref {} ".format(opts.fmin_template,opts.fmin_template)
+if opts.use_mtot_coords:
+    helper_puff_args = " --parameter mtot --parameter q --fmin {} --fref {} ".format(opts.fmin_template,opts.fmin_template)
 if opts.assume_eccentric:
     helper_puff_args += " --parameter eccentricity "
 if opts.assume_hyperbolic:
@@ -1407,7 +1425,9 @@ def lambda_m_estimate(m):
 
 if event_dict["MChirp"] >25:
     # at high mass, mc/eta correlation weak, don't want to have eta coordinate degeneracy at q=1 to reduce puff proposals  near there
-    helper_puff_args = " --parameter mc --parameter delta_mc "  
+    helper_puff_args = " --parameter mc --parameter delta_mc "
+    if opts.use_mtot_coords:
+        helper_puff_args = " --parameter mtot --parameter delta_mc "
 if opts.propose_fit_strategy:
     puff_max_it= 0
     # Strategy: One iteration of low-dimensional, followed by other dimensions of high-dimensional
@@ -1415,13 +1435,19 @@ def lambda_m_estimate(m):
     lnLoffset_late = 15 # default
     helper_cip_args += ' --no-plots --fit-method {}  '.format(fit_method)
     if not opts.internal_use_aligned_phase_coordinates:
-        helper_cip_args += '   --parameter mc --parameter delta_mc '
+        if not(opts.use_mtot_coords):
+            helper_cip_args += '   --parameter mc --parameter delta_mc '
+        else:
+            helper_cip_args += '   --parameter mtot --parameter delta_mc '
     else:
         helper_cip_args += " --parameter-implied mu1 --parameter-implied mu2 --parameter-nofit mc --parameter delta_mc "  
     if 'gp' in fit_method:
         helper_cip_args += " --cap-points 12000 "
     if not opts.no_propose_limits:
-        helper_cip_args += mc_range_str_cip + eta_range_str_cip
+        if not(opts.use_mtot_coords):
+            helper_cip_args += mc_range_str_cip + eta_range_str_cip
+        else:
+            helper_cip_args += mtot_range_str_cip + eta_range_str_cip
     if opts.force_chi_max:
         helper_cip_args += " --chi-max {} ".format(opts.force_chi_max)
     if opts.force_chi_small_max:
diff --git a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py
index ce09ba67e..3e6cd6f53 100755
--- a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py
+++ b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py
@@ -228,6 +228,7 @@ def unsafe_parse_arg_string_dict(my_argstr):
 parser.add_argument("--limit-mc-range",default=None,type=str,help="Pass this argumen through to the helper to set the mc range")
 parser.add_argument("--force-mc-range",default=None,type=str,help="Pass this argumen through to the helper to set the mc range")
 parser.add_argument("--force-eta-range",default=None,type=str,help="Pass this argumen through to the helper to set the eta range")
+parser.add_argument("--force-mtot-range",default=None,type=str,help="Pass this argument through to the helper to set the mtot range. Overrides mc parameter with mtot parameter broadly throughout the pipeline.")
 parser.add_argument("--force-comp-max",default=1000,type=float,help="Provde this value to override the value of component mass in CIP provided")
 parser.add_argument("--force-comp-min",default=1,type=float,help="Provde this value to override the value of component mass in CIP provided")
 parser.add_argument("--allow-subsolar", action='store_true', help="Override limits which otherwise prevent subsolar mass PE")
@@ -287,6 +288,8 @@ def unsafe_parse_arg_string_dict(my_argstr):
 parser.add_argument("--archive-pesummary-label",default=None,help="If provided, creates a 'pesummary' directory and fills it with this run's final output at the end of the run")
 parser.add_argument("--archive-pesummary-event-label",default="this_event",help="Label to use on the pesummary page itself")
 parser.add_argument("--internal-mitigate-fd-J-frame",default="L_frame",help="L_frame|rotate, choose method to deal with ChooseFDWaveform being in wrong frame. Default is to request L frame for inputs")
+parser.add_argument("--first-iteration-jumpstart",action='store_true',help="No ILE jobs the first iteration.  Assumes you already have .composite files and want to get going. Particularly helpful for subdag systems")
+parser.add_argument("--use-mtot-coords",action='store_true',help="Passed to the helper to configure CIP and PUFF for mtot instead of mc.")
 opts=  parser.parse_args()
 
 
@@ -562,6 +565,10 @@ def unsafe_parse_arg_string_dict(my_argstr):
     # default value for eccentricity is 0 for 'P'!  Only change this value from default if eccentricity is present, do NOT want to fill it with None in particular
     if not(event_dict['eccentricity'] is None):   
         P.eccentricity = event_dict["eccentricity"]
+    # same for the hyperbolic params
+    if not(event_dict['E0'] is None):
+        P.E0 = event_dict['E0']
+        P.p_phi0 = event_dict["p_phi0"]
     # Write 'target_params.xml.gz' file
     lalsimutils.ChooseWaveformParams_array_to_xml([P], "target_params")
 
@@ -680,6 +687,8 @@ def unsafe_parse_arg_string_dict(my_argstr):
     cmd += " --scale-mc-range  " + str(opts.scale_mc_range).replace(' ','')
 if not(opts.force_eta_range is None):
     cmd+= " --force-eta-range  " + str(opts.force_eta_range).replace(' ','')
+if not(opts.force_mtot_range is None):
+    cmd+= " --force-mtot-range  " + str(opts.force_mtot_range).replace(' ','')
 if opts.allow_subsolar:
     cmd += " --allow-subsolar "
 if opts.force_chi_max:
@@ -1076,7 +1085,10 @@ def unsafe_parse_arg_string_dict(my_argstr):
             line += " --ecc-min {}  ".format(ecc_min)
     if opts.assume_hyperbolic:
         if not(opts.internal_use_aligned_phase_coordinates):
-            line = line.replace('parameter mc', 'parameter mc --parameter E0 --parameter p_phi0 --use-hyperbolic')
+            if not(opts.use_mtot_coords):
+                line = line.replace('parameter mc', 'parameter mc --parameter E0 --parameter p_phi0 --use-hyperbolic')
+            else:
+                line = line.replace('parameter mtot', 'parameter mtot --parameter E0 --parameter p_phi0 --use-hyperbolic')
         else:
             line = line.replace('parameter-nofit mc', 'parameter-nofit mc --parameter E0 --parameter p_phi0 --use-hyperbolic')
         if not(opts.force_E0_max is None):
@@ -1281,6 +1293,8 @@ def unsafe_parse_arg_string_dict(my_argstr):
     cmd += " --comov-distance-reweighting --comov-distance-reweighting-exe `which make_uni_comov_skymap.py` --convert-ascii2h5-exe `which convert_output_format_ascii2h5.py` "
 if opts.use_gauss_early:
     cmd += " --cip-exe-G `which util_ConstructIntrinsicPosterior_GaussianResampling.py ` "
+if opts.first_iteration_jumpstart:
+    cmd += " --first-iteration-jumpstart "
 if opts.internal_use_amr:
     print(" AMR prototype: Using hardcoded aligned-spin settings, assembling grid, requires coinc!")
     cmd += " --cip-exe `which util_AMRGrid.py ` "

From 0f81d7f2a0b3a9226126f42cfc352f6ca43f3c80 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Fri, 26 Jul 2024 07:31:32 -0700
Subject: [PATCH 15/50] alt start and mtot edits

---
 .../test/single_inject/single_inject_RIFT.py  | 104 ++++++++++++++++--
 1 file changed, 96 insertions(+), 8 deletions(-)

diff --git a/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py b/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py
index 21687c298..b2df1c1a7 100755
--- a/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py
+++ b/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py
@@ -118,6 +118,7 @@ def add_frames(channel,input_frame, noise_frame,combined_frame):
 parser.add_argument("--bypass-frames",action='store_true',help="Skip making mdc and frame files, use with caution")
 parser.add_argument("--just-frames",action='store_true',help="Stop after making frame files, use with caution")
 parser.add_argument("--use-hyperbolic",action='store_true',help="Adds hyperbolic options, requires TEOBResumS approx. NOTE: development still in progress.")
+parser.add_argument("--alternate-start",action='store_true',help="Currently just passses first-iteration-jumpstart to start from CIP")
 opts =  parser.parse_args()
 
 config = configparser.ConfigParser(allow_no_value=True) #SafeConfigParser deprecated from py3.2
@@ -393,6 +394,13 @@ def add_frames(channel,input_frame, noise_frame,combined_frame):
 if opts.use_hyperbolic:
     print('Adding hyperbolic options')
     cmd += ' --assume-hyperbolic '
+    
+if opts.alternate_start:
+    # ensure there is first iteration jumpstart
+    if config.getboolean('rift-pseudo-pipe', 'first-iteration-jumpstart'):
+        pass
+    else:
+        cmd += ' --first-iteration-jumpstart '
 
     
 os.system(cmd)   
@@ -433,11 +441,38 @@ def add_frames(channel,input_frame, noise_frame,combined_frame):
         shutil.copy2('distance_marginalization_lookup.npz', working_dir_full)
     
 
-    
-# change name of dag
+# edit dag file + change name of dag
 os.chdir(working_dir_full+"/"+rundir)
+if opts.alternate_start:
+    with open('marginalize_intrinsic_parameters_BasicIterationWorkflow.dag', 'r') as file:
+        lines = file.readlines()
+        
+    # remove first 6 lines
+    lines = lines[6:]
+    
+    # Find the start of the parent/child section
+    parent_child_start = None
+    for i, line in enumerate(lines):
+        if line.startswith('PARENT'):
+            parent_child_start = i
+            break
+    
+    if parent_child_start is not None:
+        # Remove the first 2 lines in the parent/child section
+        lines = lines[:parent_child_start] + lines[parent_child_start+2:]
+
+    # Write the modified lines back to the file
+    with open('marginalize_intrinsic_parameters_BasicIterationWorkflow.dag', 'w') as file:
+        file.writelines(lines)
+
+
 os.rename('marginalize_intrinsic_parameters_BasicIterationWorkflow.dag', working_dir + '.dag')
 
+# For alternate start, copy over all.net file
+if opts.alternate_start:
+    net_file = os.path.join(working_dir_full, 'analysis_0/all.net')
+    shutil.copy2(net_file, working_dir_full+"/"+rundir+"/all.net")
+
 
 ## Update ILE requirements to avoid bad hosts ##
 
@@ -455,8 +490,11 @@ def add_frames(channel,input_frame, noise_frame,combined_frame):
 # List of file paths to modify
 file_paths = ["ILE.sub","ILE_puff.sub","iteration_4_cip/ILE.sub","iteration_4_cip/ILE_puff.sub"]
     
-if opts.use_hyperbolic:
-    file_paths = ["ILE.sub","ILE_puff.sub","iteration_2_cip/ILE.sub","iteration_2_cip/ILE_puff.sub"]
+if opts.use_hyperbolic:    
+    if config.getboolean('rift-pseudo-pipe', 'internal-use-aligned-phase-coordinates'):  
+        file_paths = ["ILE.sub","ILE_puff.sub","iteration_2_cip/ILE.sub","iteration_2_cip/ILE_puff.sub"]
+    else:
+        file_paths = ["ILE.sub","ILE_puff.sub","iteration_3_cip/ILE.sub","iteration_3_cip/ILE_puff.sub"]
     
 if opts.add_extrinsic:
     file_paths.append("ILE_extr.sub")
@@ -518,8 +556,14 @@ def add_frames(channel,input_frame, noise_frame,combined_frame):
     print('Modifying submit files for hyperbolics in the IGWN pool')
     
     # need to edit all ILE and CIP submit files
-    ILE_file_paths = ["ILE.sub","ILE_puff.sub","iteration_2_cip/ILE.sub","iteration_2_cip/ILE_puff.sub"]
-    CIP_file_paths = ["CIP.sub", "CIP_worker.sub", "CIP_0.sub", "CIP_worker0.sub", "CIP_1.sub", "CIP_worker1.sub", "CIP_2.sub", "CIP_worker2.sub", "iteration_2_cip/CIP.sub", "iteration_2_cip/CIP_worker.sub"]
+    if config.getboolean('rift-pseudo-pipe', 'internal-use-aligned-phase-coordinates'):
+        ILE_file_paths = ["ILE.sub","ILE_puff.sub","iteration_2_cip/ILE.sub","iteration_2_cip/ILE_puff.sub"]
+        CIP_file_paths = ["CIP.sub", "CIP_worker.sub", "CIP_0.sub", "CIP_worker0.sub", "CIP_1.sub", "CIP_worker1.sub", "CIP_2.sub", "CIP_worker2.sub", "iteration_2_cip/CIP.sub", "iteration_2_cip/CIP_worker.sub"]
+        PUFF_file_paths = ["PUFF.sub", "iteration_2_cip/PUFF.sub"]
+    else:
+        ILE_file_paths = ["ILE.sub","ILE_puff.sub","iteration_3_cip/ILE.sub","iteration_3_cip/ILE_puff.sub"]
+        CIP_file_paths = ["CIP.sub", "CIP_worker.sub", "CIP_0.sub", "CIP_worker0.sub", "CIP_1.sub", "CIP_worker1.sub", "CIP_2.sub", "CIP_worker2.sub", "iteration_3_cip/CIP.sub", "iteration_3_cip/CIP_worker.sub"]
+        PUFF_file_paths = ["PUFF.sub", "iteration_3_cip/PUFF.sub"]
     
 
     if opts.add_extrinsic:
@@ -534,7 +578,7 @@ def modify_submit_file(input_file):
         #current_singularity_image = 'rift_o4b_jl-chadhenshaw-teobresums_eccentric-2024-05-14_12-02-56.sif'
         #sif_path = 'osdf:///igwn/cit/staging/james.clark/' + current_singularity_image
         
-        current_singularity_image = 'rift_o4b_jl-2024-05-29_09-54-05.sif'
+        current_singularity_image = 'rift_o4b_jl-2024-07-16_12-21-57.sif'
         sif_path = 'osdf:///igwn/cit/staging/chad.henshaw/' + current_singularity_image
         
         with open(input_file, 'r') as file:
@@ -583,6 +627,50 @@ def modify_submit_file(input_file):
             
         return
     
+    def modify_CIP_coords(input_file, mtot_range):
+        with open(input_file, 'r') as file:
+            lines = file.readlines()
+
+        # Format the new mtot_range string
+        mtot_range_str = f"[{mtot_range[0]}, {mtot_range[1]}]"
+
+        for i, line in enumerate(lines):
+            if line.startswith("arguments ="):
+                # Replace --parameter mc with --parameter mtot
+                line = line.replace('--parameter mc', '--parameter mtot')
+                # Use regular expression to find and replace the --mc-range argument
+                line = re.sub(r'--mc-range\s*\'\[.*?\]\'', f"--mtot-range '{mtot_range_str}'", line)
+                lines[i] = line
+
+        with open(input_file, 'w') as file:
+            file.writelines(lines)
+    
+        return
+    
+    def modify_puff_coords(input_file):
+        with open(input_file, 'r') as file:
+            lines = file.readlines()
+        
+        for i, line in enumerate(lines):
+            if line.startswith("arguments ="):
+                # Replace --parameter mc with --parameter mtot
+                line = line.replace('--parameter mc', '--parameter mtot')
+                
+                lines[i] = line
+        with open(input_file, 'w') as file:
+            file.writelines(lines)
+    
+        return
     
     for submit_file in total_file_paths:
-        modify_submit_file(submit_file)
\ No newline at end of file
+        modify_submit_file(submit_file)
+        
+    if opts.alternate_start:
+        # calculate mtotal range
+        mtot_inj = float(config.get('injection-parameters','m1')) + float(config.get('injection-parameters','m2'))
+        mtot_range = [0.5*mtot_inj, 1.5*mtot_inj]
+        for submit_file in CIP_file_paths:
+            modify_CIP_coords(submit_file, mtot_range)
+            
+        for submit_file in PUFF_file_paths:
+            modify_puff_coords(submit_file)
\ No newline at end of file

From 1bd3e0141ae1f77c5de02e6ca72a2b05a6c2b36d Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Fri, 26 Jul 2024 08:11:00 -0700
Subject: [PATCH 16/50] mtot implementation 2

---
 MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py
index 3e6cd6f53..1053253ee 100755
--- a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py
+++ b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py
@@ -689,6 +689,8 @@ def unsafe_parse_arg_string_dict(my_argstr):
     cmd+= " --force-eta-range  " + str(opts.force_eta_range).replace(' ','')
 if not(opts.force_mtot_range is None):
     cmd+= " --force-mtot-range  " + str(opts.force_mtot_range).replace(' ','')
+if opts.use_mtot_coords:
+    cmd+= " --use-mtot-coords "
 if opts.allow_subsolar:
     cmd += " --allow-subsolar "
 if opts.force_chi_max:
@@ -1088,7 +1090,7 @@ def unsafe_parse_arg_string_dict(my_argstr):
             if not(opts.use_mtot_coords):
                 line = line.replace('parameter mc', 'parameter mc --parameter E0 --parameter p_phi0 --use-hyperbolic')
             else:
-                line = line.replace('parameter mtot', 'parameter mtot --parameter E0 --parameter p_phi0 --use-hyperbolic')
+                line = line.replace('parameter mc', 'parameter mtot --parameter E0 --parameter p_phi0 --use-hyperbolic')
         else:
             line = line.replace('parameter-nofit mc', 'parameter-nofit mc --parameter E0 --parameter p_phi0 --use-hyperbolic')
         if not(opts.force_E0_max is None):

From 21ea3d6279a626719c0719b3db0a3660953aa943 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Fri, 26 Jul 2024 11:25:32 -0700
Subject: [PATCH 17/50] mtot implementation 3

---
 MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py | 5 +----
 1 file changed, 1 insertion(+), 4 deletions(-)

diff --git a/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py b/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
index c7fee81af..27fa87fdf 100755
--- a/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
+++ b/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
@@ -1435,10 +1435,7 @@ def lambda_m_estimate(m):
     lnLoffset_late = 15 # default
     helper_cip_args += ' --no-plots --fit-method {}  '.format(fit_method)
     if not opts.internal_use_aligned_phase_coordinates:
-        if not(opts.use_mtot_coords):
-            helper_cip_args += '   --parameter mc --parameter delta_mc '
-        else:
-            helper_cip_args += '   --parameter mtot --parameter delta_mc '
+        helper_cip_args += '   --parameter mc --parameter delta_mc '
     else:
         helper_cip_args += " --parameter-implied mu1 --parameter-implied mu2 --parameter-nofit mc --parameter delta_mc "  
     if 'gp' in fit_method:

From 90849b9ed40020f88c97753f1655d49a5172bd58 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Mon, 12 Aug 2024 06:04:20 -0700
Subject: [PATCH 18/50] mtot coord minor fix

---
 .../Code/test/single_inject/single_inject_RIFT.py              | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py b/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py
index b2df1c1a7..d1080a182 100755
--- a/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py
+++ b/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py
@@ -401,7 +401,6 @@ def add_frames(channel,input_frame, noise_frame,combined_frame):
         pass
     else:
         cmd += ' --first-iteration-jumpstart '
-
     
 os.system(cmd)   
           
@@ -673,4 +672,4 @@ def modify_puff_coords(input_file):
             modify_CIP_coords(submit_file, mtot_range)
             
         for submit_file in PUFF_file_paths:
-            modify_puff_coords(submit_file)
\ No newline at end of file
+            modify_puff_coords(submit_file)

From 1f7c8718e16b2ac4cc7de087be4e7d4b1ded300f Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Tue, 27 Aug 2024 05:42:15 -0700
Subject: [PATCH 19/50] puff + mtot fix

---
 .../Code/bin/helper_LDG_Events.py             |  2 +-
 .../Code/bin/util_RIFT_pseudo_pipe.py         | 11 ++++---
 .../test/single_inject/single_inject_RIFT.py  | 30 +++++++++++--------
 3 files changed, 26 insertions(+), 17 deletions(-)

diff --git a/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py b/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
index 27fa87fdf..6aeb51693 100755
--- a/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
+++ b/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
@@ -1444,7 +1444,7 @@ def lambda_m_estimate(m):
         if not(opts.use_mtot_coords):
             helper_cip_args += mc_range_str_cip + eta_range_str_cip
         else:
-            helper_cip_args += mtot_range_str_cip + eta_range_str_cip
+            helper_cip_args += mtot_range_str_cip + eta_range_str_cip + mc_range_str_cip
     if opts.force_chi_max:
         helper_cip_args += " --chi-max {} ".format(opts.force_chi_max)
     if opts.force_chi_small_max:
diff --git a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py
index 1053253ee..908a8c958 100755
--- a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py
+++ b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py
@@ -1090,7 +1090,9 @@ def unsafe_parse_arg_string_dict(my_argstr):
             if not(opts.use_mtot_coords):
                 line = line.replace('parameter mc', 'parameter mc --parameter E0 --parameter p_phi0 --use-hyperbolic')
             else:
-                line = line.replace('parameter mc', 'parameter mtot --parameter E0 --parameter p_phi0 --use-hyperbolic')
+                #line = line.replace('parameter mc', 'parameter mtot --parameter E0 --parameter p_phi0 --use-hyperbolic')
+                line = line.replace('parameter mc', 'parameter-implied mc --parameter-nofit mtot --parameter-nofit q --parameter E0 --parameter p_phi0 --use-hyperbolic')
+                line = line.replace('parameter delta_mc', 'parameter-implied delta_mc')
         else:
             line = line.replace('parameter-nofit mc', 'parameter-nofit mc --parameter E0 --parameter p_phi0 --use-hyperbolic')
         if not(opts.force_E0_max is None):
@@ -1168,15 +1170,16 @@ def unsafe_parse_arg_string_dict(my_argstr):
 if opts.assume_eccentric:
         puff_params += " --parameter eccentricity --downselect-parameter eccentricity --downselect-parameter-range '[{},{}]' ".format(ecc_min,ecc_max)  
 if opts.assume_hyperbolic:
-        puff_params += " --parameter E0 --downselect-parameter E0 --downselect-parameter-range '[1.0,1.060]' --parameter E0 --downselect-parameter p_phi0 --downselect-parameter-range '[3.8,5.4]' "
+        puff_params += " --parameter E0 "
         if not(opts.force_E0_max is None and opts.force_E0_min is None):
             E0_max = opts.force_E0_max
             E0_min = opts.force_E0_min
-            puff_params.replace("--parameter E0 --downselect-parameter E0 --downselect-parameter-range '[1.0,1.060]'", f"--parameter E0 --downselect-parameter E0 --downselect-parameter-range '[{E0_min},{E0_max}]'")
+            puff_params += f" --downselect-parameter E0 --downselect-parameter-range [{E0_min},{E0_max}] "
+        puff_params += " --parameter p_phi0 "
         if not(opts.force_pphi0_max is None and opts.force_pphi0_min is None):
             pphi0_max = opts.force_pphi0_max
             pphi0_min = opts.force_pphi0_min
-            puff_params.replace("--parameter pphi0 --downselect-parameter pphi0 --downselect-parameter-range '[3.8,5.4]'", f"--parameter pphi0 --downselect-parameter pphi0 --downselect-parameter-range '[{pphi0_min},{pphi0_max}]'")
+            puff_params += f" --downselect-parameter p_phi0 --downselect-parameter-range [{pphi0_min},{pphi0_max}]"
 if opts.assume_highq:
     puff_params = puff_params.replace(' delta_mc ', ' eta ')  # use natural coordinates in the high q strategy. May want to do this always
     puff_max_it +=3
diff --git a/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py b/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py
index d1080a182..547cccf56 100755
--- a/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py
+++ b/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py
@@ -357,12 +357,12 @@ def add_frames(channel,input_frame, noise_frame,combined_frame):
 
 os.chdir(working_dir_full)
 
-
-# write coinc file
-cmd = "util_SimInspiralToCoinc.py --sim-xml mdc.xml.gz --event {}".format(indx)
-for ifo in ifos:
-    cmd += "  --ifo {} ".format(ifo)
-os.system(cmd)
+if not(bypass_frames):
+    # write coinc file
+    cmd = "util_SimInspiralToCoinc.py --sim-xml mdc.xml.gz --event {}".format(indx)
+    for ifo in ifos:
+        cmd += "  --ifo {} ".format(ifo)
+    os.system(cmd)
 
 # Designate rundir 
 rundir = config.get('init','rundir')
@@ -577,15 +577,20 @@ def modify_submit_file(input_file):
         #current_singularity_image = 'rift_o4b_jl-chadhenshaw-teobresums_eccentric-2024-05-14_12-02-56.sif'
         #sif_path = 'osdf:///igwn/cit/staging/james.clark/' + current_singularity_image
         
-        current_singularity_image = 'rift_o4b_jl-2024-07-16_12-21-57.sif'
-        sif_path = 'osdf:///igwn/cit/staging/chad.henshaw/' + current_singularity_image
+        #current_singularity_image = 'rift_o4b_jl-2024-07-16_12-21-57.sif'
+       # sif_path = 'osdf:///igwn/cit/staging/chad.henshaw/' + current_singularity_image
+    
+        current_singularity_image = '/cvmfs/singularity.opensciencegrid.org/james-clark/research-projects-rit/containers-rift_o4b_jl-chadhenshaw-teobresums_eccentric:latest'
         
         with open(input_file, 'r') as file:
             lines = file.readlines()
     
-        lines_to_remove = ['getenv', '+SingularityBindCVMFS', '+SingularityImage']
-        lines_to_add = [f'MY.SingularityImage = "{current_singularity_image}"', 'use_oauth_services = scitokens']
-        modifications = {'request_disk': 'request_disk = 3G', 'requirements': 'requirements = (HAS_SINGULARITY=?=TRUE)&&(IS_GLIDEIN=?=TRUE)', 'transfer_input_files': None}
+        #lines_to_remove = ['getenv', '+SingularityBindCVMFS', '+SingularityImage']
+        lines_to_remove = ['getenv']
+        #lines_to_add = [f'MY.SingularityImage = "{current_singularity_image}"', 'use_oauth_services = scitokens']
+        lines_to_add = []
+        #modifications = {'request_disk': 'request_disk = 3G', 'requirements': 'requirements = (HAS_SINGULARITY=?=TRUE)&&(IS_GLIDEIN=?=TRUE)', 'transfer_input_files': None}
+        modifications = {'request_disk': 'request_disk = 50M', 'requirements': 'requirements = (HAS_SINGULARITY=?=TRUE)&&(IS_GLIDEIN=?=TRUE)&&(TARGET.Machine =!= "deepclean.ldas.cit")'}
 
         modified_lines = []
         queue_line = None
@@ -602,7 +607,8 @@ def modify_submit_file(input_file):
             for key, new_value in modifications.items():
                 if stripped_line.startswith(key):
                     if key == 'transfer_input_files':
-                        line = stripped_line + ',' + sif_path + '\n'
+                        #line = stripped_line + ',' + sif_path + '\n'
+                        print('skipping this for cvfms')
                     else:
                         line = new_value + '\n'
                     modified_lines.append(line)

From c752f6872f955a4c0dd39e73b7551e2d2134798c Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Wed, 4 Sep 2024 08:42:01 -0700
Subject: [PATCH 20/50] tapering change in hoft, NOT hlmoft

---
 MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
index e4fe77889..19400dba9 100644
--- a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
+++ b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
@@ -2885,7 +2885,7 @@ def hoft(P, Fp=None, Fc=None,**kwargs):
                 'q'                  : M1/M2,
                 'H_hyp'              : P.E0, # energy at initial separation
                 'j_hyp'              : P.p_phi0, # angular momentum at initial separation
-                'r_hyp'              : 6000.0,
+                'r_hyp'              : 12000.0,
                 'LambdaAl2'            : P.lambda1,
                 'LambdaBl2'            : P.lambda2,
                 'chi1'              : P.s1z, #note that there are no transverse spins
@@ -2902,7 +2902,7 @@ def hoft(P, Fp=None, Fc=None,**kwargs):
                 'use_geometric_units': 0,
                 'interp_uniform_grid': 1,
                 'initial_frequency'  : P.fmin,
-                'ode_tmax'           : 3e4,
+                'ode_tmax'           : 4.05e4,
                 'distance'           : P.dist/(lal.PC_SI*1e6),
                 'inclination'        : P.incl,
                 'output_hpc'         : 0 # output plus and cross polarizations, 0=no
@@ -2946,7 +2946,7 @@ def hoft(P, Fp=None, Fc=None,**kwargs):
             amp_max_ind = np.argmax(amp)        
             amp_norm = amp / amp[amp_max_ind] # normalize amplitude for peak finding                
             # peak finding to determine system type
-            peaks, props = signal.find_peaks(amp_norm, height = 0.25)    
+            peaks, props = signal.find_peaks(amp_norm, height = 0.25, prominence = 0.1)    
             peak_heights = props['peak_heights']
             # filtering out peaks so we only keep the local maxima
             indices_to_keep = set()
@@ -3048,7 +3048,7 @@ def hoft(P, Fp=None, Fc=None,**kwargs):
                 if count == 0:
                     continue
                 if np.abs(value-ht.data.data[0]) > 0.01 * np.abs(ht.data.data[0]):
-                    n_samp=int(count/2)
+                    n_samp=int(count/16)
                     break
             vectaper= 0.5 + 0.5*np.cos(np.pi* (1-np.arange(n_samp)/(1.*n_samp)))
             nmax = np.argmax(ht.data.data)
@@ -3057,11 +3057,13 @@ def hoft(P, Fp=None, Fc=None,**kwargs):
             # If Capture waveform, no end taper
             
             # peak finding to determine system type
-            peaks, props = signal.find_peaks(ht.data.data, height = 0.25*np.abs(np.amax(ht.data.data)))
+            height_thresh = 0.25*np.abs(np.amax(ht.data.data))
+            prom_thresh = 0.1*np.abs(np.amax(ht.data.data))
+            peaks, props = signal.find_peaks(ht.data.data, height = height_thresh, prominence = prom_thresh)
             peak_heights = props['peak_heights']
             indices_to_keep = set()
             sorted_indices = np.argsort(peak_heights)[::-1]
-            tol = 300 # hardcoded peak spacing tolerance. if spacing is less than tol, discard the peak
+            tol = int(pars['srate_interp'] / 13.65) # 300 samples at srate of 4096 - minimum distance between peaks.
             for i in sorted_indices:
                 peak = peaks[i]                
                 keep = True

From a4eb5502dbe7473cb4fbed1562e0f22cda9372af Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Wed, 11 Sep 2024 05:57:40 -0700
Subject: [PATCH 21/50] hoft epoch peak fix, no hlmoft

---
 MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
index 19400dba9..ac2555de6 100644
--- a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
+++ b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
@@ -2946,7 +2946,9 @@ def hoft(P, Fp=None, Fc=None,**kwargs):
             amp_max_ind = np.argmax(amp)        
             amp_norm = amp / amp[amp_max_ind] # normalize amplitude for peak finding                
             # peak finding to determine system type
-            peaks, props = signal.find_peaks(amp_norm, height = 0.25, prominence = 0.1)    
+            height_thresh = 0.25*np.abs(amp_norm)
+            prom_thresh = 0.1*np.abs(amp_norm)
+            peaks, props = signal.find_peaks(amp_norm, height = height_thresh, prominence = prom_thresh)    
             peak_heights = props['peak_heights']
             # filtering out peaks so we only keep the local maxima
             indices_to_keep = set()

From df163bd32f2415efbfe005e1a7a2e4753b74e297 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Wed, 25 Sep 2024 08:37:32 -0700
Subject: [PATCH 22/50] reverting hoft change in lalsimutils, adding hlmoft to
 WriteFrame

---
 MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py       | 9 +++++----
 MonteCarloMarginalizeCode/Code/bin/util_LALWriteFrame.py | 7 ++++++-
 2 files changed, 11 insertions(+), 5 deletions(-)

diff --git a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
index ac2555de6..1389a186c 100644
--- a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
+++ b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
@@ -2885,7 +2885,7 @@ def hoft(P, Fp=None, Fc=None,**kwargs):
                 'q'                  : M1/M2,
                 'H_hyp'              : P.E0, # energy at initial separation
                 'j_hyp'              : P.p_phi0, # angular momentum at initial separation
-                'r_hyp'              : 12000.0,
+                'r_hyp'              : 6000.0,
                 'LambdaAl2'            : P.lambda1,
                 'LambdaBl2'            : P.lambda2,
                 'chi1'              : P.s1z, #note that there are no transverse spins
@@ -2902,7 +2902,7 @@ def hoft(P, Fp=None, Fc=None,**kwargs):
                 'use_geometric_units': 0,
                 'interp_uniform_grid': 1,
                 'initial_frequency'  : P.fmin,
-                'ode_tmax'           : 4.05e4,
+                'ode_tmax'           : 3e4,
                 'distance'           : P.dist/(lal.PC_SI*1e6),
                 'inclination'        : P.incl,
                 'output_hpc'         : 0 # output plus and cross polarizations, 0=no
@@ -3050,7 +3050,7 @@ def hoft(P, Fp=None, Fc=None,**kwargs):
                 if count == 0:
                     continue
                 if np.abs(value-ht.data.data[0]) > 0.01 * np.abs(ht.data.data[0]):
-                    n_samp=int(count/16)
+                    n_samp=int(count/2)
                     break
             vectaper= 0.5 + 0.5*np.cos(np.pi* (1-np.arange(n_samp)/(1.*n_samp)))
             nmax = np.argmax(ht.data.data)
@@ -3697,7 +3697,8 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil
                     modes_used_new.append((4,0))
                     hlmtmp2[(4,0)]=np.array(hlmtmp[k])
             modes_used=modes_used_new
-            print(modes_used,hlmtmp,hlmtmp2)
+            if not(hyp_wav):
+                print(modes_used,hlmtmp,hlmtmp2)
 #        for count,mode in enumerate(modes_used):
 #            hlmtmp2[mode]=np.array(hlmtmp[str(count)])
         for mode in modes_used:
diff --git a/MonteCarloMarginalizeCode/Code/bin/util_LALWriteFrame.py b/MonteCarloMarginalizeCode/Code/bin/util_LALWriteFrame.py
index 66bae5ed5..e1158f4c2 100755
--- a/MonteCarloMarginalizeCode/Code/bin/util_LALWriteFrame.py
+++ b/MonteCarloMarginalizeCode/Code/bin/util_LALWriteFrame.py
@@ -40,6 +40,7 @@
 parser.add_argument("--mass2",default=1.4,type=float,help='Mass 2 (solar masses)')
 parser.add_argument("--verbose", action="store_true",default=False)
 parser.add_argument("--l-max",default=4,type=float,help='Lmax number of modes')
+parser.add_argument('--gen-hlmoft', action='store_true', help='Creates hoft from hlmoft')
 opts=  parser.parse_args()
 
 
@@ -86,7 +87,11 @@
 
 
 # Generate signal
-hoft = lalsimutils.hoft(P,Lmax=opts.l_max)   # include translation of source, but NOT interpolation onto regular time grid
+if opts.gen_hlmoft:
+    hlmT = lalsimutils.hlmoft(P, Lmax=opts.l_max)
+    hoft = lalsimutils.hoft_from_hlm(hlmT, P, return_complex=False)
+else:
+    hoft = lalsimutils.hoft(P,Lmax=opts.l_max)   # include translation of source, but NOT interpolation onto regular time grid
 # zero pad to be opts.seglen long, if necessary
 if opts.seglen/hoft.deltaT > hoft.data.length:
     TDlenGoal = int(opts.seglen/hoft.deltaT)

From 06cce41978e867e4227cf5367dfe5c553d9b0f0d Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Wed, 25 Sep 2024 08:42:33 -0700
Subject: [PATCH 23/50] one more hyp_wav instance

---
 MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
index 1389a186c..bf0971b38 100644
--- a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
+++ b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
@@ -2852,6 +2852,7 @@ def hoft(P, Fp=None, Fc=None,**kwargs):
         M1=P.m1/lal.MSUN_SI
         M2=P.m2/lal.MSUN_SI
         nu=M1*M2/((M1+M2)**2)
+        hyp_wav = False
         print(P.eccentricity, P.E0)
         if (P.eccentricity == 0.0 and P.E0 == 0.0):
             print("Using ResumS master; not eccentric")

From 6b2bc52fa9e56321461198668b14a0b446b16636 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Wed, 6 Nov 2024 05:55:20 -0800
Subject: [PATCH 24/50] adjusted peak finder to deal with non-interacting cases

---
 MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
index bf0971b38..404f668e9 100644
--- a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
+++ b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
@@ -2969,6 +2969,11 @@ def hoft(P, Fp=None, Fc=None,**kwargs):
             if len(filtered_peaks) == 1:
                 # scatter case OR plunge case, we can set the epoch normally
                 hpepoch = -P.deltaT*np.argmax(np.abs(hptmp)**2+np.abs(hctmp)**2)
+            elif len(filtered_peaks) == 0:
+                # meaningless waveform, essentially a non-interacting case. Set the epoch normally
+                # These points should return very low likelihood and not interfere with the analysis
+                print('WARNING: no peak detected; non-interacting hyperbolic case')
+                hpepoch = -P.deltaT*np.argmax(np.abs(hptmp)**2+np.abs(hctmp)**2)
             else:
                 # capture case, we need to force the epoch to be the last peak
                 hpepoch = -P.deltaT*filtered_peaks[-1]
@@ -3623,6 +3628,11 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil
             if len(filtered_peaks) == 1:
                 # scatter case OR plunge case, we can set the epoch normally
                 hpepoch = -P.deltaT*np.argmax(np.abs(hptmp)**2+np.abs(hctmp)**2)
+            elif len(filtered_peaks) == 0:
+                # meaningless waveform, essentially a non-interacting case. Set the epoch normally
+                # These points should return very low likelihood and not interfere with the analysis
+                print('WARNING: no peak detected; non-interacting hyperbolic case')
+                hpepoch = -P.deltaT*np.argmax(np.abs(hptmp)**2+np.abs(hctmp)**2)                
             else:
                 # capture case, we need to force the epoch to be the last peak
                 hpepoch = -P.deltaT*filtered_peaks[-1]

From a18b8b85b318882b1c5fef4042f973d7159ce1bd Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Fri, 22 Nov 2024 09:38:09 -0800
Subject: [PATCH 25/50] add prominence to hlmoft epoch find

---
 MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
index 404f668e9..a0a021121 100644
--- a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
+++ b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
@@ -3608,7 +3608,9 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil
             amp_max_ind = np.argmax(amp)        
             amp_norm = amp / amp[amp_max_ind] # normalize amplitude for peak finding                
             # peak finding to determine system type
-            peaks, props = signal.find_peaks(amp_norm, height = 0.25)    
+            height_thresh = 0.25*np.abs(amp_norm)
+            prom_thresh = 0.1*np.abs(amp_norm)
+            peaks, props = signal.find_peaks(amp_norm, height = height_thresh, prominence = prom_thresh)    
             peak_heights = props['peak_heights']
             # filtering out peaks so we only keep the local maxima
             indices_to_keep = set()

From 09e03b91c7ac0dd80b31b1d707647baca9f2ba65 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Fri, 10 Jan 2025 08:34:14 -0800
Subject: [PATCH 26/50] tapering order change and Lmax to kwaargs

---
 .../Code/RIFT/lalsimutils.py                  | 126 ++++++++++++------
 1 file changed, 85 insertions(+), 41 deletions(-)

diff --git a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
index a0a021121..47b99c82e 100644
--- a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
+++ b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
@@ -1853,14 +1853,16 @@ def copy_lsctables_sim_inspiral(self, row):
         # Call the function to read lalmetaio.SimInspiral format
         self.copy_sim_inspiral(swigrow)
 
-    def scale_to_snr(self,new_SNR,psd, ifo_list,analyticPSD_Q=True):
+    def scale_to_snr(self,new_SNR,psd, ifo_list,analyticPSD_Q=True, **kwargs):
         """
         scale_to_snr
           - evaluates network SNR in the ifo list provided (assuming *constant* psd for all..may change)
           - uses network SNR to rescale the distance of the source, so the SNR is now  new_SNR
           - returns current_SNR, for sanity
         """
-        deltaF=findDeltaF(self)
+        Lmax = kwargs.get('Lmax', 4)  # Default to 4 if not specified in kwargs
+        
+        deltaF=findDeltaF(self, Lmax=Lmax)
         det_orig = self.detector
         IP = Overlap(fLow=self.fmin, fNyq=1./self.deltaT/2., deltaF=deltaF, psd=psd, full_output=True,analyticPSD_Q=analyticPSD_Q)
 
@@ -1869,7 +1871,7 @@ def scale_to_snr(self,new_SNR,psd, ifo_list,analyticPSD_Q=True):
         for det in ifo_list:
             self.detector = det
             self.radec = True
-            h=hoff(self)
+            h=hoff(self, Lmax=Lmax)
             rho_ifo[det] = IP.norm(h)
             current_SNR_squared +=rho_ifo[det]*rho_ifo[det]
         current_SNR = np.sqrt(current_SNR_squared)
@@ -2750,7 +2752,7 @@ def nextPow2(length):
     """
     return int(2**np.ceil(np.log2(length)))
 
-def findDeltaF(P):
+def findDeltaF(P, **kwargs):
     """
     Given ChooseWaveformParams P, generate the TD waveform,
     round the length to the next power of 2,
@@ -2758,7 +2760,9 @@ def findDeltaF(P):
     This is useful b/c deltaF is needed to define an inner product
     which is needed for norm_hoft and norm_hoff functions
     """
-    h = hoft(P)
+    Lmax = kwargs.get('Lmax', 4)  # Default to 4 if not specified in kwargs
+    
+    h = hoft(P, Lmax=Lmax)
     return 1./(nextPow2(h.data.length) * P.deltaT)
 
 def estimateWaveformDuration(P,LmaxEff=2):
@@ -2840,7 +2844,7 @@ def hoft(P, Fp=None, Fc=None,**kwargs):
         extra_waveform_args.update(kwargs['extra_waveform_args'])
     extra_params = P.to_lal_dict_extended(extra_args_dict=extra_waveform_args)
     if P.approx==lalsim.TEOBResumS and has_external_teobresum and info_use_ext:
-        Lmax=4
+        Lmax = kwargs.get('Lmax', 4)  # Default to 4 if not specified in kwargs
         modes_used = []
         distance_s = P.dist/lal.C_SI
         m_total_s = MsunInSec*(P.m1+P.m2)/lal.MSUN_SI
@@ -2942,10 +2946,8 @@ def hoft(P, Fp=None, Fc=None,**kwargs):
         else:
             ## custom epoch for the hyperbolic case ##            
             # wf amplitude
-            amp = np.abs(hptmp)**2+np.abs(hctmp)**2
-            amp_times = P.deltaT*np.arange(len(amp))
-            amp_max_ind = np.argmax(amp)        
-            amp_norm = amp / amp[amp_max_ind] # normalize amplitude for peak finding                
+            amp = np.sqrt(np.abs(hptmp)**2+np.abs(hctmp)**2)
+            amp_norm = amp / np.amax(amp) # normalize amplitude for peak finding                
             # peak finding to determine system type
             height_thresh = 0.25*np.abs(amp_norm)
             prom_thresh = 0.1*np.abs(amp_norm)
@@ -3058,16 +3060,14 @@ def hoft(P, Fp=None, Fc=None,**kwargs):
                 if np.abs(value-ht.data.data[0]) > 0.01 * np.abs(ht.data.data[0]):
                     n_samp=int(count/2)
                     break
-            vectaper= 0.5 + 0.5*np.cos(np.pi* (1-np.arange(n_samp)/(1.*n_samp)))
-            nmax = np.argmax(ht.data.data)
-            ht.data.data[0:n_samp] *= vectaper
-            # If Scattering waveform, tapers same amount as early taper
-            # If Capture waveform, no end taper
             
             # peak finding to determine system type
-            height_thresh = 0.25*np.abs(np.amax(ht.data.data))
-            prom_thresh = 0.1*np.abs(np.amax(ht.data.data))
-            peaks, props = signal.find_peaks(ht.data.data, height = height_thresh, prominence = prom_thresh)
+            ht_norm = ht.data.data/np.amax(ht.data.data)
+            height_thresh = 0.25
+            prom_thresh = 0.1
+            
+            peaks, props = signal.find_peaks(ht_norm, height = height_thresh, prominence = prom_thresh)
+            
             peak_heights = props['peak_heights']
             indices_to_keep = set()
             sorted_indices = np.argsort(peak_heights)[::-1]
@@ -3083,12 +3083,33 @@ def hoft(P, Fp=None, Fc=None,**kwargs):
                     indices_to_keep.add(i)                    
             filtered_peaks = peaks[list(indices_to_keep)]
             
+            vectaper= 0.5 + 0.5*np.cos(np.pi* (1-np.arange(n_samp)/(1.*n_samp))) # this tapers the start
+            
+            nmax = np.argmax(ht.data.data)
+            ht.data.data[0:n_samp] *= vectaper
+            
             if len(filtered_peaks) == 1:
                 # check if a scatter or a plunge            
                 if np.abs(ht.data.length-nmax) > 3e3:
+                    #scatter
+                    print('Tapering for scatter waveform')
                     n_samp2=n_samp
-                    vectaper2= 0.5 - 0.5*np.cos(np.pi* (1-np.arange(n_samp2+1)/(1.*n_samp2)))
+                    vectaper2 = 0.5 + 0.5 * np.cos(np.pi * np.arange(n_samp2 + 1) / (1. * n_samp2)) # end taper
                     ht.data.data[-(n_samp2+1):] *= vectaper2
+                else:
+                    # plunge, only need to taper the start
+                    print('Plunge waveform, only start taper')
+                    
+            elif len(filtered_peaks) == 0:
+                print('Non-interactive hyperbolic waveform, tapering both ends')
+                n_samp2=n_samp
+                vectaper2 = 0.5 + 0.5 * np.cos(np.pi * np.arange(n_samp2 + 1) / (1. * n_samp2))
+                ht.data.data[-(n_samp2+1):] *= vectaper2
+            
+            else:
+                # zoom-whirl case, taper just the start
+                print('Zoom-whirl waveform, only start taper')
+                 
                 
     if P.deltaF is not None:
         TDlen = int(1./P.deltaF * 1./P.deltaT)
@@ -3096,7 +3117,7 @@ def hoft(P, Fp=None, Fc=None,**kwargs):
         ht = lal.ResizeREAL8TimeSeries(ht, 0, TDlen)
     return ht
 
-def hoff(P, Fp=None, Fc=None, fwdplan=None):
+def hoff(P, Fp=None, Fc=None, fwdplan=None, **kwargs):
     """
     Generate a FD waveform from ChooseWaveformParams P.
     Will return a COMPLEX16FrequencySeries object.
@@ -3112,6 +3133,8 @@ def hoff(P, Fp=None, Fc=None, fwdplan=None):
         If P.deltaF == None, the TD waveform will be zero-padded
         to the next power of 2.
     """
+    Lmax = kwargs.get('Lmax', 4)  # Default to 4 if not specified in kwargs
+    
     # For FD approximants, use the ChooseFDWaveform path = hoff_FD
     if lalsim.SimInspiralImplementedFDApproximants(P.approx)==1:
         # Raise exception if unused arguments were specified
@@ -3121,11 +3144,11 @@ def hoff(P, Fp=None, Fc=None, fwdplan=None):
 
     # For TD approximants, do ChooseTDWaveform + FFT path = hoff_TD
     else:
-        hf = hoff_TD(P, Fp, Fc, fwdplan)
+        hf = hoff_TD(P, Fp, Fc, fwdplan, Lmax=Lmax)
 
     return hf
 
-def hoff_TD(P, Fp=None, Fc=None, fwdplan=None):
+def hoff_TD(P, Fp=None, Fc=None, fwdplan=None, **kwargs):
     """
     Generate a FD waveform from ChooseWaveformParams P
     by creating a TD waveform, zero-padding and
@@ -3145,7 +3168,9 @@ def hoff_TD(P, Fp=None, Fc=None, fwdplan=None):
 
     Returns a COMPLEX16FrequencySeries object
     """
-    ht = hoft(P, Fp, Fc)
+    Lmax = kwargs.get('Lmax', 4)  # Default to 4 if not specified in kwargs
+    
+    ht = hoft(P, Fp, Fc, Lmax=Lmax)
 
     if P.deltaF == None: # h(t) was not zero-padded, so do it now
         TDlen = nextPow2(ht.data.length)
@@ -3603,14 +3628,13 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil
         else:
             ## custom epoch for the hyperbolic case ##            
             # wf amplitude
-            amp = np.abs(hptmp)**2+np.abs(hctmp)**2
-            amp_times = P.deltaT*np.arange(len(amp))
-            amp_max_ind = np.argmax(amp)        
-            amp_norm = amp / amp[amp_max_ind] # normalize amplitude for peak finding                
+            amp = np.sqrt(np.abs(hptmp)**2+np.abs(hctmp)**2)
+            amp_norm = amp / np.amax(amp) # normalize amplitude for peak finding
+            
             # peak finding to determine system type
-            height_thresh = 0.25*np.abs(amp_norm)
-            prom_thresh = 0.1*np.abs(amp_norm)
-            peaks, props = signal.find_peaks(amp_norm, height = height_thresh, prominence = prom_thresh)    
+            height_thresh = 0.25
+            prom_thresh = 0.1
+            peaks, props = signal.find_peaks(amp_norm, height = height_thresh, prominence = prom_thresh) 
             peak_heights = props['peak_heights']
             # filtering out peaks so we only keep the local maxima
             indices_to_keep = set()
@@ -3758,18 +3782,17 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil
                     if np.abs(np.real(value)-np.real(hlm[(2,2)].data.data[0])) > 0.01 * np.abs(np.real(hlm[(2,2)].data.data[0])):
                         n_samp=int(count/2)
                         break
-            vectaper= 0.5 + 0.5*np.cos(np.pi* (1-np.arange(n_samp)/(1.*n_samp)))
-            nmax = np.argmax(hlm[(2,2)].data.data)
-            for mode in modes_used_new2:
-                hlm[mode].data.data[0:n_samp] *= vectaper
-            # If Capture waveform, no end taper
-            # If Scattering waveform, tapers end same amount as early taper
-            # peak finding to determine system type
-            peaks, props = signal.find_peaks(hlm[(2,2)].data.data, height = 0.25*np.amax(hlm[(2,2)].data.data))
+                        
+            # peak finding to determine system type based on 22 mode
+            wf_22 = hlm[(2,2)].data.data
+            wf_22_norm = wf_22 / np.amax(wf_22)
+            height_thresh = 0.25
+            prom_thresh = 0.1
+            peaks, props = signal.find_peaks(wf_22_norm, height = height_thresh, prominence = prom_thresh)
             peak_heights = props['peak_heights']
             indices_to_keep = set()
             sorted_indices = np.argsort(peak_heights)[::-1]
-            tol = 300 # hardcoded tolerance
+            tol = int(pars['srate_interp'] / 13.65) # 300 samples at srate of 4096 - minimum distance between peaks.
             for i in sorted_indices:
                 peak = peaks[i]                
                 keep = True
@@ -3781,13 +3804,34 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil
                     indices_to_keep.add(i)                    
             filtered_peaks = peaks[list(indices_to_keep)]
             
+            # Always taper the start
+            
+            vectaper= 0.5 + 0.5*np.cos(np.pi* (1-np.arange(n_samp)/(1.*n_samp)))
+            nmax = np.argmax(hlm[(2,2)].data.data)
+            for mode in modes_used_new2:
+                hlm[mode].data.data[0:n_samp] *= vectaper
+            
             if len(filtered_peaks) == 1:
-                #check for scatter waveform            
+                #check if scatter or plunge            
                 if np.abs(hlm[(2,2)].data.length-nmax) > 3e3:
+                    print('Scatter waveform, tapering both ends')
                     n_samp2=n_samp
-                    vectaper2= 0.5 - 0.5*np.cos(np.pi* (1-np.arange(n_samp2+1)/(1.*n_samp2)))
+                    vectaper2= 0.5 + 0.5 * np.cos(np.pi * np.arange(n_samp2 + 1) / (1. * n_samp2))
                     for mode in modes_used_new2:
                         hlm[mode].data.data[-(n_samp2+1):] *= vectaper2
+                else:
+                    print('Plunge waveform, only start taper')
+                    
+            elif len(filtered_peaks) ==0:
+                print('Non-interactive hyperbolic waveform, tapering both ends')
+                n_samp2=n_samp
+                vectaper2 = 0.5 + 0.5 * np.cos(np.pi * np.arange(n_samp2 + 1) / (1. * n_samp2))
+                for mode in modes_used_new2:
+                        hlm[mode].data.data[-(n_samp2+1):] *= vectaper2
+                        
+            else:
+                print('Zoom-whirl waveform, only start taper')                        
+            
         for mode in modes_used_new2:
             if not (P.deltaF is None):
                 TDlen = int(1./P.deltaF * 1./P.deltaT)

From 1c48bbe1025a42c1287a916d3d009e392cb0b7e5 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Fri, 10 Jan 2025 08:40:22 -0800
Subject: [PATCH 27/50] hypclass initial commit

---
 .../Code/RIFT/lalsimutils.py                  | 101 +++++++++++++++++-
 1 file changed, 100 insertions(+), 1 deletion(-)

diff --git a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
index 47b99c82e..ad9c92d2d 100644
--- a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
+++ b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
@@ -332,7 +332,7 @@ def lsu_StringFromPNOrder(order):
 #
 # Class to hold arguments of ChooseWaveform functions
 #
-valid_params = ['m1', 'm2', 's1x', 's1y', 's1z', 's2x', 's2y', 's2z', 'chi1_perp', 'chi2_perp', 'chi1_perp_bar', 'chi2_perp_bar','chi1_perp_u', 'chi2_perp_u', 's1z_bar', 's2z_bar', 'lambda1', 'lambda2', 'theta','phi', 'phiref',  'psi', 'incl', 'tref', 'dist', 'mc', 'mc_ecc', 'eta', 'delta_mc', 'chi1', 'chi2', 'thetaJN', 'phiJL', 'theta1', 'theta2', 'cos_theta1', 'cos_theta2',  'theta1_Jfix', 'theta2_Jfix', 'psiJ', 'beta', 'cos_beta', 'sin_phiJL', 'cos_phiJL', 'phi12', 'phi1', 'phi2', 'LambdaTilde', 'DeltaLambdaTilde', 'lambda_plus', 'lambda_minus', 'q', 'mtot','xi','chiz_plus', 'chiz_minus', 'chieff_aligned','fmin','fref', "SOverM2_perp", "SOverM2_L", "DeltaOverM2_perp", "DeltaOverM2_L", "shu","ampO", "phaseO",'eccentricity','chi_pavg','mu1','mu2','eos_table_index', 'E0', 'p_phi0']
+valid_params = ['m1', 'm2', 's1x', 's1y', 's1z', 's2x', 's2y', 's2z', 'chi1_perp', 'chi2_perp', 'chi1_perp_bar', 'chi2_perp_bar','chi1_perp_u', 'chi2_perp_u', 's1z_bar', 's2z_bar', 'lambda1', 'lambda2', 'theta','phi', 'phiref',  'psi', 'incl', 'tref', 'dist', 'mc', 'mc_ecc', 'eta', 'delta_mc', 'chi1', 'chi2', 'thetaJN', 'phiJL', 'theta1', 'theta2', 'cos_theta1', 'cos_theta2',  'theta1_Jfix', 'theta2_Jfix', 'psiJ', 'beta', 'cos_beta', 'sin_phiJL', 'cos_phiJL', 'phi12', 'phi1', 'phi2', 'LambdaTilde', 'DeltaLambdaTilde', 'lambda_plus', 'lambda_minus', 'q', 'mtot','xi','chiz_plus', 'chiz_minus', 'chieff_aligned','fmin','fref', "SOverM2_perp", "SOverM2_L", "DeltaOverM2_perp", "DeltaOverM2_L", "shu","ampO", "phaseO",'eccentricity','chi_pavg','mu1','mu2','eos_table_index', 'E0', 'p_phi0', 'hypclass']
 
 tex_dictionary  = {
  "mtot": '$M$',
@@ -886,6 +886,105 @@ def extract_param(self,p):
             return (self.m2+self.m1)
         if p == 'q':
             return self.m2/self.m1
+        
+        ####################################
+        # EXPERIMENTING WITH HYPCLASS HERE #
+        
+        if p == 'hypclass':
+            
+            
+            # check if valid
+            if self.E0 == 0.0:
+                print('Invalid use of hypclass: non-hyperbolic configuration')
+                return None
+            
+            # run the classifier
+            pars = {
+                'M'                  : (self.m1+self.m2)/lal.MSUN_SI,
+                'q'                  : self.m1/self.m2,
+                'H_hyp'              : self.E0, # energy at initial separation
+                'j_hyp'              : self.p_phi0, # angular momentum at initial separation
+                'r_hyp'              : 6000.0,
+                'LambdaAl2'            : self.lambda1,
+                'LambdaBl2'            : self.lambda2,
+                'chi1'              : self.s1z, #note that there are no transverse spins
+                'chi2'              : self.s2z,
+                'dt'                : self.deltaT,
+                'domain'             : 0, # 0 sets time domain
+                'arg_out'            : 1, # Request multipoles and dynamics as output of the function call - 1=yes
+                'nqc'                : 2, # sets the NQCs, 2=no
+                'nqc_coefs_hlm'      : 0, # Option for the NQC model used in the waveform. 0=none
+                'nqc_coefs_flx'      : 0, # Option for the NQC model used in the flux. 0=none
+                'use_mode_lm'        : [1,-1], # 2\pm 2 modes
+                'output_lm'          : [1,-1],
+                'srate_interp'       : 1./self.deltaT,
+                'use_geometric_units': 0,
+                'interp_uniform_grid': 1,
+                'initial_frequency'  : self.fmin,
+                'ode_tmax'           : 3e4,
+                'distance'           : self.dist/(lal.PC_SI*1e6),
+                'inclination'        : self.incl,
+                'output_hpc'         : 0 # output plus and cross polarizations, 0=no
+            }
+            
+            #print('Classifying hyperbolic waveform...')
+            #print(pars)
+            
+            t, hptmp, hctmp, hlmtmp, dym = EOBRun_module.EOBRunPy(pars)
+            
+            # peak finding to classifiy
+            
+            # wf amplitude
+            amp = np.sqrt(np.abs(hptmp)**2+np.abs(hctmp)**2)
+            amp_norm = amp / np.amax(amp) # normalize amplitude for peak finding                
+            # peak finding to determine system type
+            height_thresh = 0.25
+            prom_thresh = 0.1
+            peaks, props = signal.find_peaks(amp_norm, height = height_thresh, prominence = prom_thresh)    
+            peak_heights = props['peak_heights']
+            # filtering out peaks so we only keep the local maxima
+            indices_to_keep = set()
+            sorted_indices = np.argsort(peak_heights)[::-1]
+            tol = int(pars['srate_interp'] / 13.65) # 300 samples at srate of 4096 - minimum distance between peaks.
+            for i in sorted_indices:
+                peak = peaks[i]                
+                keep = True
+                for kept_index in indices_to_keep:
+                    if abs(peaks[kept_index] - peak) <= tol:
+                        keep = False
+                        break                
+                if keep:
+                    indices_to_keep.add(i)                    
+            filtered_peaks = peaks[list(indices_to_keep)]
+            
+            # parsing number of peaks after filtering against distance tolerance
+            if len(filtered_peaks) == 1:
+                # scatter case OR plunge case
+                
+                if np.abs(len(hptmp) - np.argmax(hptmp)) > pars['srate_interp']/5.46133: #3000 samples at 16384 [Hz]
+                    # scatter waveform
+                    print('Scatter')
+                    return 'scatter'
+                
+                else:
+                    # plunge waveform
+                    print('Plunge')
+                    return 'plunge'
+                
+            elif len(filtered_peaks) == 0:
+                # meaningless waveform
+                print('MEANINGLESS')
+                return None
+            
+            else:
+                # zoom whirl waveform
+                print('Zoom Whirl')
+                return 'zoomwhirl'                
+        
+        
+        ####################################
+        
+        
         if p == 'delta' or p=='delta_mc':  # Same access routine
             return (self.m1-self.m2)/(self.m1+self.m2)
         if p == 'mc':

From 9df5344ecce26dcb6551b61b7e17698fadb759ec Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Fri, 10 Jan 2025 08:45:05 -0800
Subject: [PATCH 28/50] force-scatter option in overlapgrid

---
 .../Code/bin/util_ManualOverlapGrid.py                 | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/MonteCarloMarginalizeCode/Code/bin/util_ManualOverlapGrid.py b/MonteCarloMarginalizeCode/Code/bin/util_ManualOverlapGrid.py
index 42e184b66..6a2759ab6 100755
--- a/MonteCarloMarginalizeCode/Code/bin/util_ManualOverlapGrid.py
+++ b/MonteCarloMarginalizeCode/Code/bin/util_ManualOverlapGrid.py
@@ -245,6 +245,7 @@
 parser.add_argument("--verbose", action="store_true",default=False, help="Extra warnings")
 parser.add_argument("--extra-verbose", action="store_true",default=False, help="Lots of messages")
 parser.add_argument("--save-plots",default=False,action='store_true', help="Write plots to file (only useful for OSX, where interactive is default")
+parser.add_argument('--force-scatter',default=False,action='store_true',help='For hyperbolic analyses forces only scatter grid points.')
 opts=  parser.parse_args()
 if opts.inj_file_out:
     opts.fname = opts.inj_file_out.replace(".xml.gz","")
@@ -366,6 +367,15 @@ def evaluate_overlap_on_grid(hfbase,param_names, grid):
         for param in downselect_dict:
             if Pgrid.extract_param(param) < downselect_dict[param][0] or Pgrid.extract_param(param) > downselect_dict[param][1]:
                 include_item =False
+                
+        if opts.force_scatter:
+            # removes non-scatter points from the hyperbolic grid
+            hypclass = Pgrid.extract_param('hypclass')
+            if hypclass == 'scatter':
+                include_item = True
+            else:
+                include_item = False
+                
         if include_item:
          grid_revised.append(line)
          if Pgrid.m2 <= Pgrid.m1:  # do not add grid elements with m2> m1, to avoid possible code pathologies !

From 8ad3d68c671fdc992f6bf97c54e45e564f6ef92a Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Mon, 20 Jan 2025 09:00:44 -0800
Subject: [PATCH 29/50] removing print statements

---
 MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
index ad9c92d2d..3fcd5f8f8 100644
--- a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
+++ b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
@@ -963,22 +963,18 @@ def extract_param(self,p):
                 
                 if np.abs(len(hptmp) - np.argmax(hptmp)) > pars['srate_interp']/5.46133: #3000 samples at 16384 [Hz]
                     # scatter waveform
-                    print('Scatter')
                     return 'scatter'
                 
                 else:
                     # plunge waveform
-                    print('Plunge')
                     return 'plunge'
                 
             elif len(filtered_peaks) == 0:
                 # meaningless waveform
-                print('MEANINGLESS')
                 return None
             
             else:
                 # zoom whirl waveform
-                print('Zoom Whirl')
                 return 'zoomwhirl'                
         
         

From 5efd90edc9add62477113470cd1ded58e4b3a536 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Mon, 20 Jan 2025 09:01:31 -0800
Subject: [PATCH 30/50] add --force-scatter-grids option

---
 MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py b/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
index 6aeb51693..f2da3e33f 100755
--- a/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
+++ b/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
@@ -252,8 +252,13 @@ def get_observing_run(t):
 parser.add_argument("--use-cvmfs-frames",action='store_true',help="If true, require LIGO frames are present (usually via CVMFS). User is responsible for generating cache file compatible with it.  This option insures that the cache file is properly transferred (because you have generated it)")
 parser.add_argument("--use-ini",default=None,type=str,help="Attempt to parse LI ini file to set corresponding options. WARNING: MAY OVERRIDE SOME OTHER COMMAND-LINE OPTIONS")
 parser.add_argument("--verbose",action='store_true')
+parser.add_argument("--force-scatter-grids",action='store_true',help="Eliminates all non-scatter intrinsic points from hyperbolic grids throughout the workflow.")
 opts=  parser.parse_args()
 
+# need --assume-hyperbolic when using --force-scatter-grids
+if opts.force_scatter_grids and not opts.assume_hyperbolic:
+    parser.error("--force-scatter-grids requires --assume-hyperbolic!")
+
 if opts.use_mtot_coords:
     if opts.force_mtot_range is None:
         print('Using the mtot coords requires a specified range!')
@@ -1272,6 +1277,8 @@ def crit_m2(delta):
         cmd += " --random-parameter eccentricity --random-parameter-range " + ecc_range_str
     if opts.assume_hyperbolic:
         cmd += f" --random-parameter E0 --random-parameter-range [{opts.E0_min},{opts.E0_max}] --random-parameter p_phi0 --random-parameter-range [{opts.pphi0_min},{opts.pphi0_max}] "
+        if opts.force_scatter_grids:
+            cmd += " --force-scatter "
     if opts.internal_tabular_eos_file:
         cmd += " --tabular-eos-file {} ".format(opts.internal_tabular_eos_file)
         grid_size *=2  # larger grids needed for discrete realization scenarios
@@ -1422,6 +1429,8 @@ def lambda_m_estimate(m):
     helper_puff_args += " --parameter eccentricity "
 if opts.assume_hyperbolic:
     helper_puff_args += " --parameter E0 --parameter p_phi0 "
+    if opts.force_scatter_grids:
+        helper_puff_args += " --force-scatter "
 
 if event_dict["MChirp"] >25:
     # at high mass, mc/eta correlation weak, don't want to have eta coordinate degeneracy at q=1 to reduce puff proposals  near there

From cf7bc1ae01ed6876adbde6b304198c8cd0a51ec0 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Mon, 20 Jan 2025 09:16:02 -0800
Subject: [PATCH 31/50] add --force-scatter to puffball

---
 .../Code/bin/util_ParameterPuffball.py             | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/MonteCarloMarginalizeCode/Code/bin/util_ParameterPuffball.py b/MonteCarloMarginalizeCode/Code/bin/util_ParameterPuffball.py
index ca979919c..d169b32f2 100755
--- a/MonteCarloMarginalizeCode/Code/bin/util_ParameterPuffball.py
+++ b/MonteCarloMarginalizeCode/Code/bin/util_ParameterPuffball.py
@@ -54,6 +54,7 @@
 parser.add_argument("--downselect-parameter-range",action='append',type=str)
 parser.add_argument("--enforce-duration-bound",default=None,type=float,help="If present, enforce a duration bound. Used to prevent grid placement for obscenely long signals, when the window size is prescribed")
 parser.add_argument("--regularize",action='store_true',help="Add some ad-hoc terms based on priors, to help with nearly-singular matricies")
+parser.add_argument('--force-scatter',default=False,action='store_true',help='For hyperbolic analyses forces only scatter grid points.')
 opts=  parser.parse_args()
 
 if opts.random_parameter is None:
@@ -261,6 +262,19 @@ def test_index_distance(indx,thresh=opts.force_away):
             val = val/ lal.MSUN_SI
         if val < downselect_dict[param][0] or val > downselect_dict[param][1]:
             include_item =False
+            
+    if opts.force_scatter:
+        if include_item==False:
+            # no need to evaluate if the point is already downselected out
+            pass            
+        else:        
+            # removes non-scatter points from the hyperbolic grid
+            hypclass = P.extract_param('hypclass')
+            if hypclass == 'scatter':
+                    include_item = True
+            else:
+                include_item = False
+    
     if include_item:
         P_out.append(P)
 

From 54d82555dd6ca12f31501535aacdf621eb8be76c Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Mon, 20 Jan 2025 09:16:37 -0800
Subject: [PATCH 32/50] fix conditional in --force-scatter

---
 .../Code/bin/util_ManualOverlapGrid.py           | 16 ++++++++++------
 1 file changed, 10 insertions(+), 6 deletions(-)

diff --git a/MonteCarloMarginalizeCode/Code/bin/util_ManualOverlapGrid.py b/MonteCarloMarginalizeCode/Code/bin/util_ManualOverlapGrid.py
index 6a2759ab6..2ccd2e68a 100755
--- a/MonteCarloMarginalizeCode/Code/bin/util_ManualOverlapGrid.py
+++ b/MonteCarloMarginalizeCode/Code/bin/util_ManualOverlapGrid.py
@@ -369,12 +369,16 @@ def evaluate_overlap_on_grid(hfbase,param_names, grid):
                 include_item =False
                 
         if opts.force_scatter:
-            # removes non-scatter points from the hyperbolic grid
-            hypclass = Pgrid.extract_param('hypclass')
-            if hypclass == 'scatter':
-                include_item = True
-            else:
-                include_item = False
+            if include_item==False:
+                # no need to evaluate if the point is already downselected out
+                pass
+            else:                
+                # removes non-scatter points from the hyperbolic grid
+                hypclass = Pgrid.extract_param('hypclass')
+                if hypclass == 'scatter':
+                    include_item = True
+                else:
+                    include_item = False
                 
         if include_item:
          grid_revised.append(line)

From 3835376266c200c8adecfeeaf208c0e0bd4655e4 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Mon, 20 Jan 2025 09:22:09 -0800
Subject: [PATCH 33/50] add --force-scatter to CIP

---
 ...nstructIntrinsicPosterior_GenericCoordinates.py | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py
index a5840dde4..7a57ec1eb 100755
--- a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py
+++ b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py
@@ -343,6 +343,7 @@ def extract_combination_from_LI(samples_LI, p):
 parser.add_argument("--assume-eos-but-primary-bh",action='store_true',help="Special case of known EOS, but primary is a BH")
 parser.add_argument("--use-eccentricity", action="store_true")
 parser.add_argument("--use-hyperbolic", action="store_true")
+parser.add_argument("--force-scatter",default=False,action='store_true', help='For hyperbolic analyses forces only scatter grid points')
 parser.add_argument("--tripwire-fraction",default=0.05,type=float,help="Fraction of nmax of iterations after which n_eff needs to be greater than 1+epsilon for a small number epsilon")
 
 # FIXME hacky options added by me (Liz) to try to get my capstone project to work.
@@ -3011,6 +3012,19 @@ def parse_corr_params(my_str):
                     if opts.verbose:
                         print(" Sample: Skipping " , line, ' due to ', p, val, downselect_dict[p])
 
+        if opts.force_scatter:
+            if include_item==False:
+                # no need to evaluate if the point is already downselected out
+                pass
+            else:
+                # removes non-scatter points from the hyperbolic grid
+                hypclass = Pgrid.extract_param('hypclass')
+                if hypclass == 'scatter':
+                    include_item = True
+                else:
+                    include_item = False
+                        
+                        
         # Set some superfluous quantities, needed only for PN approximants, so the result is generated sensibly
         Pgrid.ampO =opts.amplitude_order
         Pgrid.phaseO =opts.phase_order

From 3bfc66069a8001ee556e8c9f6e4768ec64d7ee24 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Mon, 20 Jan 2025 09:22:52 -0800
Subject: [PATCH 34/50] add --force-scatter-grids

---
 .../Code/bin/util_RIFT_pseudo_pipe.py           | 17 +++++++++++++++++
 1 file changed, 17 insertions(+)

diff --git a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py
index 908a8c958..2e1f35c9b 100755
--- a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py
+++ b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py
@@ -290,8 +290,13 @@ def unsafe_parse_arg_string_dict(my_argstr):
 parser.add_argument("--internal-mitigate-fd-J-frame",default="L_frame",help="L_frame|rotate, choose method to deal with ChooseFDWaveform being in wrong frame. Default is to request L frame for inputs")
 parser.add_argument("--first-iteration-jumpstart",action='store_true',help="No ILE jobs the first iteration.  Assumes you already have .composite files and want to get going. Particularly helpful for subdag systems")
 parser.add_argument("--use-mtot-coords",action='store_true',help="Passed to the helper to configure CIP and PUFF for mtot instead of mc.")
+parser.add_argument("--force-scatter-grids",action='store_true',help="Eliminates all non-scatter intrinsic points from hyperbolic grids throughout the workflow.") 
 opts=  parser.parse_args()
 
+# need --assume-hyperbolic when using --force-scatter-grids
+if opts.force_scatter_grids and not opts.assume_hyperbolic:
+    parser.error("--force-scatter-grids requires --assume-hyperbolic!")
+
 
 if (opts.use_ini):
     # Attempt to lazy-parse all command line arguments from ini file
@@ -658,6 +663,9 @@ def unsafe_parse_arg_string_dict(my_argstr):
     if not(opts.force_pphi0_min is None):
         pphi0_min = opts.force_pphi0_min
         cmd += " --pphi0-min {}  ".format(pphi0_min)
+        
+    if opts.force_scatter_grids:
+        cmd += " --force-scatter-grids "
 if opts.assume_highq:
     cmd+= ' --assume-highq  --force-grid-stretch-mc-factor 2'  # the mc range, tuned to equal-mass binaries, is probably too narrow. Workaround until fixed in helper
     npts_it =1000
@@ -1107,6 +1115,11 @@ def unsafe_parse_arg_string_dict(my_argstr):
         if not(opts.force_pphi0_min is None):
             pphi0_min = opts.force_pphi0_min
             line += " --pphi0-min {}  ".format(pphi0_min)
+            
+        if opts.force_scatter_grids:
+            line += " --force-scatter "
+        
+        
     if not(opts.manual_extra_cip_args is None):
         line += " {} ".format(opts.manual_extra_cip_args)  # embed with space on each side, avoid collisions
     line += "\n"
@@ -1180,6 +1193,10 @@ def unsafe_parse_arg_string_dict(my_argstr):
             pphi0_max = opts.force_pphi0_max
             pphi0_min = opts.force_pphi0_min
             puff_params += f" --downselect-parameter p_phi0 --downselect-parameter-range [{pphi0_min},{pphi0_max}]"
+            
+        if opts.force_scatter_grids:
+            puff_params += ' --force-scatter '
+            
 if opts.assume_highq:
     puff_params = puff_params.replace(' delta_mc ', ' eta ')  # use natural coordinates in the high q strategy. May want to do this always
     puff_max_it +=3

From d1c21ca771e475a48c4aa22c710d276eb350a7fb Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Mon, 20 Jan 2025 09:50:00 -0800
Subject: [PATCH 35/50] added --hyperbolic option to skip tapering

---
 .../Code/bin/util_LALWriteFrame.py                     | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/MonteCarloMarginalizeCode/Code/bin/util_LALWriteFrame.py b/MonteCarloMarginalizeCode/Code/bin/util_LALWriteFrame.py
index e1158f4c2..fdff2ef58 100755
--- a/MonteCarloMarginalizeCode/Code/bin/util_LALWriteFrame.py
+++ b/MonteCarloMarginalizeCode/Code/bin/util_LALWriteFrame.py
@@ -41,6 +41,7 @@
 parser.add_argument("--verbose", action="store_true",default=False)
 parser.add_argument("--l-max",default=4,type=float,help='Lmax number of modes')
 parser.add_argument('--gen-hlmoft', action='store_true', help='Creates hoft from hlmoft')
+parser.add_argument('--hyperbolic', action='store_true', help='skips tapering for hyperbolic waveforms')
 opts=  parser.parse_args()
 
 
@@ -52,7 +53,8 @@
     P.randomize(aligned_spin_Q=True,default_inclination=opts.incl)
     P.m1 = opts.mass1*lalsimutils.lsu_MSUN
     P.m2 = opts.mass2*lalsimutils.lsu_MSUN
-    P.taper = lalsimutils.lsu_TAPER_START
+    if not(opts.hyperbolic):
+        P.taper = lalsimutils.lsu_TAPER_START
     P.tref =1000000000  # default
     if opts.approx:
         P.approx = lalsim.GetApproximantFromString(str(opts.approx))
@@ -66,10 +68,12 @@
     xmldoc = utils.load_filename(filename, verbose = True, contenthandler =lalsimutils.cthdler)
     sim_inspiral_table = lsctables.SimInspiralTable.get_table(xmldoc)
     P.copy_sim_inspiral(sim_inspiral_table[int(event)])
-    P.taper = lalsimutils.lsu_TAPER_START
+    if not(opts.hyperbolic):
+        P.taper = lalsimutils.lsu_TAPER_START
     if opts.approx:
         P.approx = lalsim.GetApproximantFromString(str(opts.approx))
-P.taper = lalsimutils.lsu_TAPER_START  # force taper
+if not(opts.hyperbolic):
+    P.taper = lalsimutils.lsu_TAPER_START  # force taper
 P.detector = opts.instrument
 if opts.approx == "EccentricTD":
     P.phaseO = 3

From 151d47ca7cd4ab598e3a47c60bdd5d9a202c5b63 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Fri, 31 Jan 2025 06:13:56 -0800
Subject: [PATCH 36/50] update to end tapering in hlmoft

---
 MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py | 14 +++++++++++---
 1 file changed, 11 insertions(+), 3 deletions(-)

diff --git a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
index 3fcd5f8f8..dc20d7c57 100644
--- a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
+++ b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
@@ -3870,6 +3870,7 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil
             for mode in modes_used_new2:
                 hlm[mode].data.data[0:n_samp] *= vectaper
         else:
+            # determine location of start taper
             if not 'n_samp' in locals():
                 for count,value in enumerate(hlm[(2,2)].data.data):
                     if count ==0:
@@ -3877,6 +3878,15 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil
                     if np.abs(np.real(value)-np.real(hlm[(2,2)].data.data[0])) > 0.01 * np.abs(np.real(hlm[(2,2)].data.data[0])):
                         n_samp=int(count/2)
                         break
+            
+            # determine location of end taper
+            if not 'n_samp2' in locals():
+                for count, value in enumerate(reversed(hlm[(2,2)].data.data)):  # Scan backwards
+                    if count == 0:
+                        continue
+                    if np.abs(np.real(value) - np.real(hlm[(2,2)].data.data[-1])) > 0.01 * np.abs(np.real(hlm[(2,2)].data.data[-1])):
+                        n_samp2 = int(count / 2)
+                        break
                         
             # peak finding to determine system type based on 22 mode
             wf_22 = hlm[(2,2)].data.data
@@ -3908,9 +3918,8 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil
             
             if len(filtered_peaks) == 1:
                 #check if scatter or plunge            
-                if np.abs(hlm[(2,2)].data.length-nmax) > 3e3:
+                if np.abs(hlm[(2,2)].data.length-nmax) > 3e3:# NOTE - NEED TO FIX THIS!!!
                     print('Scatter waveform, tapering both ends')
-                    n_samp2=n_samp
                     vectaper2= 0.5 + 0.5 * np.cos(np.pi * np.arange(n_samp2 + 1) / (1. * n_samp2))
                     for mode in modes_used_new2:
                         hlm[mode].data.data[-(n_samp2+1):] *= vectaper2
@@ -3919,7 +3928,6 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil
                     
             elif len(filtered_peaks) ==0:
                 print('Non-interactive hyperbolic waveform, tapering both ends')
-                n_samp2=n_samp
                 vectaper2 = 0.5 + 0.5 * np.cos(np.pi * np.arange(n_samp2 + 1) / (1. * n_samp2))
                 for mode in modes_used_new2:
                         hlm[mode].data.data[-(n_samp2+1):] *= vectaper2

From fa6df13582bdb0e3715bcbebf1e8e8a4c381b756 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Tue, 4 Feb 2025 14:00:11 -0800
Subject: [PATCH 37/50] fixed conditional for opts.save_hyperbolic

---
 .../Code/bin/integrate_likelihood_extrinsic_batchmode           | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode b/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode
index 219d024c5..342d5b580 100755
--- a/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode
+++ b/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode
@@ -1634,7 +1634,7 @@ def analyze_event(P_list, indx_event, data_dict, psd_dict, fmax, opts,inv_spec_t
         if opts.save_eccentricity:
             # output format when eccentricity is being used
             numpy.savetxt(fname_output_txt, numpy.array([[event_id, m1, m2, P.s1x, P.s1y, P.s1z, P.s2x, P.s2y, P.s2z,  P.eccentricity, log_res+manual_avoid_overflow_logarithm, sqrt_var_over_res,sampler.ntotal, neff ]]))  #dict_return["convergence_test_results"]["normal_integral]"
-        if opts.save_hyperbolic:
+        elif opts.save_hyperbolic:
             # output format when hyperbolic is being used
             numpy.savetxt(fname_output_txt, numpy.array([[event_id, m1, m2, P.s1x, P.s1y, P.s1z, P.s2x, P.s2y, P.s2z,  P.E0, P.p_phi0, log_res+manual_avoid_overflow_logarithm, sqrt_var_over_res,sampler.ntotal, neff ]]))
         elif not (P.lambda1>0 or P.lambda2>0):

From 6d466f21534e91e580c8faab2b93c6c81849596e Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Fri, 7 Feb 2025 04:08:16 -0800
Subject: [PATCH 38/50] plunge/scatter class fix 1

---
 MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
index dc20d7c57..1b4ce03bb 100644
--- a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
+++ b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
@@ -961,7 +961,8 @@ def extract_param(self,p):
             if len(filtered_peaks) == 1:
                 # scatter case OR plunge case
                 
-                if np.abs(len(hptmp) - np.argmax(hptmp)) > pars['srate_interp']/5.46133: #3000 samples at 16384 [Hz]
+                #if np.abs(len(hptmp) - np.argmax(hptmp)) > pars['srate_interp']/5.46133: #3000 samples at 16384 [Hz]
+                if np.abs(hptmp)[-1] > 1e-26:
                     # scatter waveform
                     return 'scatter'
                 
@@ -3918,7 +3919,8 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil
             
             if len(filtered_peaks) == 1:
                 #check if scatter or plunge            
-                if np.abs(hlm[(2,2)].data.length-nmax) > 3e3:# NOTE - NEED TO FIX THIS!!!
+                #if np.abs(hlm[(2,2)].data.length-nmax) > 3e3:# NOTE - NEED TO FIX THIS!!!
+                if np.abs(hlm[(2,2)].data.data)[-1] > 1e-26:
                     print('Scatter waveform, tapering both ends')
                     vectaper2= 0.5 + 0.5 * np.cos(np.pi * np.arange(n_samp2 + 1) / (1. * n_samp2))
                     for mode in modes_used_new2:

From 9dab1fb779f67c397222dd349e0b56ad1904d445 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Fri, 21 Feb 2025 04:47:04 -0800
Subject: [PATCH 39/50] classifier updates

---
 .../Code/RIFT/lalsimutils.py                  | 212 +++++++++++-------
 1 file changed, 133 insertions(+), 79 deletions(-)

diff --git a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
index 1b4ce03bb..eeed7f97d 100644
--- a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
+++ b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
@@ -885,20 +885,15 @@ def extract_param(self,p):
         if p == 'mtot':
             return (self.m2+self.m1)
         if p == 'q':
-            return self.m2/self.m1
-        
-        ####################################
-        # EXPERIMENTING WITH HYPCLASS HERE #
-        
+            return self.m2/self.m1        
         if p == 'hypclass':
-            
-            
+            # Checks type of hyperbolic waveform: scatter, plunge, zoomwhirl, or meaningless
             # check if valid
             if self.E0 == 0.0:
                 print('Invalid use of hypclass: non-hyperbolic configuration')
                 return None
             
-            # run the classifier
+            # Generate waveform
             pars = {
                 'M'                  : (self.m1+self.m2)/lal.MSUN_SI,
                 'q'                  : self.m1/self.m2,
@@ -915,8 +910,8 @@ def extract_param(self,p):
                 'nqc'                : 2, # sets the NQCs, 2=no
                 'nqc_coefs_hlm'      : 0, # Option for the NQC model used in the waveform. 0=none
                 'nqc_coefs_flx'      : 0, # Option for the NQC model used in the flux. 0=none
-                'use_mode_lm'        : [1,-1], # 2\pm 2 modes
-                'output_lm'          : [1,-1],
+                'use_mode_lm'        : [1], # 22 mode
+                'output_lm'          : [1],
                 'srate_interp'       : 1./self.deltaT,
                 'use_geometric_units': 0,
                 'interp_uniform_grid': 1,
@@ -927,16 +922,22 @@ def extract_param(self,p):
                 'output_hpc'         : 0 # output plus and cross polarizations, 0=no
             }
             
-            #print('Classifying hyperbolic waveform...')
-            #print(pars)
-            
             t, hptmp, hctmp, hlmtmp, dym = EOBRun_module.EOBRunPy(pars)
             
             # peak finding to classifiy
             
-            # wf amplitude
-            amp = np.sqrt(np.abs(hptmp)**2+np.abs(hctmp)**2)
-            amp_norm = amp / np.amax(amp) # normalize amplitude for peak finding                
+            # wf amplitude - 22 mode only
+            tmp_22 = np.array(hlmtmp['1'])
+            tmp_22[0] *= (m_total_s/distance_s)*nu
+            amp = np.abs(tmp_22[0] * np.exp(-1j*(2*(np.pi/2.)+tmp_22[1])))
+            amp_norm = amp / np.amax(amp) # normalize amplitude for peak finding
+            
+            # Check for cases where the amplitude is max at the start or end
+            if np.argmax(amp_norm) == 0 or np.argmax(amp_norm) == len(amp_norm) - 1:
+                reclassify = True
+            else:
+                reclassify = False
+            
             # peak finding to determine system type
             height_thresh = 0.25
             prom_thresh = 0.1
@@ -959,29 +960,47 @@ def extract_param(self,p):
             
             # parsing number of peaks after filtering against distance tolerance
             if len(filtered_peaks) == 1:
-                # scatter case OR plunge case
-                
-                #if np.abs(len(hptmp) - np.argmax(hptmp)) > pars['srate_interp']/5.46133: #3000 samples at 16384 [Hz]
-                if np.abs(hptmp)[-1] > 1e-26:
+                if np.abs(amp)[-1] > 1e-26:
                     # scatter waveform
-                    return 'scatter'
-                
+                    return 'scatter'                
                 else:
                     # plunge waveform
                     return 'plunge'
                 
             elif len(filtered_peaks) == 0:
                 # meaningless waveform
-                return None
+                reclassify = True
             
             else:
                 # zoom whirl waveform
-                return 'zoomwhirl'                
-        
-        
-        ####################################
-        
-        
+                return 'zoomwhirl'
+            
+            if reclassify:
+                print('Running re-classification')
+                # run minimal threshold peak finder
+                all_peaks, all_props = signal.find_peaks(amp_norm, height=0.0001, prominence=0.0001)
+                
+                # checking for minima if no peaks found                
+                if len(all_peaks) == 0:
+                    print("No peaks found, checking for minima instead...")
+                    all_peaks, all_props = signal.find_peaks((-1*amp_norm + 1.0), height=0.0001, prominence=0.0001)
+                
+                if len(all_props['prominences']) > 3:
+                    print('MANY peaks detected on reclassification, evaluating...')
+                    # these can be scatter or plunge
+                    if np.abs(amp)[-1] < 1e-26:
+                        print('Reclassifying to Plunge')
+                        return 'plunge'
+                    else:
+                        print('Reclassifying to Scatter')
+                        return 'scatter'
+                elif len(all_props['prominences']) == 3 or len(all_props['prominences']) == 2 or len(all_props['prominences']) == 1:
+                    # these are always scaters
+                    print('Reclassifying to Scatter')
+                    return 'scatter'
+                else:
+                    print('No peaks detected after reclassifcation')
+                    return 'meaningless'
         if p == 'delta' or p=='delta_mc':  # Same access routine
             return (self.m1-self.m2)/(self.m1+self.m2)
         if p == 'mc':
@@ -3657,6 +3676,8 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil
         elif (P.eccentricity == 0.0):
             print("Using hyperbolic call RIFT O4b branch")
             hyp_wav = True # convenient way to know if the waveform is hyperbolic
+            #print('NOTE! Forcing just the 22 mode!')
+            #k = [1] # forcing just the 22 mode
             pars = {
                 'M'                  : M1+M2,
                 'q'                  : M1/M2,
@@ -3722,12 +3743,21 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil
             ## Set the epoch for non-hyperbolic cases ##
             hpepoch = -P.deltaT*np.argmax(np.abs(hptmp)**2+np.abs(hctmp)**2)            
         else:
-            ## custom epoch for the hyperbolic case ##            
-            # wf amplitude
-            amp = np.sqrt(np.abs(hptmp)**2+np.abs(hctmp)**2)
+            ## custom epoch for the hyperbolic case ##
+                        
+            # wf amplitude - 22 mode only
+            tmp_22 = np.array(hlmtmp['1'])
+            tmp_22[0] *= (m_total_s/distance_s)*nu
+            amp = np.abs(tmp_22[0] * np.exp(-1j*(2*(np.pi/2.)+tmp_22[1])))
             amp_norm = amp / np.amax(amp) # normalize amplitude for peak finding
             
-            # peak finding to determine system type
+            # Check for cases where the amplitude is max at the start or end
+            if np.argmax(amp_norm) == 0 or np.argmax(amp_norm) == len(amp_norm) - 1:
+                reclassify = True
+            else:
+                reclassify = False
+            
+            # initial peak finding to determine system type
             height_thresh = 0.25
             prom_thresh = 0.1
             peaks, props = signal.find_peaks(amp_norm, height = height_thresh, prominence = prom_thresh) 
@@ -3749,16 +3779,61 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil
             # parsing number of peaks after filtering against distance tolerance
             if len(filtered_peaks) == 1:
                 # scatter case OR plunge case, we can set the epoch normally
-                hpepoch = -P.deltaT*np.argmax(np.abs(hptmp)**2+np.abs(hctmp)**2)
+                hpepoch = -P.deltaT*np.argmax(amp)
+                
+                if np.abs(amp)[-1] > 1e-26:
+                    hypclass = 'scatter' # maybe should do through assign_param?
+                else:
+                    hypclass = 'plunge'
+                
+                
             elif len(filtered_peaks) == 0:
-                # meaningless waveform, essentially a non-interacting case. Set the epoch normally
-                # These points should return very low likelihood and not interfere with the analysis
-                print('WARNING: no peak detected; non-interacting hyperbolic case')
-                hpepoch = -P.deltaT*np.argmax(np.abs(hptmp)**2+np.abs(hctmp)**2)                
+                hypclass = 'meaningless'
+                hpepoch = -P.deltaT*np.argmax(amp)
+                reclassify = True
+                
             else:
                 # capture case, we need to force the epoch to be the last peak
+                hypclass = 'zoomwhirl'
                 hpepoch = -P.deltaT*filtered_peaks[-1]
                 
+            if reclassify:
+                print('Running re-classification')
+                # run minimal threshold peak finder
+                all_peaks, all_props = signal.find_peaks(amp_norm, height=0.000001, prominence=0.000001)
+                
+                # checking for minima if no peaks found
+                
+                if len(all_peaks) == 0:
+                    print("No peaks found, checking for minima instead...")
+                    all_peaks, all_props = signal.find_peaks((-1*amp_norm + 1.0), height=0.000001, prominence=0.000001)
+                
+                if len(all_props['prominences']) > 3:
+                    print('MANY peaks detected on reclassification, evaluating...')
+                    
+                    if np.abs(amp)[-1] < 1e-26:
+                        print('Reclassifying to Plunge')
+                        hpepoch = -P.deltaT*np.argmax(amp)
+                        hypclass = 'plunge' 
+                    else:
+                        print('Reclassifying to Scatter')
+                        max_peak_index = all_peaks[np.argmax(all_props['peak_heights'])]
+                        print(f"Largest peak is at index {max_peak_index} with height {amp_norm[max_peak_index]}")
+                        # setting the epoch to the largest amplitude peak
+                        hpepoch = -P.deltaT*max_peak_index
+                elif len(all_props['prominences']) == 3 or len(all_props['prominences']) == 2 or len(all_props['prominences']) == 1:
+                    print('Reclassifying to Scatter')
+                    hypclass = 'scatter'
+                    
+                    max_peak_index = all_peaks[np.argmax(all_props['peak_heights'])]
+                    print(f"Largest peak is at index {max_peak_index} with height {amp_norm[max_peak_index]}")
+                    
+                    # setting the epoch to the largest amplitude peak
+                    hpepoch = -P.deltaT*max_peak_index
+                else:
+                    print('No peaks detected after reclassifcation')
+                    hypclass='meaningless'
+                
         hlmlen = len(hptmp)
         hlm = {}
         hlmtmp2 = {}
@@ -3886,56 +3961,35 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil
                     if count == 0:
                         continue
                     if np.abs(np.real(value) - np.real(hlm[(2,2)].data.data[-1])) > 0.01 * np.abs(np.real(hlm[(2,2)].data.data[-1])):
-                        n_samp2 = int(count / 2)
+                        n_samp2 = int(count / 2) 
                         break
-                        
-            # peak finding to determine system type based on 22 mode
-            wf_22 = hlm[(2,2)].data.data
-            wf_22_norm = wf_22 / np.amax(wf_22)
-            height_thresh = 0.25
-            prom_thresh = 0.1
-            peaks, props = signal.find_peaks(wf_22_norm, height = height_thresh, prominence = prom_thresh)
-            peak_heights = props['peak_heights']
-            indices_to_keep = set()
-            sorted_indices = np.argsort(peak_heights)[::-1]
-            tol = int(pars['srate_interp'] / 13.65) # 300 samples at srate of 4096 - minimum distance between peaks.
-            for i in sorted_indices:
-                peak = peaks[i]                
-                keep = True
-                for kept_index in indices_to_keep:
-                    if abs(peaks[kept_index] - peak) <= tol:
-                        keep = False
-                        break                
-                if keep:
-                    indices_to_keep.add(i)                    
-            filtered_peaks = peaks[list(indices_to_keep)]
             
             # Always taper the start
             
             vectaper= 0.5 + 0.5*np.cos(np.pi* (1-np.arange(n_samp)/(1.*n_samp)))
             nmax = np.argmax(hlm[(2,2)].data.data)
             for mode in modes_used_new2:
+                #pass
                 hlm[mode].data.data[0:n_samp] *= vectaper
-            
-            if len(filtered_peaks) == 1:
-                #check if scatter or plunge            
-                #if np.abs(hlm[(2,2)].data.length-nmax) > 3e3:# NOTE - NEED TO FIX THIS!!!
-                if np.abs(hlm[(2,2)].data.data)[-1] > 1e-26:
-                    print('Scatter waveform, tapering both ends')
-                    vectaper2= 0.5 + 0.5 * np.cos(np.pi * np.arange(n_samp2 + 1) / (1. * n_samp2))
-                    for mode in modes_used_new2:
-                        hlm[mode].data.data[-(n_samp2+1):] *= vectaper2
-                else:
-                    print('Plunge waveform, only start taper')
-                    
-            elif len(filtered_peaks) ==0:
-                print('Non-interactive hyperbolic waveform, tapering both ends')
-                vectaper2 = 0.5 + 0.5 * np.cos(np.pi * np.arange(n_samp2 + 1) / (1. * n_samp2))
+                
+            if hypclass == 'scatter':
+                # Taper for scatter
+                vectaper2= 0.5 + 0.5 * np.cos(np.pi * np.arange(n_samp2 + 1) / (1. * n_samp2))
                 for mode in modes_used_new2:
-                        hlm[mode].data.data[-(n_samp2+1):] *= vectaper2
-                        
-            else:
-                print('Zoom-whirl waveform, only start taper')                        
+                    print(f'n_samp2 is {n_samp2} for mode {mode}')
+                    hlm[mode].data.data[-(n_samp2+1):] *= vectaper2
+            elif hypclass == 'plunge':
+                # taper for plunge
+                print('Plunge waveform, only start taper')
+            elif hypclass == 'zoomwhirl':
+                # taper for ZW
+                print('Zoom-whirl waveform, only start taper')
+            elif hypclass =='meaningless':
+                # zero out meaningless
+                for mode in modes_used_new2:
+                    hlm[mode].data.data *= 0.0
+                    
+        # end off custom hyperbolic stuff
             
         for mode in modes_used_new2:
             if not (P.deltaF is None):

From b7f70f9e0a39589b5efa126a968c815ccd69357b Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Fri, 21 Feb 2025 05:03:52 -0800
Subject: [PATCH 40/50] add verbose option

---
 MonteCarloMarginalizeCode/Code/bin/util_FrameZeroNoiseSNR.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/MonteCarloMarginalizeCode/Code/bin/util_FrameZeroNoiseSNR.py b/MonteCarloMarginalizeCode/Code/bin/util_FrameZeroNoiseSNR.py
index 371ec6548..6d95abb38 100755
--- a/MonteCarloMarginalizeCode/Code/bin/util_FrameZeroNoiseSNR.py
+++ b/MonteCarloMarginalizeCode/Code/bin/util_FrameZeroNoiseSNR.py
@@ -15,6 +15,7 @@
 parser.add_argument("--fmin-snr",  type=float, default=20)
 parser.add_argument("--fmax-snr", type=float, default=1500)
 parser.add_argument("--plot-sanity",action='store_true')
+parser.add_argument("--verbose", default=False, action='store_true')
 opts=  parser.parse_args()
 
 fminSNR = opts.fmin_snr
@@ -48,7 +49,7 @@
     else:
         analyticPSD_Q=False
         print("Reading PSD for instrument %s from %s" % (ifo, psd_name[ifo]))
-        psd_dict[ifo] = lalsimutils.load_resample_and_clean_psd(psd_name[ifo], ifo, df)
+        psd_dict[ifo] = lalsimutils.load_resample_and_clean_psd(psd_name[ifo], ifo, df, verbose=opts.verbose)
     IP = lalsimutils.ComplexIP(fLow=fminSNR, fNyq=fSample/2,deltaF=df,psd=psd_dict[ifo],fMax=fmaxSNR,analyticPSD_Q=analyticPSD_Q)
     rhoDet = rho_dict[ifo] = IP.norm(data_dict[ifo])
     print(ifo, rho_dict[ifo])

From 4a8697f1a3b4b19b6d056246e9d445f0a24b1755 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Fri, 21 Feb 2025 05:14:18 -0800
Subject: [PATCH 41/50] zero window and split seglen

---
 .../Code/bin/helper_LDG_Events.py                | 16 +++++++++++-----
 1 file changed, 11 insertions(+), 5 deletions(-)

diff --git a/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py b/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
index f2da3e33f..6a443d210 100755
--- a/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
+++ b/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
@@ -1175,11 +1175,17 @@ def crit_m2(delta):
     helper_ile_args += " --data-start-time " + str(data_start_time) + " --data-end-time " + str(data_end_time)  + " --inv-spec-trunc-time 0 --window-shape " + str(window_shape)
 elif opts.data_LI_seglen:
     seglen = opts.data_LI_seglen
-    # Use LI-style positioning of trigger relative to 2s before end of buffer
-    # Use LI-style tukey windowing
-    window_shape = 0.4*2/seglen
-    data_end_time = event_dict["tref"]+2
-    data_start_time = event_dict["tref"] +2 - seglen
+    if opts.assume_hyperbolic:
+        window_shape = 0.0 # DO NOT window at all
+        # split seglen across the event time
+        data_start_time = event_dict["tref"] - seglen/2
+        data_end_time = event_dict["tref"] + seglen/2
+    else:
+        # Use LI-style positioning of trigger relative to 2s before end of buffer
+        # Use LI-style tukey windowing
+        window_shape = 0.4*2/seglen
+        data_end_time = event_dict["tref"]+2
+        data_start_time = event_dict["tref"] +2 - seglen
     helper_ile_args += " --data-start-time " + str(data_start_time) + " --data-end-time " + str(data_end_time)  + " --inv-spec-trunc-time 0  --window-shape " + str(window_shape)
 if opts.assume_eccentric:
     helper_ile_args += " --save-eccentricity "

From 3281ce68578c1a4dc9790efacde8e70dfa5604bf Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Fri, 21 Feb 2025 07:08:44 -0800
Subject: [PATCH 42/50] hypclass unit fix

---
 MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
index eeed7f97d..f61fe1ff7 100644
--- a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
+++ b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
@@ -928,6 +928,11 @@ def extract_param(self,p):
             
             # wf amplitude - 22 mode only
             tmp_22 = np.array(hlmtmp['1'])
+            distance_s = self.dist/lal.C_SI
+            m_total_s = MsunInSec*(self.m1+self.m2)/lal.MSUN_SI
+            M1=self.m1/lal.MSUN_SI
+            M2=self.m2/lal.MSUN_SI
+            nu=M1*M2/((M1+M2)**2)
             tmp_22[0] *= (m_total_s/distance_s)*nu
             amp = np.abs(tmp_22[0] * np.exp(-1j*(2*(np.pi/2.)+tmp_22[1])))
             amp_norm = amp / np.amax(amp) # normalize amplitude for peak finding

From b27bfe8af5afe43dd1302fa803cd1183b63a5b72 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Fri, 21 Feb 2025 07:11:14 -0800
Subject: [PATCH 43/50] add option to force 22 mode only

---
 MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py     |  8 ++++++--
 .../Code/RIFT/likelihood/factored_likelihood.py        | 10 +++++++---
 .../Code/bin/helper_LDG_Events.py                      |  3 +++
 .../Code/bin/integrate_likelihood_extrinsic_batchmode  |  3 ++-
 .../Code/bin/util_LALWriteFrame.py                     |  6 +++++-
 .../Code/bin/util_RIFT_pseudo_pipe.py                  |  6 +++++-
 6 files changed, 28 insertions(+), 8 deletions(-)

diff --git a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
index f61fe1ff7..b44f2df02 100644
--- a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
+++ b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
@@ -3681,8 +3681,9 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil
         elif (P.eccentricity == 0.0):
             print("Using hyperbolic call RIFT O4b branch")
             hyp_wav = True # convenient way to know if the waveform is hyperbolic
-            #print('NOTE! Forcing just the 22 mode!')
-            #k = [1] # forcing just the 22 mode
+            if kwargs.get('force_22_mode', False):
+                print('Forcing ONLY the 22 modes')
+                k = [1]
             pars = {
                 'M'                  : M1+M2,
                 'q'                  : M1/M2,
@@ -4015,6 +4016,9 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil
 #        for mode in hlm:
 #            hlm[mode].data.data*= hp.data.data
 
+        print('Actual modes used:')
+        print(hlm.keys())
+
         return hlm
     else: # (P.approx == lalSEOBv4 or P.approx == lalsim.SEOBNRv2 or P.approx == lalsim.SEOBNRv1 or  P.approx == lalsim.EOBNRv2 
         extra_params = P.to_lal_dict_extended(extra_args_dict=extra_waveform_args)
diff --git a/MonteCarloMarginalizeCode/Code/RIFT/likelihood/factored_likelihood.py b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/factored_likelihood.py
index 6b3b12224..df25912b9 100644
--- a/MonteCarloMarginalizeCode/Code/RIFT/likelihood/factored_likelihood.py
+++ b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/factored_likelihood.py
@@ -121,7 +121,7 @@ def PrecomputeLikelihoodTerms(event_time_geo, t_window, P, data_dict,
         extra_waveform_kwargs={},
         use_gwsignal=False,
         use_gwsignal_approx=None,
-       use_external_EOB=False,nr_lookup=False,nr_lookup_valid_groups=None,no_memory=True,perturbative_extraction=False,perturbative_extraction_full=False,hybrid_use=False,hybrid_method='taper_add',use_provided_strain=False,ROM_group=None,ROM_param=None,ROM_use_basis=False,ROM_limit_basis_size=None,skip_interpolation=False):
+       use_external_EOB=False,nr_lookup=False,nr_lookup_valid_groups=None,no_memory=True,perturbative_extraction=False,perturbative_extraction_full=False,hybrid_use=False,hybrid_method='taper_add',use_provided_strain=False,ROM_group=None,ROM_param=None,ROM_use_basis=False,ROM_limit_basis_size=None,skip_interpolation=False,**kwargs):
     """
     Compute < h_lm(t) | d > and < h_lm | h_l'm' >
 
@@ -142,6 +142,10 @@ def PrecomputeLikelihoodTerms(event_time_geo, t_window, P, data_dict,
     rholms_intp = {}
     crossTerms = {}
     crossTermsV = {}
+    
+    if kwargs.get('force_22_mode', False):
+        # forcing 22 mode only for hyperbolic case
+        force_22_mode = True
 
     # Compute hlms at a reference distance, distance scaling is applied later
     P.dist = distMpcRef*1e6*lsu.lsu_PC
@@ -198,7 +202,7 @@ def PrecomputeLikelihoodTerms(event_time_geo, t_window, P, data_dict,
     elif use_gwsignal and (has_GWS):  # this MUST be called first, so the P.approx is never tested
         if not quiet:
             print( "  FACTORED LIKELIHOOD WITH hlmoff (GWsignal) " )            
-        hlms, hlms_conj = rgws.std_and_conj_hlmoff(P,Lmax,approx_string=use_gwsignal_approx,**extra_waveform_kwargs)
+        hlms, hlms_conj = rgws.std_and_conj_hlmoff(P,Lmax,approx_string=use_gwsignal_approx,force_22_mode=force_22_mode,**extra_waveform_kwargs)
 
     elif (not nr_lookup) and (not NR_group) and ( P.approx ==lalsim.SEOBNRv2 or P.approx == lalsim.SEOBNRv1 or P.approx==lalsim.SEOBNRv3 or P.approx == lsu.lalSEOBv4 or P.approx ==lsu.lalSEOBNRv4HM or P.approx == lalsim.EOBNRv2 or P.approx == lsu.lalTEOBv2 or P.approx==lsu.lalTEOBv4 ):
         # note: alternative to this branch is to call hlmoff, which will actually *work* if ChooseTDModes is propertly implemented for that model
@@ -250,7 +254,7 @@ def PrecomputeLikelihoodTerms(event_time_geo, t_window, P, data_dict,
         #         hlms_conj = lsu.SphHarmFrequencySeries_to_dict(hlms_conj_list, Lmax) # a dictionary
         # else:
         #         hlms_conj = hlms_conj_list
-        hlms, hlms_conj = lsu.std_and_conj_hlmoff(P,Lmax,**extra_waveform_kwargs)
+        hlms, hlms_conj = lsu.std_and_conj_hlmoff(P,Lmax,force_22_mode=force_22_mode,**extra_waveform_kwargs)
     elif (nr_lookup or NR_group) and useNR:
 	    # look up simulation
 	    # use nrwf to get hlmf
diff --git a/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py b/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
index 6a443d210..46f3ef5b2 100755
--- a/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
+++ b/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
@@ -253,6 +253,7 @@ def get_observing_run(t):
 parser.add_argument("--use-ini",default=None,type=str,help="Attempt to parse LI ini file to set corresponding options. WARNING: MAY OVERRIDE SOME OTHER COMMAND-LINE OPTIONS")
 parser.add_argument("--verbose",action='store_true')
 parser.add_argument("--force-scatter-grids",action='store_true',help="Eliminates all non-scatter intrinsic points from hyperbolic grids throughout the workflow.")
+parser.add_argument('--force-hyperbolic-22', action='store_true', help='Forces just the 22 modes for hyperbolic waveforms')
 opts=  parser.parse_args()
 
 # need --assume-hyperbolic when using --force-scatter-grids
@@ -1191,6 +1192,8 @@ def crit_m2(delta):
     helper_ile_args += " --save-eccentricity "
 if opts.assume_hyperbolic:
     helper_ile_args += " --save-hyperbolic "
+    if opts.force_hyperbolic_22:
+        helper_ile_args += " --force-hyperbolic-22 "
 if opts.propose_initial_grid_fisher: # and (P.extract_param('mc')/lal.MSUN_SI < 10.):
     cmd  = "util_AnalyticFisherGrid.py  --inj-file-out  proposed-grid  "
     # Add standard downselects : do not have m1, m2 be less than 1
diff --git a/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode b/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode
index 342d5b580..9829927e7 100755
--- a/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode
+++ b/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode
@@ -241,6 +241,7 @@ optp.add_option("--internal-waveform-extra-lalsuite-args",type=str,default=None)
 optp.add_option("--verbose",action='store_true')
 optp.add_option("--save-eccentricity", action="store_true")
 optp.add_option("--save-hyperbolic", action="store_true")
+optp.add_option('--force-hyperbolic-22', default=False, action='store_true', help='Forces just the 22 modes for hyperbolic waveforms')
 #
 # Add the integration options
 #
@@ -1187,7 +1188,7 @@ def analyze_event(P_list, indx_event, data_dict, psd_dict, fmax, opts,inv_spec_t
             NR_group=NR_template_group,NR_param=NR_template_param,
             use_gwsignal=opts.use_gwsignal,
             use_gwsignal_approx=opts.approximant,
-            use_external_EOB=opts.use_external_EOB,nr_lookup=opts.nr_lookup,nr_lookup_valid_groups=opts.nr_lookup_group,perturbative_extraction=opts.nr_perturbative_extraction,perturbative_extraction_full=opts.nr_perturbative_extraction_full,use_provided_strain=opts.nr_use_provided_strain,hybrid_use=opts.nr_hybrid_use,hybrid_method=opts.nr_hybrid_method,ROM_group=opts.rom_group,ROM_param=opts.rom_param,ROM_use_basis=opts.rom_use_basis,verbose=opts.verbose,quiet=not opts.verbose,ROM_limit_basis_size=opts.rom_limit_basis_size_to,no_memory=opts.no_memory,skip_interpolation=opts.vectorized, extra_waveform_kwargs=extra_waveform_kwargs)
+            use_external_EOB=opts.use_external_EOB,nr_lookup=opts.nr_lookup,nr_lookup_valid_groups=opts.nr_lookup_group,perturbative_extraction=opts.nr_perturbative_extraction,perturbative_extraction_full=opts.nr_perturbative_extraction_full,use_provided_strain=opts.nr_use_provided_strain,hybrid_use=opts.nr_hybrid_use,hybrid_method=opts.nr_hybrid_method,ROM_group=opts.rom_group,ROM_param=opts.rom_param,ROM_use_basis=opts.rom_use_basis,verbose=opts.verbose,quiet=not opts.verbose,ROM_limit_basis_size=opts.rom_limit_basis_size_to,no_memory=opts.no_memory,skip_interpolation=opts.vectorized, extra_waveform_kwargs=extra_waveform_kwargs, force_22_mode=opts.force_hyperbolic_22)
 
     if opts.auto_logarithm_offset and guess_snr:
       # important: this only impacts *this* analysis
diff --git a/MonteCarloMarginalizeCode/Code/bin/util_LALWriteFrame.py b/MonteCarloMarginalizeCode/Code/bin/util_LALWriteFrame.py
index fdff2ef58..a1a0091bd 100755
--- a/MonteCarloMarginalizeCode/Code/bin/util_LALWriteFrame.py
+++ b/MonteCarloMarginalizeCode/Code/bin/util_LALWriteFrame.py
@@ -42,6 +42,7 @@
 parser.add_argument("--l-max",default=4,type=float,help='Lmax number of modes')
 parser.add_argument('--gen-hlmoft', action='store_true', help='Creates hoft from hlmoft')
 parser.add_argument('--hyperbolic', action='store_true', help='skips tapering for hyperbolic waveforms')
+parser.add_argument('--force-hyperbolic-22', action='store_true', help='Forces just the 22 modes for hyperbolic waveforms')
 opts=  parser.parse_args()
 
 
@@ -92,7 +93,10 @@
 
 # Generate signal
 if opts.gen_hlmoft:
-    hlmT = lalsimutils.hlmoft(P, Lmax=opts.l_max)
+    if opts.hyperbolic and opts.force_hyperbolic_22:
+        hlmT = lalsimutils.hlmoft(P, Lmax=opts.l_max, force_22_mode=True)
+    else:
+        hlmT = lalsimutils.hlmoft(P, Lmax=opts.l_max)
     hoft = lalsimutils.hoft_from_hlm(hlmT, P, return_complex=False)
 else:
     hoft = lalsimutils.hoft(P,Lmax=opts.l_max)   # include translation of source, but NOT interpolation onto regular time grid
diff --git a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py
index 2e1f35c9b..490980fd7 100755
--- a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py
+++ b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py
@@ -290,7 +290,8 @@ def unsafe_parse_arg_string_dict(my_argstr):
 parser.add_argument("--internal-mitigate-fd-J-frame",default="L_frame",help="L_frame|rotate, choose method to deal with ChooseFDWaveform being in wrong frame. Default is to request L frame for inputs")
 parser.add_argument("--first-iteration-jumpstart",action='store_true',help="No ILE jobs the first iteration.  Assumes you already have .composite files and want to get going. Particularly helpful for subdag systems")
 parser.add_argument("--use-mtot-coords",action='store_true',help="Passed to the helper to configure CIP and PUFF for mtot instead of mc.")
-parser.add_argument("--force-scatter-grids",action='store_true',help="Eliminates all non-scatter intrinsic points from hyperbolic grids throughout the workflow.") 
+parser.add_argument("--force-scatter-grids",action='store_true',help="Eliminates all non-scatter intrinsic points from hyperbolic grids throughout the workflow.")
+parser.add_argument("--force-hyperbolic-22", action='store_true', help='Forces just the 22 modes for hyperbolic waveforms')
 opts=  parser.parse_args()
 
 # need --assume-hyperbolic when using --force-scatter-grids
@@ -666,6 +667,9 @@ def unsafe_parse_arg_string_dict(my_argstr):
         
     if opts.force_scatter_grids:
         cmd += " --force-scatter-grids "
+        
+    if opts.force_hyperbolic_22:
+        cmd += " --force-hyperbolic-22 "
 if opts.assume_highq:
     cmd+= ' --assume-highq  --force-grid-stretch-mc-factor 2'  # the mc range, tuned to equal-mass binaries, is probably too narrow. Workaround until fixed in helper
     npts_it =1000

From 87434b4e6ea97de501dd8ce60c0fd139ee92febe Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Fri, 21 Feb 2025 07:17:46 -0800
Subject: [PATCH 44/50] remove unnecessary print statement

---
 MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
index b44f2df02..93c3ecfba 100644
--- a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
+++ b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
@@ -3982,7 +3982,6 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil
                 # Taper for scatter
                 vectaper2= 0.5 + 0.5 * np.cos(np.pi * np.arange(n_samp2 + 1) / (1. * n_samp2))
                 for mode in modes_used_new2:
-                    print(f'n_samp2 is {n_samp2} for mode {mode}')
                     hlm[mode].data.data[-(n_samp2+1):] *= vectaper2
             elif hypclass == 'plunge':
                 # taper for plunge

From 92340e6a13269ec9eb2ca0467aa37bf22f4403f2 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Fri, 21 Feb 2025 08:54:45 -0800
Subject: [PATCH 45/50] added --force-plunge grids and --force-zoomwhirl grids

---
 .../Code/bin/helper_LDG_Events.py             | 21 +++++++++++--
 ...ctIntrinsicPosterior_GenericCoordinates.py | 26 ++++++++++++++++
 .../Code/bin/util_ManualOverlapGrid.py        | 31 +++++++++++++++++++
 .../Code/bin/util_ParameterPuffball.py        | 30 ++++++++++++++++++
 .../Code/bin/util_RIFT_pseudo_pipe.py         | 31 +++++++++++++++++--
 5 files changed, 133 insertions(+), 6 deletions(-)

diff --git a/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py b/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
index 46f3ef5b2..c53c36b15 100755
--- a/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
+++ b/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
@@ -253,12 +253,19 @@ def get_observing_run(t):
 parser.add_argument("--use-ini",default=None,type=str,help="Attempt to parse LI ini file to set corresponding options. WARNING: MAY OVERRIDE SOME OTHER COMMAND-LINE OPTIONS")
 parser.add_argument("--verbose",action='store_true')
 parser.add_argument("--force-scatter-grids",action='store_true',help="Eliminates all non-scatter intrinsic points from hyperbolic grids throughout the workflow.")
+parser.add_argument("--force-plunge-grids",action='store_true',help="Eliminates all non-plunge intrinsic points from hyperbolic grids throughout the workflow.")
+parser.add_argument("--force-zoomwhirl-grids",action='store_true',help="Eliminates all non-zoomwhirl intrinsic points from hyperbolic grids throughout the workflow.")
 parser.add_argument('--force-hyperbolic-22', action='store_true', help='Forces just the 22 modes for hyperbolic waveforms')
 opts=  parser.parse_args()
 
-# need --assume-hyperbolic when using --force-scatter-grids
-if opts.force_scatter_grids and not opts.assume_hyperbolic:
-    parser.error("--force-scatter-grids requires --assume-hyperbolic!")
+# Ensure --assume-hyperbolic is set when using any --force-X-grids option
+# Ensure only ONE of the --force-X-grids options is set
+force_grids = [opts.force_scatter_grids, opts.force_plunge_grids, opts.force_zoomwhirl_grids]
+if any(force_grids) and not opts.assume_hyperbolic:
+    parser.error("Using --force-scatter-grids, --force-plunge-grids, or --force-zoomwhirl-grids requires --assume-hyperbolic!")
+
+if sum(bool(x) for x in force_grids) > 1:
+    parser.error("CANNOT use multiple --force-X-grids options at the same time!")
 
 if opts.use_mtot_coords:
     if opts.force_mtot_range is None:
@@ -1288,6 +1295,10 @@ def crit_m2(delta):
         cmd += f" --random-parameter E0 --random-parameter-range [{opts.E0_min},{opts.E0_max}] --random-parameter p_phi0 --random-parameter-range [{opts.pphi0_min},{opts.pphi0_max}] "
         if opts.force_scatter_grids:
             cmd += " --force-scatter "
+        if opts.force_plunge_grids:
+            cmd += " --force-plunge "
+        if opts.force_zoomwhirl_grids:
+            cmd += " --force-zoomwhirl "
     if opts.internal_tabular_eos_file:
         cmd += " --tabular-eos-file {} ".format(opts.internal_tabular_eos_file)
         grid_size *=2  # larger grids needed for discrete realization scenarios
@@ -1440,6 +1451,10 @@ def lambda_m_estimate(m):
     helper_puff_args += " --parameter E0 --parameter p_phi0 "
     if opts.force_scatter_grids:
         helper_puff_args += " --force-scatter "
+    if opts.force_plunge_grids:
+        helper_puff_args += " --force-plunge "
+    if opts.force_zoomwhirl_grids:
+        helper_puff_args += " --force-zoomwhirl "
 
 if event_dict["MChirp"] >25:
     # at high mass, mc/eta correlation weak, don't want to have eta coordinate degeneracy at q=1 to reduce puff proposals  near there
diff --git a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py
index 7a57ec1eb..7d2303356 100755
--- a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py
+++ b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py
@@ -344,6 +344,8 @@ def extract_combination_from_LI(samples_LI, p):
 parser.add_argument("--use-eccentricity", action="store_true")
 parser.add_argument("--use-hyperbolic", action="store_true")
 parser.add_argument("--force-scatter",default=False,action='store_true', help='For hyperbolic analyses forces only scatter grid points')
+parser.add_argument("--force-plunge",default=False,action='store_true', help='For hyperbolic analyses forces only plunge grid points')
+parser.add_argument("--force-zoomwhirl",default=False,action='store_true', help='For hyperbolic analyses forces only zoomwhirl grid points')
 parser.add_argument("--tripwire-fraction",default=0.05,type=float,help="Fraction of nmax of iterations after which n_eff needs to be greater than 1+epsilon for a small number epsilon")
 
 # FIXME hacky options added by me (Liz) to try to get my capstone project to work.
@@ -3023,6 +3025,30 @@ def parse_corr_params(my_str):
                     include_item = True
                 else:
                     include_item = False
+
+        if opts.force_plunge:
+            if include_item==False:
+                # no need to evaluate if the point is already downselected out
+                pass
+            else:
+                # removes non-plunge points from the hyperbolic grid
+                hypclass = Pgrid.extract_param('hypclass')
+                if hypclass == 'plunge':
+                    include_item = True
+                else:
+                    include_item = False
+
+        if opts.force_zoomwhirl:
+            if include_item==False:
+                # no need to evaluate if the point is already downselected out
+                pass
+            else:
+                # removes non-zoomwhirl points from the hyperbolic grid
+                hypclass = Pgrid.extract_param('hypclass')
+                if hypclass == 'zoomwhirl':
+                    include_item = True
+                else:
+                    include_item = False
                         
                         
         # Set some superfluous quantities, needed only for PN approximants, so the result is generated sensibly
diff --git a/MonteCarloMarginalizeCode/Code/bin/util_ManualOverlapGrid.py b/MonteCarloMarginalizeCode/Code/bin/util_ManualOverlapGrid.py
index 2ccd2e68a..d09ab4168 100755
--- a/MonteCarloMarginalizeCode/Code/bin/util_ManualOverlapGrid.py
+++ b/MonteCarloMarginalizeCode/Code/bin/util_ManualOverlapGrid.py
@@ -246,7 +246,14 @@
 parser.add_argument("--extra-verbose", action="store_true",default=False, help="Lots of messages")
 parser.add_argument("--save-plots",default=False,action='store_true', help="Write plots to file (only useful for OSX, where interactive is default")
 parser.add_argument('--force-scatter',default=False,action='store_true',help='For hyperbolic analyses forces only scatter grid points.')
+parser.add_argument('--force-plunge',default=False,action='store_true',help='For hyperbolic analyses forces only plunge grid points.')
+parser.add_argument('--force-zoomwhirl',default=False,action='store_true',help='For hyperbolic analyses forces only zoomwhirl grid points.')
 opts=  parser.parse_args()
+
+force_options = [opts.force_scatter, opts.force_plunge, opts.force_zoomwhirl]  # Add more if needed
+if sum(bool(x) for x in force_options) > 1:
+    parser.error("CANNOT use multiple --force-X options at the same time!")
+
 if opts.inj_file_out:
     opts.fname = opts.inj_file_out.replace(".xml.gz","")
 
@@ -379,6 +386,30 @@ def evaluate_overlap_on_grid(hfbase,param_names, grid):
                     include_item = True
                 else:
                     include_item = False
+        
+        if opts.force_plunge:
+            if include_item==False:
+                # no need to evaluate if the point is already downselected out
+                pass
+            else:                
+                # removes non-plunge points from the hyperbolic grid
+                hypclass = Pgrid.extract_param('hypclass')
+                if hypclass == 'plunge':
+                    include_item = True
+                else:
+                    include_item = False
+                    
+        if opts.force_zoomwhirl:
+            if include_item==False:
+                # no need to evaluate if the point is already downselected out
+                pass
+            else:                
+                # removes non-zoomwhirl points from the hyperbolic grid
+                hypclass = Pgrid.extract_param('hypclass')
+                if hypclass == 'zoomwhirl':
+                    include_item = True
+                else:
+                    include_item = False
                 
         if include_item:
          grid_revised.append(line)
diff --git a/MonteCarloMarginalizeCode/Code/bin/util_ParameterPuffball.py b/MonteCarloMarginalizeCode/Code/bin/util_ParameterPuffball.py
index d169b32f2..cb22b9a80 100755
--- a/MonteCarloMarginalizeCode/Code/bin/util_ParameterPuffball.py
+++ b/MonteCarloMarginalizeCode/Code/bin/util_ParameterPuffball.py
@@ -55,8 +55,14 @@
 parser.add_argument("--enforce-duration-bound",default=None,type=float,help="If present, enforce a duration bound. Used to prevent grid placement for obscenely long signals, when the window size is prescribed")
 parser.add_argument("--regularize",action='store_true',help="Add some ad-hoc terms based on priors, to help with nearly-singular matricies")
 parser.add_argument('--force-scatter',default=False,action='store_true',help='For hyperbolic analyses forces only scatter grid points.')
+parser.add_argument('--force-plunge',default=False,action='store_true',help='For hyperbolic analyses forces only plunge grid points.')
+parser.add_argument('--force-zoomwhirl',default=False,action='store_true',help='For hyperbolic analyses forces only zoomwhirl grid points.')
 opts=  parser.parse_args()
 
+force_options = [opts.force_scatter, opts.force_plunge, opts.force_zoomwhirl]  # Add more if needed
+if sum(bool(x) for x in force_options) > 1:
+    parser.error("CANNOT use multiple --force-X options at the same time!")
+
 if opts.random_parameter is None:
     opts.random_parameter = []
 
@@ -275,6 +281,30 @@ def test_index_distance(indx,thresh=opts.force_away):
             else:
                 include_item = False
     
+    if opts.force_plunge:
+        if include_item==False:
+            # no need to evaluate if the point is already downselected out
+            pass            
+        else:        
+            # removes non-plunge points from the hyperbolic grid
+            hypclass = P.extract_param('hypclass')
+            if hypclass == 'plunge':
+                    include_item = True
+            else:
+                include_item = False
+                
+    if opts.force_zoomwhirl:
+        if include_item==False:
+            # no need to evaluate if the point is already downselected out
+            pass            
+        else:        
+            # removes non-zoomwhirl points from the hyperbolic grid
+            hypclass = P.extract_param('hypclass')
+            if hypclass == 'zoomwhirl':
+                    include_item = True
+            else:
+                include_item = False
+    
     if include_item:
         P_out.append(P)
 
diff --git a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py
index 490980fd7..44168db28 100755
--- a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py
+++ b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py
@@ -291,12 +291,19 @@ def unsafe_parse_arg_string_dict(my_argstr):
 parser.add_argument("--first-iteration-jumpstart",action='store_true',help="No ILE jobs the first iteration.  Assumes you already have .composite files and want to get going. Particularly helpful for subdag systems")
 parser.add_argument("--use-mtot-coords",action='store_true',help="Passed to the helper to configure CIP and PUFF for mtot instead of mc.")
 parser.add_argument("--force-scatter-grids",action='store_true',help="Eliminates all non-scatter intrinsic points from hyperbolic grids throughout the workflow.")
+parser.add_argument("--force-plunge-grids",action='store_true',help="Eliminates all non-plunge intrinsic points from hyperbolic grids throughout the workflow.")
+parser.add_argument("--force-zoomwhirl-grids",action='store_true',help="Eliminates all non-zoomwhirl intrinsic points from hyperbolic grids throughout the workflow.")
 parser.add_argument("--force-hyperbolic-22", action='store_true', help='Forces just the 22 modes for hyperbolic waveforms')
 opts=  parser.parse_args()
 
-# need --assume-hyperbolic when using --force-scatter-grids
-if opts.force_scatter_grids and not opts.assume_hyperbolic:
-    parser.error("--force-scatter-grids requires --assume-hyperbolic!")
+# Ensure --assume-hyperbolic is set when using any --force-X-grids option
+# Ensure only ONE of the --force-X-grids options is set
+force_grids = [opts.force_scatter_grids, opts.force_plunge_grids, opts.force_zoomwhirl_grids]
+if any(force_grids) and not opts.assume_hyperbolic:
+    parser.error("Using --force-scatter-grids, --force-plunge-grids, or --force-zoomwhirl-grids requires --assume-hyperbolic!")
+
+if sum(bool(x) for x in force_grids) > 1:
+    parser.error("CANNOT use multiple --force-X-grids options at the same time!")
 
 
 if (opts.use_ini):
@@ -668,6 +675,12 @@ def unsafe_parse_arg_string_dict(my_argstr):
     if opts.force_scatter_grids:
         cmd += " --force-scatter-grids "
         
+    if opts.force_plunge_grids:
+        cmd += " --force-plunge-grids "
+        
+    if opts.force_zoomwhirl_grids:
+        cmd += " --force-zoomwhirl-grids "
+        
     if opts.force_hyperbolic_22:
         cmd += " --force-hyperbolic-22 "
 if opts.assume_highq:
@@ -1122,6 +1135,12 @@ def unsafe_parse_arg_string_dict(my_argstr):
             
         if opts.force_scatter_grids:
             line += " --force-scatter "
+            
+        if opts.force_plunge_grids:
+            line += " --force-plunge "
+            
+        if opts.force_zoomwhirl_grids:
+            line += " --force-zoomwhirl "
         
         
     if not(opts.manual_extra_cip_args is None):
@@ -1201,6 +1220,12 @@ def unsafe_parse_arg_string_dict(my_argstr):
         if opts.force_scatter_grids:
             puff_params += ' --force-scatter '
             
+        if opts.force_plunge_grids:
+            puff_params += ' --force-plunge '
+            
+        if opts.force_zoomwhirl_grids:
+            puff_params += ' --force-zoomwhirl '
+            
 if opts.assume_highq:
     puff_params = puff_params.replace(' delta_mc ', ' eta ')  # use natural coordinates in the high q strategy. May want to do this always
     puff_max_it +=3

From 6f61267f80a8d76740e88833cfb96d46c243a0e6 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Tue, 8 Apr 2025 06:22:57 -0700
Subject: [PATCH 46/50] removed mc_range_str_cip from use_mtot_coords logic

---
 MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py b/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
index c53c36b15..36ed9bd5e 100755
--- a/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
+++ b/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py
@@ -1477,7 +1477,7 @@ def lambda_m_estimate(m):
         if not(opts.use_mtot_coords):
             helper_cip_args += mc_range_str_cip + eta_range_str_cip
         else:
-            helper_cip_args += mtot_range_str_cip + eta_range_str_cip + mc_range_str_cip
+            helper_cip_args += mtot_range_str_cip + eta_range_str_cip
     if opts.force_chi_max:
         helper_cip_args += " --chi-max {} ".format(opts.force_chi_max)
     if opts.force_chi_small_max:

From dff4f65e16480c6bfe57202ccb3b2a44f94c3ed2 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Mon, 21 Apr 2025 13:25:07 -0700
Subject: [PATCH 47/50] added E0 and p_phi0

---
 .../Code/bin/convert_output_format_inference2ile              | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/MonteCarloMarginalizeCode/Code/bin/convert_output_format_inference2ile b/MonteCarloMarginalizeCode/Code/bin/convert_output_format_inference2ile
index f84696cbe..077fe5904 100755
--- a/MonteCarloMarginalizeCode/Code/bin/convert_output_format_inference2ile
+++ b/MonteCarloMarginalizeCode/Code/bin/convert_output_format_inference2ile
@@ -84,6 +84,10 @@ for indx in np.arange(opts.target_size):
     P = lalsimutils.ChooseWaveformParams(m1=m1,m2=m2,dist=d)
     if "eccentricity" in samples_in.dtype.names:
      P.eccentricity = samples_in["eccentricity"][fac_reduce*indx]
+    if "E0" in samples_in.dtype.names:
+     P.E0 = samples_in["E0"][fac_reduce*indx]
+    if "p_phi0" in samples_in.dtype.names:
+     P.p_phi0 = samples_in["p_phi0"][fac_reduce*indx]
     if "time_maxl" in samples_in.dtype.names:
      P.tref = samples_in["time_maxl"][fac_reduce*indx]
     elif 'time' in samples_in.dtype.names:

From 7dbd9671ae18d68024648019df0a4454ca60a832 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Fri, 4 Jul 2025 05:21:45 -0700
Subject: [PATCH 48/50] impact parameter and scattering angle

---
 .../Code/RIFT/lalsimutils.py                  | 70 ++++++++++++++++++-
 .../Code/RIFT/misc/samples_utils.py           | 42 +++++++++++
 2 files changed, 111 insertions(+), 1 deletion(-)

diff --git a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
index 93c3ecfba..e7aa2cd3f 100644
--- a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
+++ b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
@@ -332,7 +332,7 @@ def lsu_StringFromPNOrder(order):
 #
 # Class to hold arguments of ChooseWaveform functions
 #
-valid_params = ['m1', 'm2', 's1x', 's1y', 's1z', 's2x', 's2y', 's2z', 'chi1_perp', 'chi2_perp', 'chi1_perp_bar', 'chi2_perp_bar','chi1_perp_u', 'chi2_perp_u', 's1z_bar', 's2z_bar', 'lambda1', 'lambda2', 'theta','phi', 'phiref',  'psi', 'incl', 'tref', 'dist', 'mc', 'mc_ecc', 'eta', 'delta_mc', 'chi1', 'chi2', 'thetaJN', 'phiJL', 'theta1', 'theta2', 'cos_theta1', 'cos_theta2',  'theta1_Jfix', 'theta2_Jfix', 'psiJ', 'beta', 'cos_beta', 'sin_phiJL', 'cos_phiJL', 'phi12', 'phi1', 'phi2', 'LambdaTilde', 'DeltaLambdaTilde', 'lambda_plus', 'lambda_minus', 'q', 'mtot','xi','chiz_plus', 'chiz_minus', 'chieff_aligned','fmin','fref', "SOverM2_perp", "SOverM2_L", "DeltaOverM2_perp", "DeltaOverM2_L", "shu","ampO", "phaseO",'eccentricity','chi_pavg','mu1','mu2','eos_table_index', 'E0', 'p_phi0', 'hypclass']
+valid_params = ['m1', 'm2', 's1x', 's1y', 's1z', 's2x', 's2y', 's2z', 'chi1_perp', 'chi2_perp', 'chi1_perp_bar', 'chi2_perp_bar','chi1_perp_u', 'chi2_perp_u', 's1z_bar', 's2z_bar', 'lambda1', 'lambda2', 'theta','phi', 'phiref',  'psi', 'incl', 'tref', 'dist', 'mc', 'mc_ecc', 'eta', 'delta_mc', 'chi1', 'chi2', 'thetaJN', 'phiJL', 'theta1', 'theta2', 'cos_theta1', 'cos_theta2',  'theta1_Jfix', 'theta2_Jfix', 'psiJ', 'beta', 'cos_beta', 'sin_phiJL', 'cos_phiJL', 'phi12', 'phi1', 'phi2', 'LambdaTilde', 'DeltaLambdaTilde', 'lambda_plus', 'lambda_minus', 'q', 'mtot','xi','chiz_plus', 'chiz_minus', 'chieff_aligned','fmin','fref', "SOverM2_perp", "SOverM2_L", "DeltaOverM2_perp", "DeltaOverM2_L", "shu","ampO", "phaseO",'eccentricity','chi_pavg','mu1','mu2','eos_table_index', 'E0', 'p_phi0', 'hypclass', 'b_hyp', 'phi_scattter']
 
 tex_dictionary  = {
  "mtot": '$M$',
@@ -1006,6 +1006,74 @@ def extract_param(self,p):
                 else:
                     print('No peaks detected after reclassifcation')
                     return 'meaningless'
+                
+        if p == 'b_hyp':
+            # Calculates impact parameter for hyperbolic cases. 
+            # See https://arxiv.org/pdf/2309.07228 for details.
+            # check if valid use
+            if self.E0 == 0.0:
+                print('Invalid use of b_hyp: non-hyperbolic configuration')
+                return None
+                        
+            h = self.E0 # mass normalized local system energy at initial separation
+            j = self.p_phi0 # mass normalized local system angular momentum at initial separation
+            
+            M1=self.m1/lal.MSUN_SI
+            M2=self.m2/lal.MSUN_SI
+            nu=M1*M2/((M1+M2)**2)
+            
+            E_eff = 1. + (h**2 - 1)/(2*nu) # effective energy of the one-body particle            
+            b = j*h / (np.sqrt(E_eff**2 - 1)) # dimensionless impact parameter normalized by total mass
+            
+            return b
+        
+        if p == 'phi_scatter':
+            # Calculates scattering angle for hyperbolic scatters 
+            # See https://arxiv.org/pdf/1402.7307,
+            # https://teobresums.bitbucket.io/gallery/scattering-angles/
+            
+            # NOTE that you can calculate this for capture/plunge cases, but it is physically meaningless
+            
+            # check if valid use
+            if self.E0 == 0.0:
+                print('Invalid use of phi_scatter: non-hyperbolic configuration')
+                return None
+            
+            # Generate waveform
+            pars = {
+                'M'                  : (self.m1+self.m2)/lal.MSUN_SI,
+                'q'                  : self.m1/self.m2,
+                'H_hyp'              : self.E0, # energy at initial separation
+                'j_hyp'              : self.p_phi0, # angular momentum at initial separation
+                'r_hyp'              : 6000.0,
+                'LambdaAl2'            : self.lambda1,
+                'LambdaBl2'            : self.lambda2,
+                'chi1'              : self.s1z, #note that there are no transverse spins
+                'chi2'              : self.s2z,
+                'dt'                : self.deltaT,
+                'domain'             : 0, # 0 sets time domain
+                'arg_out'            : 1, # Request multipoles and dynamics as output of the function call - 1=yes
+                'nqc'                : 2, # sets the NQCs, 2=no
+                'nqc_coefs_hlm'      : 0, # Option for the NQC model used in the waveform. 0=none
+                'nqc_coefs_flx'      : 0, # Option for the NQC model used in the flux. 0=none
+                'use_mode_lm'        : [1], # 22 mode
+                'output_lm'          : [1],
+                'srate_interp'       : 1./self.deltaT,
+                'use_geometric_units': 0,
+                'interp_uniform_grid': 1,
+                'initial_frequency'  : self.fmin,
+                'ode_tmax'           : 3e4,
+                'distance'           : self.dist/(lal.PC_SI*1e6),
+                'inclination'        : self.incl,
+                'output_hpc'         : 0 # output plus and cross polarizations, 0=no
+            }
+            
+            t, hptmp, hctmp, hlmtmp, dyn = EOBRun_module.EOBRunPy(pars)
+            
+            angle = dyn['phi'][-1] - dyn['phi'][0] - np.pi
+            
+            return angle
+            
         if p == 'delta' or p=='delta_mc':  # Same access routine
             return (self.m1-self.m2)/(self.m1+self.m2)
         if p == 'mc':
diff --git a/MonteCarloMarginalizeCode/Code/RIFT/misc/samples_utils.py b/MonteCarloMarginalizeCode/Code/RIFT/misc/samples_utils.py
index 4d96e57f2..064285311 100644
--- a/MonteCarloMarginalizeCode/Code/RIFT/misc/samples_utils.py
+++ b/MonteCarloMarginalizeCode/Code/RIFT/misc/samples_utils.py
@@ -137,6 +137,18 @@ def extract_combination_from_LI(samples_LI, p):
         a = a[8:]
         if a =='q' and 'm1' in samples_LI.dtype.names:
             return samples_LI["m1"]/samples_LI["m2"]
+        
+    if p == 'b_hyp':
+        samples = np.array([samples_LI["m1"], samples_LI["m2"], samples_LI["a1x"], samples_LI["a1y"], samples_LI["a1z"], samples_LI["a2x"], samples_LI["a2y"], samples_LI["a2z"], samples_LI["E0"], samples_LI["p_phi0"]]).T
+        with Pool(12) as pool:   
+            b_hyp = np.array(pool.map(f_b_hyp, samples))          
+        return b_hyp
+    
+    if p == 'phi_scatter':
+        samples = np.array([samples_LI["m1"], samples_LI["m2"], samples_LI["a1x"], samples_LI["a1y"], samples_LI["a1z"], samples_LI["a2x"], samples_LI["a2y"], samples_LI["a2z"], samples_LI["E0"], samples_LI["p_phi0"]]).T
+        with Pool(12) as pool:   
+            phi_scatter = np.array(pool.map(f_phi_scatter, samples))          
+        return phi_scatter
 
     print(" No access for parameter ", p)
     return np.zeros(len(samples_LI['m1']))  # to avoid causing a hard failure
@@ -272,3 +284,33 @@ def fchip(sample):
             P.s2z = sample[7]
             chip = P.extract_param('chi_p')
             return chip  
+
+def f_b_hyp(sample):
+            P=lalsimutils.ChooseWaveformParams()
+            P.m1 = sample[0]
+            P.m2 = sample[1]
+            P.s1x = sample[2]
+            P.s1y = sample[3]
+            P.s1z = sample[4]
+            P.s2x = sample[5]
+            P.s2y = sample[6]
+            P.s2z = sample[7]
+            P.E0 = sample[8]
+            P.p_phi0 = sample[9]
+            b_hyp = P.extract_param('b_hyp')
+            return b_hyp
+        
+def f_phi_scatter(sample):
+            P=lalsimutils.ChooseWaveformParams()
+            P.m1 = sample[0]
+            P.m2 = sample[1]
+            P.s1x = sample[2]
+            P.s1y = sample[3]
+            P.s1z = sample[4]
+            P.s2x = sample[5]
+            P.s2y = sample[6]
+            P.s2z = sample[7]
+            P.E0 = sample[8]
+            P.p_phi0 = sample[9]
+            phi_scatter = P.extract_param('phi_scatter')
+            return phi_scatter
\ No newline at end of file

From 802a86d7a79567df0f2c7305a00c1d35e046dab8 Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Fri, 4 Jul 2025 05:50:32 -0700
Subject: [PATCH 49/50] fixed latex translation for b_hyp, phi_scatter

---
 MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
index e7aa2cd3f..b7a92aff4 100644
--- a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
+++ b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py
@@ -375,6 +375,8 @@ def lsu_StringFromPNOrder(order):
   "eccentricity":"$e$",
   "E0": "$E_0 / M$",
   "p_phi0": r"$p_{\phi}^0$",
+  "b_hyp" : r"$b_{hyp}$",
+  "phi_scatter" : r"$\phi_{sc}$",  
   # tex labels for inherited LI names
  "a1z": r'$\chi_{1,z}$',
  "a2z": r'$\chi_{2,z}$',

From b4682cf0af36b11ff32875d03690f1903d0c8d8a Mon Sep 17 00:00:00 2001
From: Chad Henshaw <cgh3@gatech.edu>
Date: Fri, 4 Jul 2025 07:17:44 -0700
Subject: [PATCH 50/50] added b_hyp

---
 .../Code/bin/util_CharacterizePosterior.py                 | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/MonteCarloMarginalizeCode/Code/bin/util_CharacterizePosterior.py b/MonteCarloMarginalizeCode/Code/bin/util_CharacterizePosterior.py
index 78ae72b1b..f73681e8e 100755
--- a/MonteCarloMarginalizeCode/Code/bin/util_CharacterizePosterior.py
+++ b/MonteCarloMarginalizeCode/Code/bin/util_CharacterizePosterior.py
@@ -36,6 +36,13 @@
         dat_1d = (dat['m2']*dat['a2z']+dat['m1']*dat['a1z'])/(dat['m1']+dat['m2'])
     elif param=='mtot':
         dat_1d = dat['m1']+dat['m2']
+    elif param=='b_hyp':
+        h = dat['E0']
+        j = dat['p_phi0']
+        nu = dat['m1']*dat['m2']/((dat['m1']+dat['m2'])**2)
+        E_eff = 1. + (h**2 - 1)/(2*nu)
+        b = j*h / (np.sqrt(E_eff**2 - 1))
+        dat_1d = b
     else:
         dat_1d = dat[param]
     quant_here  = np.percentile(dat_1d,100*quantile_list)