Commit 2b919c51 authored by Salvatore Vitale's avatar Salvatore Vitale
Browse files

working on postproc

parent e7feb9d7
......@@ -804,6 +804,7 @@ class Posterior(object):
'end_time': lambda inj:float(inj.get_end()),
'phi0':lambda inj:inj.phi0,
'phi_orb': lambda inj: inj.coa_phase,
'phase': lambda inj: inj.coa_phase,
'dist':lambda inj:inj.distance,
'distance':lambda inj:inj.distance,
'ra':self._inj_longitude,
......@@ -1935,8 +1936,11 @@ class Posterior(object):
return maplong
else:
return inj.longitude
def _inj_spins(self, inj, frame='OrbitalL'):
from lalsimulation import SimInspiralTransformPrecessingWvf2PE
spins = {}
f_ref = self._injFref
......@@ -1945,16 +1949,90 @@ class Posterior(object):
else:
axis = lalsim.SimInspiralGetFrameAxisFromString(frame)
s1x=inj.spin1x
s1y=inj.spin1y
s1z=inj.spin1z
s2x=inj.spin2x
s2y=inj.spin2y
s2z=inj.spin2z
iota=inj.inclination
m1, m2 = inj.mass1, inj.mass2
mc, eta = inj.mchirp, inj.eta
#iota, s1x, s1y, s1z, s2x, s2y, s2z = \
# lalsim.SimInspiralInitialConditionsPrecessingApproxs(inj.inclination,
# inj.spin1x, inj.spin1y, inj.spin1z,
# inj.spin2x, inj.spin2y, inj.spin2z,
# m1*lal.MSUN_SI, m2*lal.MSUN_SI, f_ref, inj.coa_phase, axis)
m1, m2 = inj.mass1, inj.mass2
mc, eta = inj.mchirp, inj.eta
a1, theta1, phi1 = cart2sph(s1x, s1y, s1z)
a2, theta2, phi2 = cart2sph(s2x, s2y, s2z)
# Convert to radiation frame
spins = {'a1':a1, 'theta1':theta1, 'phi1':phi1,
'a2':a2, 'theta2':theta2, 'phi2':phi2,
'iota':iota}
# If spins are aligned, save the sign of the z-component
if inj.spin1x == inj.spin1y == inj.spin2x == inj.spin2y == 0.:
spins['a1z'] = inj.spin1z
spins['a2z'] = inj.spin2z
L = orbital_momentum(f_ref, m1,m2, iota)
S1 = np.hstack((s1x, s1y, s1z))
S2 = np.hstack((s2x, s2y, s2z))
zhat = np.array([0., 0., 1.])
aligned_comp_spin1 = array_dot(S1, zhat)
aligned_comp_spin2 = array_dot(S2, zhat)
chi = aligned_comp_spin1 + aligned_comp_spin2 + \
np.sqrt(1. - 4.*eta) * (aligned_comp_spin1 - aligned_comp_spin2)
S1 *= m1**2
S2 *= m2**2
J = L + S1 + S2
tilt1 = array_ang_sep(L, S1)
tilt2 = array_ang_sep(L, S2)
beta = array_ang_sep(J, L)
spins['beta'] = beta
spins['spinchi'] = chi
theta_jn,phi_jl,tilt1,tilt2,phi12,chi1,chi2=SimInspiralTransformPrecessingWvf2PE(inj.inclination,inj.spin1x, inj.spin1y, inj.spin1z,inj.spin2x, inj.spin2y, inj.spin2z, m1, m2, f_ref, inj.coa_phase)
print 'fref ',f_ref
#theta_jn,phi_jl,tilt1,tilt2,phi12,chi1,chi2=SimInspiralTransformPrecessingWvf2PE(iota, s1x, s1y, s1z, s2x, s2y, s2z,m1*lal.MSUN_SI, m2*lal.MSUN_SI, f_ref, inj.coa_phase)
print inj.inclination
print inj.spin1x, inj.spin1y, inj.spin1z
print inj.spin2x, inj.spin2y, inj.spin2z
print m1, m2, f_ref, inj.coa_phase
print "in polar" , theta_jn,phi_jl,tilt1,tilt2,phi12,chi1,chi2
spins['theta_jn']=theta_jn
spins['phi12']=phi12
spins['tilt1']=tilt1
spins['tilt2']=tilt2
spins['phi_jl']=phi_jl
## print trying to inverst
print "Inverting "
iota_back,a1x_back,a1y_back,a1z_back,a2x_back,a2y_back,a2z_back = \
lalsim.SimInspiralTransformPrecessingNewInitialConditions(theta_jn,phi_jl,tilt1,tilt2,phi12,chi1,chi2,m1*lal.MSUN_SI,m2*lal.MSUN_SI,f_ref,inj.coa_phase)
print a1x_back,a1y_back,a1z_back
print a2x_back,a2y_back,a2z_back
print iota_back
print "ODL WAY"
iota, s1x, s1y, s1z, s2x, s2y, s2z = \
lalsim.SimInspiralInitialConditionsPrecessingApproxs(inj.inclination,
inj.spin1x, inj.spin1y, inj.spin1z,
inj.spin2x, inj.spin2y, inj.spin2z,
m1*lal.MSUN_SI, m2*lal.MSUN_SI, f_ref, inj.coa_phase, axis)
theta_jn,phi_jl,tilt1,tilt2,phi12,chi1,chi2=SimInspiralTransformPrecessingWvf2PE(iota, s1x, s1y, s1z, s2x, s2y, s2z,m1, m2, f_ref, inj.coa_phase)
print theta_jn,phi_jl,tilt1,tilt2,phi12,chi1,chi2
# Convert to radiation frame
"""iota, s1x, s1y, s1z, s2x, s2y, s2z = \
lalsim.SimInspiralInitialConditionsPrecessingApproxs(inj.inclination,
inj.spin1x, inj.spin1y, inj.spin1z,
inj.spin2x, inj.spin2y, inj.spin2z,
m1*lal.MSUN_SI, m2*lal.MSUN_SI, f_ref, inj.coa_phase, axis)
a1, theta1, phi1 = cart2sph(s1x, s1y, s1z)
a2, theta2, phi2 = cart2sph(s2x, s2y, s2z)
......@@ -2013,7 +2091,9 @@ class Posterior(object):
# J should now be || z and L should have a azimuthal angle phi_jl
phi_jl = np.arctan2(L[1], L[0])
spins['phi_jl'] = np.pi + phi_jl
print "pure phijl from atan before manupulations" , phi_jl
print "-----> inj.coa_phase ",inj.coa_phase
spins['phi_jl'] = np.pi/2. + phi_jl
# bring L in the Z-X plane, with negative x
J = ROTATEZ(np.pi-phi_jl, J[0], J[1], J[2])
......@@ -2035,7 +2115,9 @@ class Posterior(object):
else:
phi12 = phi2 - phi1
spins['phi12'] = phi12
"""
return spins
class BurstPosterior(Posterior):
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment