Merge branch 'lalinference-point-release' into 'o3-release'

Lalinference point release

See merge request !840

configure.ac

 AC_PREREQ([2.63])
-AC_INIT([LALSuite],[6.57],[lal-discuss@ligo.org])
+AC_INIT([LALSuite],[6.58],[lal-discuss@ligo.org])
 AC_CONFIG_SRCDIR([configure.ac])
 AC_CONFIG_AUX_DIR([gnuscripts])
 AC_CONFIG_MACRO_DIR([gnuscripts])

lalinference/configure.ac

 AC_PREREQ([2.63])
-AC_INIT([LALInference],[1.11.0.1],[lal-discuss@ligo.org])
+AC_INIT([LALInference],[1.11.1],[lal-discuss@ligo.org])
 AC_CONFIG_HEADERS([src/config.h src/LALInferenceConfig.h])
 AC_CONFIG_SRCDIR([src/LALInference.h])
 AC_CONFIG_AUX_DIR([gnuscripts])
@@ -83,9 +83,9 @@ LALSUITE_DISTCHECK_CONFIGURE_FLAGS
 #    then increment age.
 # 6. if any interfaces have been removed since the last public release,
 #    then set age to 0.
-AC_SUBST([LIBCURRENT],[17])
+AC_SUBST([LIBCURRENT],[18])
 AC_SUBST([LIBREVISION],[0])
-AC_SUBST([LIBAGE],[0])
+AC_SUBST([LIBAGE],[1])
 AC_SUBST([LIBVERSION],[${LIBCURRENT}:${LIBREVISION}:${LIBAGE}])
 
 # nightly build

lalinference/debian/changelog

+lalinference (1.11.1-1) unstable; urgency=low
+
+  * O3 point release
+
+ -- Adam Mercer <adam.mercer@ligo.org>  Fri, 21 Jun 2019 11:06:11 -0700
+
 lalinference (1.11.0-1) unstable; urgency=low
 
   * O3 point release

lalinference/lalinference.spec.in

@@ -225,6 +225,9 @@ rm -Rf ${RPM_BUILD_DIR}/%{name}-%{version}%{?nightly:-%{nightly}}
 
 # dates should be formatted using: 'date +"%a %b %d %Y"'
 %changelog
+* Fri Jun 21 2019 Adam Mercer <adam.mercer@ligo.org> 1.11.1-1
+- O3 point release
+
 * Tue May 21 2019 Adam Mercer <adam.mercer@ligo.org> 1.11.0-1
 - O3 point release
 

lalinference/python/cbcBayesPosToSimInspiral.py

@@ -102,7 +102,7 @@ def compute_mass_parameterizations(samples):
         mc = samples['mchirp']
         if not has_eta:
             if has_q:
-                eta = bppu.q2eta(mc, samples['q'])
+                eta = bppu.q2eta(samples['q'])
             else:
                 raise ValueError("Chirp mass given with no mass ratio.")
         else:

lalinference/python/lalinference/bayespputils.py

@@ -136,15 +136,15 @@ spinParamsAli=['spin1','spin2','a1z','a2z']
 spinParamsEff=['chi','effectivespin','chi_eff','chi_tot','chi_p']
 spinParams=spinParamsPrec+spinParamsEff+spinParamsAli
 # Source frame params
-cosmoParam=['m1_source','m2_source','mtotal_source','mc_source','redshift','mf_source']
+cosmoParam=['m1_source','m2_source','mtotal_source','mc_source','redshift','mf_source','m1_source_maxldist','m2_source_maxldist','mtotal_source_maxldist','mc_source_maxldist','redshift_maxldist','mf_source_maxldist']
 #Strong Field
 ppEParams=['ppEalpha','ppElowera','ppEupperA','ppEbeta','ppElowerb','ppEupperB','alphaPPE','aPPE','betaPPE','bPPE']
 tigerParams=['dchi%i'%(i) for i in range(8)] + ['dchi%il'%(i) for i in [5,6] ] + ['dxi%d'%(i+1) for i in range(6)] + ['dalpha%i'%(i+1) for i in range(5)] + ['dbeta%i'%(i+1) for i in range(3)] + ['dsigma%i'%(i+1) for i in range(4)]
 bransDickeParams=['omegaBD','ScalarCharge1','ScalarCharge2']
 massiveGravitonParams=['lambdaG']
-tidalParams=['lambda1','lambda2','lam_tilde','dlam_tilde','lambdat','dlambdat']
+tidalParams=['lambda1','lambda2','lam_tilde','dlam_tilde','lambdat','dlambdat','lambdas','bluni']
 eosParams=['logp1','gamma1','gamma2','gamma3']
-energyParams=['e_rad', 'l_peak']
+energyParams=['e_rad', 'l_peak','e_rad_maxldist']
 strongFieldParams=ppEParams+tigerParams+bransDickeParams+massiveGravitonParams+tidalParams+eosParams+energyParams
 
 #Extrinsic
@@ -163,7 +163,7 @@ calibParams=['calpha_l1','calpha_h1','calpha_v1','calamp_l1','calamp_h1','calamp
 ## Greedy bin sizes for cbcBPP and confidence leves used for the greedy bin intervals
 confidenceLevels=[0.67,0.9,0.95,0.99]
 
-greedyBinSizes={'mc':0.025,'m1':0.1,'m2':0.1,'mass1':0.1,'mass2':0.1,'mtotal':0.1,'mc_source':0.025,'m1_source':0.1,'m2_source':0.1,'mtotal_source':0.1,'eta':0.001,'q':0.01,'asym_massratio':0.01,'iota':0.01,'cosiota':0.02,'time':1e-4,'time_mean':1e-4,'distance':1.0,'dist':1.0,'redshift':0.01,'mchirp':0.025,'chirpmass':0.025,'spin1':0.04,'spin2':0.04,'a1z':0.04,'a2z':0.04,'a1':0.02,'a2':0.02,'phi1':0.05,'phi2':0.05,'theta1':0.05,'theta2':0.05,'ra':0.05,'dec':0.05,'chi':0.05,'chi_eff':0.05,'chi_tot':0.05,'chi_p':0.05,'costilt1':0.02,'costilt2':0.02,'thatas':0.05,'costheta_jn':0.02,'beta':0.05,'omega':0.05,'cosbeta':0.02,'ppealpha':1.0,'ppebeta':1.0,'ppelowera':0.01,'ppelowerb':0.01,'ppeuppera':0.01,'ppeupperb':0.01,'polarisation':0.04,'rightascension':0.05,'declination':0.05,'massratio':0.001,'inclination':0.01,'phase':0.05,'tilt1':0.05,'tilt2':0.05,'phi_jl':0.05,'theta_jn':0.05,'phi12':0.05,'flow':1.0,'phase_maxl':0.05,'calamp_l1':0.01,'calamp_h1':0.01,'calamp_v1':0.01,'calpha_h1':0.01,'calpha_l1':0.01,'calpha_v1':0.01,'logdistance':0.1,'psi':0.1,'costheta_jn':0.1,'mf':0.1,'mf_source':0.1,'af':0.02,'e_rad':0.1,'l_peak':0.1}
+greedyBinSizes={'mc':0.025,'m1':0.1,'m2':0.1,'mass1':0.1,'mass2':0.1,'mtotal':0.1,'mc_source':0.025,'m1_source':0.1,'m2_source':0.1,'mtotal_source':0.1,'mc_source_maxldist':0.025,'m1_source_maxldist':0.1,'m2_source_maxldist':0.1,'mtotal_source_maxldist':0.1,'eta':0.001,'q':0.01,'asym_massratio':0.01,'iota':0.01,'cosiota':0.02,'time':1e-4,'time_mean':1e-4,'distance':1.0,'dist':1.0,'distance_maxl':1.0,'redshift':0.01,'redshift_maxldist':0.01,'mchirp':0.025,'chirpmass':0.025,'spin1':0.04,'spin2':0.04,'a1z':0.04,'a2z':0.04,'a1':0.02,'a2':0.02,'phi1':0.05,'phi2':0.05,'theta1':0.05,'theta2':0.05,'ra':0.05,'dec':0.05,'chi':0.05,'chi_eff':0.05,'chi_tot':0.05,'chi_p':0.05,'costilt1':0.02,'costilt2':0.02,'thatas':0.05,'costheta_jn':0.02,'beta':0.05,'omega':0.05,'cosbeta':0.02,'ppealpha':1.0,'ppebeta':1.0,'ppelowera':0.01,'ppelowerb':0.01,'ppeuppera':0.01,'ppeupperb':0.01,'polarisation':0.04,'rightascension':0.05,'declination':0.05,'massratio':0.001,'inclination':0.01,'phase':0.05,'tilt1':0.05,'tilt2':0.05,'phi_jl':0.05,'theta_jn':0.05,'phi12':0.05,'flow':1.0,'phase_maxl':0.05,'calamp_l1':0.01,'calamp_h1':0.01,'calamp_v1':0.01,'calpha_h1':0.01,'calpha_l1':0.01,'calpha_v1':0.01,'logdistance':0.1,'psi':0.1,'costheta_jn':0.1,'mf':0.1,'mf_source':0.1,'mf_source_maxldist':0.1,'af':0.02,'e_rad':0.1,'e_rad_maxldist':0.1,'l_peak':0.1}
 for s in snrParams:
     greedyBinSizes[s]=0.05
 for derived_time in ['h1_end_time','l1_end_time','v1_end_time','h1l1_delay','l1v1_delay','h1v1_delay']:
@@ -293,10 +293,17 @@ def get_prior(name):
       'mtotal_source':None,
       'mc_source':None,
       'redshift':None,
+      'm1_source_maxldist':None,
+      'm2_source_maxldist':None,
+      'mtotal_source_maxldist':None,
+      'mc_source_maxldist':None,
+      'redshift_maxldist':None,
       'mf':None,
       'mf_source':None,
+      'mf_source_maxldist':None,
       'af':None,
       'e_rad':None,
+      'e_rad_maxldist':None,
       'l_peak':None,
       'spin1':'uniform',
       'spin2':'uniform',
@@ -320,6 +327,7 @@ def get_prior(name):
       'time':'uniform',
       'time_mean':'uniform',
       'dist':'x**2',
+      'distance_maxl':'x**2',
       'ra':'uniform',
       'dec':'np.cos',
       'phase':'uniform',
@@ -343,6 +351,8 @@ def get_prior(name):
       'lambda2': 'uniform',
       'lam_tilde' : None,
       'dlam_tilde': None,
+      'lambdas':'uniform',
+      'bluni':'uniform',
       'logp1':None,
       'gamma1':None,
       'gamma2':None,
@@ -391,10 +401,17 @@ def plot_label(param):
         'mtotal_source':r'$M_\mathrm{total}^\mathrm{source}\,(\mathrm{M}_\odot)$',
         'mc_source':r'$\mathcal{M}^\mathrm{source}\,(\mathrm{M}_\odot)$',
         'redshift':r'$z$',
+        'm1_source_maxldist':r'$m_{1}^\mathrm{source - maxLdist}\,(\mathrm{M}_\odot)$',
+        'm2_source_maxldist':r'$m_{2}^\mathrm{source - maxLdist}\,(\mathrm{M}_\odot)$',
+        'mtotal_source_maxldist':r'$M_\mathrm{total}^\mathrm{source - maxLdist}\,(\mathrm{M}_\odot)$',
+        'mc_source_maxldist':r'$\mathcal{M}^\mathrm{source - maxLdist}\,(\mathrm{M}_\odot)$',
+        'redshift_maxldist':r'$z^\mathrm{maxLdist}$',
         'mf':r'$M_\mathrm{final}\,(\mathrm{M}_\odot)$',
         'mf_source':r'$M_\mathrm{final}^\mathrm{source}\,(\mathrm{M}_\odot)$',
+        'mf_source_maxldist':r'$M_\mathrm{final}^\mathrm{source - maxLdist}\,(\mathrm{M}_\odot)$',
         'af':r'$a_\mathrm{final}$',
         'e_rad':r'$E_\mathrm{rad}\,(\mathrm{M}_\odot)$',
+        'e_rad_maxldist':r'$E_\mathrm{rad}^\mathrm{maxLdist}\,(\mathrm{M}_\odot)$',
         'l_peak':r'$L_\mathrm{peak}\,(10^{56}\,\mathrm{ergs}\,\mathrm{s}^{-1})$',
         'spin1':r'$S_1$',
         'spin2':r'$S_2$',
@@ -418,9 +435,11 @@ def plot_label(param):
         'time':r'$t_\mathrm{c}\,(\mathrm{s})$',
         'time_mean':r'$<t>\,(\mathrm{s})$',
         'dist':r'$d_\mathrm{L}\,(\mathrm{Mpc})$',
+        'distance_maxl':r'$d_\mathrm{L}^\mathrm{maxL}\,(\mathrm{Mpc})$',
         'ra':r'$\alpha$',
         'dec':r'$\delta$',
         'phase':r'$\phi\,(\mathrm{rad})$',
+        'phase_maxl':r'$\phi^\mathrm{maxL}\,(\mathrm{rad})$',
         'psi':r'$\psi\,(\mathrm{rad})$',
         'theta_jn':r'$\theta_\mathrm{JN}\,(\mathrm{rad})$',
         'costheta_jn':r'$\mathrm{cos}(\theta_\mathrm{JN})$',
@@ -484,7 +503,9 @@ def plot_label(param):
         'h1_optimal_snr':r'$\rho^{opt}_{H1}$',
         'l1_optimal_snr':r'$\rho^{opt}_{L1}$',
         'v1_optimal_snr':r'$\rho^{opt}_{V1}$',
-        'matched_filter_snr':r'$\rho^{MF}$'
+        'matched_filter_snr':r'$\rho^{MF}$',
+        'lambdas':r'$\Lambda_S$',
+        'bluni' : r'$BL_{uniform}$'
       }
 
     # Handle cases where multiple names have been used
@@ -917,7 +938,7 @@ class Posterior(object):
         ('mass2' not in pos.names or 'm2' not in pos.names):
 
             pos.append_mapping(('m1','m2'),q2ms,(mchirp_name,q_name))
-            pos.append_mapping('eta',q2eta,(mchirp_name,q_name))
+            pos.append_mapping('eta',q2eta,q_name)
 
         if ('m1' in pos.names and 'm2' in pos.names and not 'mtotal' in pos.names ):
             pos.append_mapping('mtotal', lambda m1,m2: m1+m2, ('m1','m2') )
@@ -995,31 +1016,6 @@ class Posterior(object):
             except KeyError:
                 print("Warning: Cannot find spin parameters.  Skipping spin angle calculations.")
 
-        #Calculate effective precessing spin magnitude
-        if ('a1' in pos.names and 'tilt1' in pos.names and 'm1' in pos.names ) and ('a2' in pos.names and 'tilt2' in pos.names and 'm2' in pos.names):
-            pos.append_mapping('chi_p', chi_precessing, ['a1', 'tilt1', 'm1', 'a2', 'tilt2', 'm2'])
-
-
-        # Calculate redshift from luminosity distance measurements
-        if('distance' in pos.names):
-            pos.append_mapping('redshift', calculate_redshift, 'distance')
-        elif('dist' in pos.names):
-            pos.append_mapping('redshift', calculate_redshift, 'dist')
-
-
-        # Calculate source mass parameters
-        if ('m1' in pos.names) and ('redshift' in pos.names):
-            pos.append_mapping('m1_source', source_mass, ['m1', 'redshift'])
-
-        if ('m2' in pos.names) and ('redshift' in pos.names):
-            pos.append_mapping('m2_source', source_mass, ['m2', 'redshift'])
-
-        if ('mtotal' in pos.names) and ('redshift' in pos.names):
-            pos.append_mapping('mtotal_source', source_mass, ['mtotal', 'redshift'])
-
-        if ('mc' in pos.names) and ('redshift' in pos.names):
-            pos.append_mapping('mc_source', source_mass, ['mc', 'redshift'])
-
         #Store signed spin magnitudes in separate parameters and make a1,a2 magnitudes
         if 'a1' in pos.names:
             if 'tilt1' in pos.names:
@@ -1060,6 +1056,11 @@ class Posterior(object):
         # Calculate redshift from luminosity distance measurements
         if('distance' in pos.names):
             pos.append_mapping('redshift', calculate_redshift, 'distance')
+        elif('dist' in pos.names):
+            pos.append_mapping('redshift', calculate_redshift, 'dist')
+        # If using the DistanceMarginalisation, compute the maxL redshift distribution from the maxL d_L
+        elif('distance_maxl' in pos.names):
+            pos.append_mapping('redshift_maxldist', calculate_redshift, 'distance_maxl')
 
         # Calculate source mass parameters
         if ('m1' in pos.names) and ('redshift' in pos.names):
@@ -1074,9 +1075,22 @@ class Posterior(object):
         if ('mc' in pos.names) and ('redshift' in pos.names):
             pos.append_mapping('mc_source', source_mass, ['mc', 'redshift'])
 
+        # Calculate source mass parameters if DistanceMarginalisation was used, using the maxL distance and redshift
+        if ('m1' in pos.names) and ('redshift_maxldist' in pos.names):
+            pos.append_mapping('m1_source_maxldist', source_mass, ['m1', 'redshift_maxldist'])
+
+        if ('m2' in pos.names) and ('redshift_maxldist' in pos.names):
+            pos.append_mapping('m2_source_maxldist', source_mass, ['m2', 'redshift_maxldist'])
+
+        if ('mtotal' in pos.names) and ('redshift_maxldist' in pos.names):
+            pos.append_mapping('mtotal_source_maxldist', source_mass, ['mtotal', 'redshift_maxldist'])
+
+        if ('mc' in pos.names) and ('redshift_maxldist' in pos.names):
+            pos.append_mapping('mc_source_maxldist', source_mass, ['mc', 'redshift_maxldist'])
+
         #Calculate new tidal parameters
         new_tidal_params = ['lam_tilde','dlam_tilde']
-        old_tidal_params = ['lambda1','lambda2','eta']
+        old_tidal_params = ['lambda1','lambda2','q']
         if 'lambda1' in pos.names or 'lambda2' in pos.names:
             try:
                 pos.append_mapping(new_tidal_params, symm_tidal_params, old_tidal_params)
@@ -1156,6 +1170,12 @@ class Posterior(object):
             except Exception as e:
                 print("Could not calculate final source frame mass. The error was: %s"%(str(e)))
 
+        if ('mf' in pos.names) and ('redshift_maxldist' in pos.names):
+            try:
+                pos.append_mapping('mf_source_maxldist', source_mass, ['mf', 'redshift_maxldist'])
+            except Exception as e:
+                print("Could not calculate final source frame mass. The error was: %s"%(str(e)))
+
         # Calculate radiated energy and peak luminosity
         if ('mtotal_source' in pos.names) and ('mf_source' in pos.names):
             try:
@@ -1163,6 +1183,12 @@ class Posterior(object):
             except Exception as e:
                 print("Could not calculate radiated energy. The error was: %s"%(str(e)))
 
+        if ('mtotal_source_maxldist' in pos.names) and ('mf_source_maxldist' in pos.names):
+            try:
+                pos.append_mapping('e_rad_maxldist', lambda mtot_s, mf_s: mtot_s-mf_s, ['mtotal_source_maxldist', 'mf_source_maxldist'])
+            except Exception as e:
+                print("Could not calculate radiated energy. The error was: %s"%(str(e)))
+
         if ('q' in pos.names) and ('a1z' in pos.names) and ('a2z' in pos.names):
             try:
                 pos.append_mapping('l_peak', bbh_aligned_Lpeak_6mode_SHXJDK, ['q', 'a1z', 'a2z'])
@@ -3567,14 +3593,13 @@ def q2ms(mc,q):
 #
 #
 
-def q2eta(mc,q):
+def q2eta(q):
     """
-    Utility function for converting mchirp,q to eta. The
+    Utility function for converting q to eta. The
     rvalue is eta.
     """
-    m1,m2 = q2ms(mc,q)
-    eta = m1*m2/( (m1+m2)*(m1+m2) )
-    return eta
+    eta = q/((1+q)*(1+q))
+    return np.clip(eta,0,0.25) # Explicitly cap eta at 0.25, in case it exceeds this slightly due to floating-point issues
 #
 #
 
@@ -3721,12 +3746,24 @@ def component_momentum(m, a, theta, phi):
 #
 #
 
-def symm_tidal_params(lambda1,lambda2,eta):
+def symm_tidal_params(lambda1,lambda2,q):
     """
-    Calculate best tidal parameters
+    Calculate best tidal parameters [Eqs. (5) and (6) in Wade et al. PRD 89, 103012 (2014)]
+    Requires q <= 1
     """
-    lam_tilde = (8./13.)*((1.+7.*eta-31.*eta*eta)*(lambda1+lambda2) + np.sqrt(1.-4.*eta)*(1.+9.*eta-11.*eta*eta)*(lambda1-lambda2))
-    dlam_tilde = (1./2.)*(np.sqrt(1.-4.*eta)*(1.-13272.*eta/1319.+8944.*eta*eta/1319.)*(lambda1+lambda2) + (1.-15910.*eta/1319.+32850.*eta*eta/1319.+3380.*eta*eta*eta/1319.)*(lambda1-lambda2))
+    lambdap = lambda1 + lambda2
+    lambdam = lambda1 - lambda2
+
+    # Check that q <= 1, as expected
+    if np.any(q > 1):
+        raise ValueError("q > 1, while this function requires q <= 1.")
+
+    dmbym = (1. - q)/(1. + q) # Equivalent to sqrt(1 - 4*eta) for q <= 1
+
+    eta = q2eta(q)
+
+    lam_tilde = (8./13.)*((1.+7.*eta-31.*eta*eta)*lambdap + dmbym*(1.+9.*eta-11.*eta*eta)*lambdam)
+    dlam_tilde = (1./2.)*(dmbym*(1.-13272.*eta/1319.+8944.*eta*eta/1319.)*lambdap + (1.-15910.*eta/1319.+32850.*eta*eta/1319.+3380.*eta*eta*eta/1319.)*lambdam)
     return lam_tilde, dlam_tilde
 
 def spin_angles(fref,mc,eta,incl,a1,theta1,phi1,a2=None,theta2=None,phi2=None):
@@ -5408,15 +5445,15 @@ def find_ndownsample(samples, nDownsample):
         splineParams=["spcal_npts", "spcal_active","constantcal_active"]
         for i in np.arange(25):
             for k in lal.cached_detector_by_prefix:
-                splineParams.append(k+'_spcal_freq_'+str(i))
-                splineParams.append(k+'_spcal_logfreq_'+str(i))
+                splineParams.append(k.lower()+'_spcal_freq_'+str(i))
+                splineParams.append(k.lower()+'_spcal_logfreq_'+str(i))
 
         nonParams = ["logpost", "post", "cycle", "timestamp", "snrh1", "snrl1", "snrv1",
                      "margtime","margtimephi","margtime","time_max","time_min",
                      "time_mean", "time_maxl","sky_frame","psdscaleflag","logdeltaf","flow","f_ref",
                      "lal_amporder","lal_pnorder","lal_approximant","tideo","spino","signalmodelflag",
                      "temperature","nifo","nlocaltemps","ntemps","randomseed","samplerate","segmentlength","segmentstart",
-                     "t0", "phase_maxl", "azimuth", "cosalpha", "lal_amporder"] + logParams + snrParams + splineParams
+                     "t0", "phase_maxl", "azimuth", "cosalpha", "lal_amporder", "bluni"] + logParams + snrParams + splineParams
         fixedParams = [p for p in samples.colnames if all(x==samples[p][0] for x in samples[p])]
         print("Fixed parameters: "+str(fixedParams))
         nonParams.extend(fixedParams)

lalinference/src/LALInference.c

@@ -4314,3 +4314,102 @@ void LALInferencePrintSplineCalibration(FILE *output, LALInferenceThreadState *t
         }
     }
 }
+
+void LALInferenceBinaryLove(LALInferenceVariables *vars, REAL8 *lambda1, REAL8 *lambda2){
+
+    /**
+    *  Implelentation of the parametrization from Chatziioannou, Haster, Zimmerman (2018)
+    *  https://doi.org/10.1103/PhysRevD.97.104036
+    */
+
+    REAL8 lambdaS = *(REAL8*) LALInferenceGetVariable(vars, "lambdaS");
+    REAL8 lambdaSm1o5 = pow(lambdaS,-1.0/5.0);
+    REAL8 lambdaSm2o5 = lambdaSm1o5*lambdaSm1o5;
+    REAL8 lambdaSm3o5 = lambdaSm2o5*lambdaSm1o5;
+    REAL8 lambdaSsqrt = sqrt(lambdaS);
+
+    REAL8 q=1.0;
+
+    if (LALInferenceCheckVariable(vars,"q")) {
+        q = *(REAL8 *)LALInferenceGetVariable(vars,"q");
+    }
+    else{
+        REAL8 m1 = *(REAL8 *)LALInferenceGetVariable(vars,"mass1");
+        REAL8 m2 = *(REAL8 *)LALInferenceGetVariable(vars,"mass2");
+        q = m2/m1;
+    }
+    REAL8 q2 = q*q;
+
+    /* Eqn 2 , incorporating the dependance on mass ratio */
+
+    REAL8 BL_n = 0.743;
+    REAL8 q_for_Fnofq = pow(q,10.0/(3.0-BL_n));
+    REAL8 Fnofq = (1.0-q_for_Fnofq)/(1.0+q_for_Fnofq);
+
+    /* bXX and cXX coefficients are given in Table I */
+
+    REAL8 b11 = -27.7408;
+    REAL8 b12 = 8.42358;
+    REAL8 b21 = 122.686;
+    REAL8 b22 = -19.7551;
+    REAL8 b31 = -175.496;
+    REAL8 b32 = 133.708;
+
+    REAL8 c11 = -25.5593;
+    REAL8 c12 = 5.58527;
+    REAL8 c21 = 92.0337;
+    REAL8 c22 = 26.8586;
+    REAL8 c31 = -70.247;
+    REAL8 c32 = -56.3076i;
+
+    /* Eqn 1, giving the lambdaA_fitOnly, not yet accounting for the uncertainty in the fit */
+
+    REAL8 numerator = 1.0 + (b11*q*lambdaSm1o5) + (b12*q2*lambdaSm1o5) + (b21*q*lambdaSm2o5) + (b22*q2*lambdaSm2o5) + (b31*q*lambdaSm3o5) + (b32*q2*lambdaSm3o5);
+
+    REAL8 denominator = 1.0 + (c11*q*lambdaSm1o5) + (c12*q2*lambdaSm1o5) + (c21*q*lambdaSm2o5) + (c22*q2*lambdaSm2o5) + (c31*q*lambdaSm3o5) + (c32*q2*lambdaSm3o5);
+
+    REAL8 lambdaA_fitOnly = Fnofq*lambdaS*numerator/denominator;
+
+    /* Eqn 6, correction on fit for lambdaA caused by uncertainty in the mean of the lambdaS residual fit, 
+    *         using coefficients mu_1, mu_2 and mu_3 from Table II */
+
+    REAL8 lambdaA_lambdaS_meanCorr = (137.1252739/(lambdaS*lambdaS)) - (32.8026613/lambdaS) + 0.5168637;
+
+    /* Eqn 8, correction on fit for lambdaA caused by uncertainty in the standard deviation of lambdaS residual fit,
+    *         using coefficients sigma_1, sigma_2 sigma_3 and  sigma_4 from Table II */
+
+    REAL8 lambdaA_lambdaS_stdCorr = (-0.0000739*lambdaS*lambdaSsqrt) + (0.0103778*lambdaS) + (0.4581717*lambdaSsqrt)  - 0.8341913;
+
+    /* Eqn 7, correction on fit for lambdaA caused by uncertainty in the mean of the q residual fit, 
+    *         using coefficients mu_4 and mu_5 from Table II */
+
+    REAL8 lambdaA_q_meanCorr = (-11.2765281*q2) + (14.9499544*q) - 4.6638851;
+
+    /* Eqn 9, correction on fit for lambdaA caused by uncertainty in the standard deviation of the q residual fit, 
+    *         using coefficients sigma_5, sigma_6 and sigma_7 from Table II */
+
+    REAL8 lambdaA_q_stdCorr = (-201.4323962*q2) + (273.9268276*q) - 71.2342246;
+
+    /* Eqn 4, averaging the corrections from the mean of the residual fits */
+
+    REAL8 lambdaA_meanCorr = (lambdaA_lambdaS_meanCorr + lambdaA_q_meanCorr)/2.0;
+
+    /* Eqn 5, averaging the corrections from the standard deviations of the residual fits */
+
+    REAL8 lambdaA_stdCorr = sqrt((lambdaA_lambdaS_stdCorr*lambdaA_lambdaS_stdCorr) + (lambdaA_q_stdCorr*lambdaA_q_stdCorr));
+
+    /* Draw a correction on the fit from a Gaussian distribution wiht width lambdaA_stdCorr 
+    *  this is done by sampling a inverse cdf through a U{0,1} variable called BLuni */
+
+    REAL8 BLuni = *(REAL8*) LALInferenceGetVariable(vars, "BLuni");
+
+    REAL8 lambdaA_scatter = gsl_cdf_gaussian_Pinv(BLuni, lambdaA_stdCorr);
+
+    /* Add the correction of the residual mean and the Gaussian scatter to the lambdaA_fitOnly value */
+
+    REAL8 lambdaA = lambdaA_fitOnly + (lambdaA_meanCorr + lambdaA_scatter);
+
+    *lambda1 = lambdaS - lambdaA;
+    *lambda2 = lambdaS + lambdaA;
+    return;
+}

lalinference/src/LALInference.h

@@ -86,6 +86,8 @@
 #include <gsl/gsl_randist.h>