zebra/pyzebra
3
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

apply black formating

Browse Source 
Co-authored-by: Ivan Usov <ivan.a.usov@gmail.com>
This commit is contained in:
zaharko 2020-04-23 10:46:00 +02:00
parent a287b78245
commit 6a3ca53de1
3 changed files with 231 additions and 177 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
pyzebra/app.py
View File

@ -77,6 +77,7 @@ def update_image():
image_glyph.color_mapper.low = im_min image_glyph.color_mapper.low = im_min
image_glyph.color_mapper.high = im_max image_glyph.color_mapper.high = im_max
def calculate_hkl(setup_type="nb_bi"): def calculate_hkl(setup_type="nb_bi"):
h = np.empty(shape=(IMAGE_H, IMAGE_W)) h = np.empty(shape=(IMAGE_H, IMAGE_W))
k = np.empty(shape=(IMAGE_H, IMAGE_W)) k = np.empty(shape=(IMAGE_H, IMAGE_W))
@ -435,6 +436,7 @@ def auto_toggle_callback(state):
update_image() update_image()
auto_toggle = Toggle(label="Auto Range", active=True, button_type="default") auto_toggle = Toggle(label="Auto Range", active=True, button_type="default")
auto_toggle.on_click(auto_toggle_callback) auto_toggle.on_click(auto_toggle_callback)

16
pyzebra/h5.py
View File

@ -1,5 +1,6 @@
import h5py import h5py
def read_h5meta(filepath): def read_h5meta(filepath):
"""Read and parse content of a h5meta file. """Read and parse content of a h5meta file.
@ -46,17 +47,18 @@ def read_detector_data(filepath):
det_data = {"data": data} det_data = {"data": data}
det_data["rot_angle"] = h5f["/entry1/area_detector2/rotation_angle"][:] # om, sometimes ph det_data["rot_angle"] = h5f["/entry1/area_detector2/rotation_angle"][:] # om, sometimes ph
det_data["pol_angle"] = h5f["/entry1/ZEBRA/area_detector2/polar_angle"][:] # gammad det_data["pol_angle"] = h5f["/entry1/ZEBRA/area_detector2/polar_angle"][:] # gammad
det_data["tlt_angle"] = h5f["/entry1/ZEBRA/area_detector2/tilt_angle"][:] # nud det_data["tlt_angle"] = h5f["/entry1/ZEBRA/area_detector2/tilt_angle"][:] # nud
det_data["ddist"] = h5f["/entry1/ZEBRA/area_detector2/distance"][:] det_data["ddist"] = h5f["/entry1/ZEBRA/area_detector2/distance"][:]
det_data["wave"] = h5f["/entry1/ZEBRA/monochromator/wavelength"][:] det_data["wave"] = h5f["/entry1/ZEBRA/monochromator/wavelength"][:]
det_data["chi_angle"] = h5f["/entry1/sample/chi"][:] # ch det_data["chi_angle"] = h5f["/entry1/sample/chi"][:] # ch
det_data["phi_angle"] = h5f["/entry1/sample/phi"][:] # ph det_data["phi_angle"] = h5f["/entry1/sample/phi"][:] # ph
det_data["UB"] = h5f["/entry1/sample/UB"][:].reshape(3,3) det_data["UB"] = h5f["/entry1/sample/UB"][:].reshape(3, 3)
return det_data return det_data
def open_h5meta(filepath): def open_h5meta(filepath):
"""Open h5meta file like *.cami """Open h5meta file like *.cami

288
pyzebra/xtal.py
View File

@ -4,7 +4,8 @@ from scipy.optimize import curve_fit
from matplotlib import pyplot as plt from matplotlib import pyplot as plt
import pyzebra import pyzebra
def z4frgn(wave,ga,nu):
def z4frgn(wave, ga, nu):
"""CALCULATES DIFFRACTION VECTOR IN LAB SYSTEM FROM GA AND NU """CALCULATES DIFFRACTION VECTOR IN LAB SYSTEM FROM GA AND NU
Args: Args:
@ -15,16 +16,17 @@ def z4frgn(wave,ga,nu):
""" """
sin = np.sin sin = np.sin
cos = np.cos cos = np.cos
pir = 180/np.pi pir = 180 / np.pi
gar = ga/pir gar = ga / pir
nur = nu/pir nur = nu / pir
z4 = [0., 0., 0.] z4 = [0.0, 0.0, 0.0]
z4[0]=( sin(gar)*cos(nur) )/wave z4[0] = (sin(gar) * cos(nur)) / wave
z4[1]=( cos(gar)*cos(nur)-1. )/wave z4[1] = (cos(gar) * cos(nur) - 1.0) / wave
z4[2]=( sin(nur) )/wave z4[2] = (sin(nur)) / wave
return z4 return z4
def phimat(phi): def phimat(phi):
"""BUSING AND LEVY CONVENTION ROTATION MATRIX FOR PHI OR OMEGA """BUSING AND LEVY CONVENTION ROTATION MATRIX FOR PHI OR OMEGA
@ -36,19 +38,20 @@ def phimat(phi):
""" """
sin = np.sin sin = np.sin
cos = np.cos cos = np.cos
pir = 180/np.pi pir = 180 / np.pi
phr = phi/pir phr = phi / pir
dum = np.zeros(9).reshape(3,3) dum = np.zeros(9).reshape(3, 3)
dum[0,0] = cos(phr) dum[0, 0] = cos(phr)
dum[0,1] = sin(phr) dum[0, 1] = sin(phr)
dum[1,0] = -dum[0,1] dum[1, 0] = -dum[0, 1]
dum[1,1] = dum[0,0] dum[1, 1] = dum[0, 0]
dum[2,2] = 1 dum[2, 2] = 1
return dum return dum
def z1frnb(wave,ga,nu,om):
def z1frnb(wave, ga, nu, om):
"""CALCULATE DIFFRACTION VECTOR Z1 FROM GA, OM, NU, ASSUMING CH=PH=0 """CALCULATE DIFFRACTION VECTOR Z1 FROM GA, OM, NU, ASSUMING CH=PH=0
Args: Args:
@ -58,13 +61,14 @@ def z1frnb(wave,ga,nu,om):
Z1 Z1
""" """
z4 = z4frgn(wave,ga,nu) z4 = z4frgn(wave, ga, nu)
dum = phimat(phi=om) dum = phimat(phi=om)
dumt = np.transpose(dum) dumt = np.transpose(dum)
z3 = dumt.dot(z4) z3 = dumt.dot(z4)
return z3 return z3
def chimat(chi): def chimat(chi):
"""BUSING AND LEVY CONVENTION ROTATION MATRIX FOR CHI """BUSING AND LEVY CONVENTION ROTATION MATRIX FOR CHI
@ -76,19 +80,20 @@ def chimat(chi):
""" """
sin = np.sin sin = np.sin
cos = np.cos cos = np.cos
pir = 180/np.pi pir = 180 / np.pi
chr = chi/pir chr = chi / pir
dum = np.zeros(9).reshape(3,3) dum = np.zeros(9).reshape(3, 3)
dum[0,0] = cos(chr) dum[0, 0] = cos(chr)
dum[0,2] = sin(chr) dum[0, 2] = sin(chr)
dum[1,1] = 1 dum[1, 1] = 1
dum[2,0] = -dum[0,2] dum[2, 0] = -dum[0, 2]
dum[2,2] = dum[0,0] dum[2, 2] = dum[0, 0]
return dum return dum
def z1frz3(z3,chi,phi):
def z1frz3(z3, chi, phi):
"""CALCULATE Z1 = [PHI]T.[CHI]T.Z3 """CALCULATE Z1 = [PHI]T.[CHI]T.Z3
Args: Args:
@ -108,7 +113,8 @@ def z1frz3(z3,chi,phi):
return z1 return z1
def z1frmd(wave,ga,om,chi,phi,nu):
def z1frmd(wave, ga, om, chi, phi, nu):
"""CALCULATE DIFFRACTION VECTOR Z1 FROM CH, PH, GA, OM, NU """CALCULATE DIFFRACTION VECTOR Z1 FROM CH, PH, GA, OM, NU
Args: Args:
@ -117,12 +123,13 @@ def z1frmd(wave,ga,om,chi,phi,nu):
Returns: Returns:
Z1 Z1
""" """
z3 = z1frnb(wave,ga,nu,om) z3 = z1frnb(wave, ga, nu, om)
z1 = z1frz3(z3,chi,phi) z1 = z1frz3(z3, chi, phi)
return z1 return z1
def det2pol(ddist,gammad,nud,x,y):
def det2pol(ddist, gammad, nud, x, y):
"""CONVERTS FROM DETECTOR COORDINATES TO POLAR COORDINATES """CONVERTS FROM DETECTOR COORDINATES TO POLAR COORDINATES
Args: Args:
@ -137,19 +144,20 @@ def det2pol(ddist,gammad,nud,x,y):
xpix = 0.734 xpix = 0.734
ypix = 1.4809 ypix = 1.4809
xobs = (x - xnorm)*xpix xobs = (x - xnorm) * xpix
yobs = (y - ynorm)*ypix yobs = (y - ynorm) * ypix
a = xobs a = xobs
b = ddist * np.cos(yobs/ddist) b = ddist * np.cos(yobs / ddist)
z = ddist * np.sin(yobs/ddist) z = ddist * np.sin(yobs / ddist)
d = np.sqrt(a*a+b*b) d = np.sqrt(a * a + b * b)
pir = 180/np.pi pir = 180 / np.pi
gamma = gammad + np.arctan2(a,b)*pir gamma = gammad + np.arctan2(a, b) * pir
nu = nud + np.arctan2(z,d)*pir nu = nud + np.arctan2(z, d) * pir
return gamma, nu return gamma, nu
def eqchph(z1): def eqchph(z1):
"""CALCULATE CHI, PHI TO PUT THE VECTOR Z1 IN THE EQUATORIAL PLANE """CALCULATE CHI, PHI TO PUT THE VECTOR Z1 IN THE EQUATORIAL PLANE
@ -159,13 +167,13 @@ def eqchph(z1):
Returns: Returns:
chi, phi chi, phi
""" """
pir = 180/np.pi pir = 180 / np.pi
if z1[0] != 0 or z1[1] != 0: if z1[0] != 0 or z1[1] != 0:
ph = np.arctan2(z1[1],z1[0]) ph = np.arctan2(z1[1], z1[0])
ph = ph * pir ph = ph * pir
d = np.sqrt(z1[0]*z1[0]+z1[1]*z1[1]) d = np.sqrt(z1[0] * z1[0] + z1[1] * z1[1])
ch = np.arctan2(z1[2],d) ch = np.arctan2(z1[2], d)
ch = ch * pir ch = ch * pir
else: else:
ph = 0 ph = 0
@ -179,7 +187,7 @@ def eqchph(z1):
return ch, ph return ch, ph
def dandth(wave,z1): def dandth(wave, z1):
"""CALCULATE D-SPACING (REAL SPACE) AND THETA FROM LENGTH OF Z """CALCULATE D-SPACING (REAL SPACE) AND THETA FROM LENGTH OF Z
Args: Args:
@ -188,16 +196,16 @@ def dandth(wave,z1):
Returns: Returns:
ds, th ds, th
""" """
pir = 180/np.pi pir = 180 / np.pi
ierr = 0 ierr = 0
dstar = np.sqrt(z1[0]*z1[0]+z1[1]*z1[1]+z1[2]*z1[2]) dstar = np.sqrt(z1[0] * z1[0] + z1[1] * z1[1] + z1[2] * z1[2])
if dstar > 0.0001: if dstar > 0.0001:
ds = 1/dstar ds = 1 / dstar
sint = wave * dstar/2 sint = wave * dstar / 2
if np.abs(sint) <= 1: if np.abs(sint) <= 1:
th = np.arcsin(sint)*pir th = np.arcsin(sint) * pir
else: else:
ierr = 2 ierr = 2
th = 0 th = 0
@ -209,7 +217,7 @@ def dandth(wave,z1):
return ds, th, ierr return ds, th, ierr
def angs4c(wave,z1,ch2,ph2): def angs4c(wave, z1, ch2, ph2):
"""CALCULATE 2-THETA, OMEGA (=THETA), CHI, PHI TO PUT THE """CALCULATE 2-THETA, OMEGA (=THETA), CHI, PHI TO PUT THE
VECTOR Z1 IN THE BISECTING DIFFRACTION CONDITION VECTOR Z1 IN THE BISECTING DIFFRACTION CONDITION
@ -222,10 +230,10 @@ def angs4c(wave,z1,ch2,ph2):
ch2, ph2 = eqchph(z1) ch2, ph2 = eqchph(z1)
ch = ch2 ch = ch2
ph = ph2 ph = ph2
ds, th, ierr = dandth(wave,z1) ds, th, ierr = dandth(wave, z1)
if ierr == 0: if ierr == 0:
om = th om = th
tth = th*2 tth = th * 2
else: else:
tth = 0 tth = 0
om = 0 om = 0
@ -235,7 +243,7 @@ def angs4c(wave,z1,ch2,ph2):
return tth, om, ch, ph, ierr return tth, om, ch, ph, ierr
def fixdnu(wave,z1,ch2,ph2,nu): def fixdnu(wave, z1, ch2, ph2, nu):
"""CALCULATE A SETTING CH,PH,GA,OM TO PUT THE DIFFRACTED BEAM AT NU. """CALCULATE A SETTING CH,PH,GA,OM TO PUT THE DIFFRACTED BEAM AT NU.
PH PUTS THE DIFFRACTION VECTOR Z1 INTO THE CHI CIRCLE (AS FOR PH PUTS THE DIFFRACTION VECTOR Z1 INTO THE CHI CIRCLE (AS FOR
BISECTING GEOMETRY), CH BRINGS THE VECTOR TO THE APPROPRIATE NU BISECTING GEOMETRY), CH BRINGS THE VECTOR TO THE APPROPRIATE NU
@ -247,9 +255,9 @@ def fixdnu(wave,z1,ch2,ph2,nu):
Returns: Returns:
tth, om, ch, ph tth, om, ch, ph
""" """
pir = 180/np.pi pir = 180 / np.pi
tth, om, ch, ph, ierr = angs4c(wave,z1,ch2,ph2) tth, om, ch, ph, ierr = angs4c(wave, z1, ch2, ph2)
theta = om theta = om
if ierr != 0: if ierr != 0:
ch = 0 ch = 0
@ -257,15 +265,15 @@ def fixdnu(wave,z1,ch2,ph2,nu):
ga = 0 ga = 0
om = 0 om = 0
else: else:
if np.abs(np.cos(nu/pir)) > 0.0001: if np.abs(np.cos(nu / pir)) > 0.0001:
cosga = np.cos(tth/pir)/np.cos(nu/pir) cosga = np.cos(tth / pir) / np.cos(nu / pir)
if np.abs(cosga) <= 1: if np.abs(cosga) <= 1:
ga = np.arccos(cosga)*pir ga = np.arccos(cosga) * pir
z4 = z4frgn(wave,ga,nu) z4 = z4frgn(wave, ga, nu)
om = np.arctan2(-z4[1],z4[0])*pir om = np.arctan2(-z4[1], z4[0]) * pir
ch2 = np.arcsin(z4[2]*wave/(2*np.sin(theta/pir)))*pir ch2 = np.arcsin(z4[2] * wave / (2 * np.sin(theta / pir))) * pir
ch = ch - ch2 ch = ch - ch2
ch = ch - 360 * np.trunc((np.sign(ch)*180+ch)/360) ch = ch - 360 * np.trunc((np.sign(ch) * 180 + ch) / 360)
else: else:
ierr = -2 ierr = -2
ch = 0 ch = 0
@ -276,9 +284,9 @@ def fixdnu(wave,z1,ch2,ph2,nu):
if theta > 44.99 and theta < 45.01: if theta > 44.99 and theta < 45.01:
ga = 90 ga = 90
om = 90 om = 90
ch2 = np.sign(nu)*45 ch2 = np.sign(nu) * 45
ch = ch - ch2 ch = ch - ch2
ch = ch - 360 * np.trunc((np.sign(ch)*180+ch)/360) ch = ch - 360 * np.trunc((np.sign(ch) * 180 + ch) / 360)
else: else:
ierr = -1 ierr = -1
ch = 0 ch = 0
@ -289,10 +297,11 @@ def fixdnu(wave,z1,ch2,ph2,nu):
return ch, ph, ga, om return ch, ph, ga, om
#for test run: # for test run:
# angtohkl(wave=1.18,ddist=616,gammad=48.66,om=-22.80,ch=0,ph=0,nud=0,x=128,y=64) # angtohkl(wave=1.18,ddist=616,gammad=48.66,om=-22.80,ch=0,ph=0,nud=0,x=128,y=64)
def angtohkl(wave,ddist,gammad,om,ch,ph,nud,x,y):
def angtohkl(wave, ddist, gammad, om, ch, ph, nud, x, y):
"""finds hkl-indices of a reflection from its position (x,y,angles) at the 2d-detector """finds hkl-indices of a reflection from its position (x,y,angles) at the 2d-detector
Args: Args:
@ -301,33 +310,72 @@ def angtohkl(wave,ddist,gammad,om,ch,ph,nud,x,y):
Returns: Returns:
""" """
pir = 180/np.pi pir = 180 / np.pi
#define ub matrix if testing angtohkl(wave=1.18,ddist=616,gammad=48.66,om=-22.80,ch=0,ph=0,nud=0,x=128,y=64) against f90: # define ub matrix if testing angtohkl(wave=1.18,ddist=616,gammad=48.66,om=-22.80,ch=0,ph=0,nud=0,x=128,y=64) against f90:
# ub = np.array([-0.0178803,-0.0749231,0.0282804,-0.0070082,-0.0368001,-0.0577467,0.1609116,-0.0099281,0.0006274]).reshape(3,3) # ub = np.array([-0.0178803,-0.0749231,0.0282804,-0.0070082,-0.0368001,-0.0577467,0.1609116,-0.0099281,0.0006274]).reshape(3,3)
ub = np.array([0.04489,0.02045,-0.2334,-0.06447,0.00129,-0.16356,-0.00328,0.2542,0.0196]).reshape(3,3) ub = np.array(
print('The input values are: ga=', gammad, ', om=', om, ', ch=', ch, ', ph=', ph, ', nu=', nud, ', x=', x, ', y=', y) [0.04489, 0.02045, -0.2334, -0.06447, 0.00129, -0.16356, -0.00328, 0.2542, 0.0196]
).reshape(3, 3)
print(
"The input values are: ga=",
gammad,
", om=",
om,
", ch=",
ch,
", ph=",
ph,
", nu=",
nud,
", x=",
x,
", y=",
y,
)
ga, nu = det2pol(ddist,gammad,nud,x,y) ga, nu = det2pol(ddist, gammad, nud, x, y)
print('The calculated actual angles are: ga=', ga, ', om=', om, ', ch=', ch, ', ph=', ph, ', nu=', nu) print(
"The calculated actual angles are: ga=",
ga,
", om=",
om,
", ch=",
ch,
", ph=",
ph,
", nu=",
nu,
)
z1 = z1frmd(wave,ga,om,ch,ph,nu) z1 = z1frmd(wave, ga, om, ch, ph, nu)
print('The diffraction vector is:', z1[0],z1[1],z1[2]) print("The diffraction vector is:", z1[0], z1[1], z1[2])
ubinv = np.linalg.inv(ub) ubinv = np.linalg.inv(ub)
h = ubinv[0,0]*z1[0]+ubinv[0,1]*z1[1]+ubinv[0,2]*z1[2] h = ubinv[0, 0] * z1[0] + ubinv[0, 1] * z1[1] + ubinv[0, 2] * z1[2]
k = ubinv[1,0]*z1[0]+ubinv[1,1]*z1[1]+ubinv[1,2]*z1[2] k = ubinv[1, 0] * z1[0] + ubinv[1, 1] * z1[1] + ubinv[1, 2] * z1[2]
l = ubinv[2,0]*z1[0]+ubinv[2,1]*z1[1]+ubinv[2,2]*z1[2] l = ubinv[2, 0] * z1[0] + ubinv[2, 1] * z1[1] + ubinv[2, 2] * z1[2]
print('The Miller indexes are:', h,k,l) print("The Miller indexes are:", h, k, l)
ch2, ph2 = eqchph(z1) ch2, ph2 = eqchph(z1)
ch, ph, ga, om = fixdnu(wave,z1,ch2,ph2,nu) ch, ph, ga, om = fixdnu(wave, z1, ch2, ph2, nu)
print('Bisecting angles to put reflection into the detector center: ga=', ga, ', om=', om, ', ch=', ch, ', ph=', ph, ', nu=', nu) print(
"Bisecting angles to put reflection into the detector center: ga=",
ga,
", om=",
om,
", ch=",
ch,
", ph=",
ph,
", nu=",
nu,
)
def ang2hkl(wave, ddist, gammad, om, ch, ph, nud, ub, x, y): def ang2hkl(wave, ddist, gammad, om, ch, ph, nud, ub, x, y):
@ -340,6 +388,7 @@ def ang2hkl(wave, ddist, gammad, om, ch, ph, nud, ub, x, y):
return hkl return hkl
def gauss(x, *p): def gauss(x, *p):
"""Defines Gaussian function """Defines Gaussian function
@ -350,9 +399,10 @@ def gauss(x, *p):
Gaussian function Gaussian function
""" """
A, mu, sigma = p A, mu, sigma = p
return A*np.exp(-(x-mu)**2/(2.*sigma**2)) return A * np.exp(-((x - mu) ** 2) / (2.0 * sigma ** 2))
def box_int(file,box):
def box_int(file, box):
"""Calculates center of the peak in the NB-geometry angles and Intensity of the peak """Calculates center of the peak in the NB-geometry angles and Intensity of the peak
Args: Args:
@ -362,72 +412,72 @@ def box_int(file,box):
gamma, omPeak, nu polar angles, Int and data for 3 fit plots gamma, omPeak, nu polar angles, Int and data for 3 fit plots
""" """
dat=pyzebra.read_detector_data(file) dat = pyzebra.read_detector_data(file)
sttC=dat["pol_angle"][0] sttC = dat["pol_angle"][0]
om=dat["rot_angle"] om = dat["rot_angle"]
nuC=dat["tlt_angle"][0] nuC = dat["tlt_angle"][0]
ddist=dat["ddist"] ddist = dat["ddist"]
# defining indices # defining indices
i0, j0, iN, jN, fr0, frN = box i0, j0, iN, jN, fr0, frN = box
# omega fit # omega fit
om=dat["rot_angle"][fr0:frN] om = dat["rot_angle"][fr0:frN]
cnts = np.sum(dat["data"][fr0:frN,j0:jN,i0:iN], axis=(1,2)) cnts = np.sum(dat["data"][fr0:frN, j0:jN, i0:iN], axis=(1, 2))
p0 = [1., 0., 1.] p0 = [1.0, 0.0, 1.0]
coeff, var_matrix = curve_fit(gauss, range(len(cnts)), cnts, p0=p0) coeff, var_matrix = curve_fit(gauss, range(len(cnts)), cnts, p0=p0)
frC = fr0+coeff[1] frC = fr0 + coeff[1]
omF = dat["rot_angle"][math.floor(frC)] omF = dat["rot_angle"][math.floor(frC)]
omC = dat["rot_angle"][math.ceil(frC)] omC = dat["rot_angle"][math.ceil(frC)]
frStep = frC-math.floor(frC) frStep = frC - math.floor(frC)
omStep = omC-omF omStep = omC - omF
omP = omF + omStep*frStep omP = omF + omStep * frStep
Int = coeff[1]*abs(coeff[2]*omStep)*math.sqrt(2)*math.sqrt(np.pi) Int = coeff[1] * abs(coeff[2] * omStep) * math.sqrt(2) * math.sqrt(np.pi)
# omega plot # omega plot
x_fit = np.linspace(0, len(cnts), 100) x_fit = np.linspace(0, len(cnts), 100)
y_fit = gauss(x_fit, *coeff) y_fit = gauss(x_fit, *coeff)
plt.figure() plt.figure()
plt.subplot(131) plt.subplot(131)
plt.plot(range(len(cnts)), cnts) plt.plot(range(len(cnts)), cnts)
plt.plot(x_fit, y_fit) plt.plot(x_fit, y_fit)
plt.ylabel('Intensity in the box') plt.ylabel("Intensity in the box")
plt.xlabel('Frame N of the box') plt.xlabel("Frame N of the box")
label='om' label = "om"
# gamma fit # gamma fit
sliceXY = dat["data"][fr0:frN,j0:jN,i0:iN] sliceXY = dat["data"][fr0:frN, j0:jN, i0:iN]
sliceXZ = np.sum(sliceXY,axis=1) sliceXZ = np.sum(sliceXY, axis=1)
sliceYZ = np.sum(sliceXY,axis=2) sliceYZ = np.sum(sliceXY, axis=2)
projX = np.sum(sliceXZ, axis=0) projX = np.sum(sliceXZ, axis=0)
p0 = [1., 0., 1.] p0 = [1.0, 0.0, 1.0]
coeff, var_matrix = curve_fit(gauss, range(len(projX)), projX, p0=p0) coeff, var_matrix = curve_fit(gauss, range(len(projX)), projX, p0=p0)
x= i0+coeff[1] x = i0 + coeff[1]
# gamma plot # gamma plot
x_fit = np.linspace(0, len(projX), 100) x_fit = np.linspace(0, len(projX), 100)
y_fit = gauss(x_fit, *coeff) y_fit = gauss(x_fit, *coeff)
plt.subplot(132) plt.subplot(132)
plt.plot(range(len(projX)), projX) plt.plot(range(len(projX)), projX)
plt.plot(x_fit, y_fit) plt.plot(x_fit, y_fit)
plt.ylabel('Intensity in the box') plt.ylabel("Intensity in the box")
plt.xlabel('X-pixel of the box') plt.xlabel("X-pixel of the box")
# nu fit # nu fit
projY = np.sum(sliceYZ, axis=0) projY = np.sum(sliceYZ, axis=0)
p0 = [1., 0., 1.] p0 = [1.0, 0.0, 1.0]
coeff, var_matrix = curve_fit(gauss, range(len(projY)), projY, p0=p0) coeff, var_matrix = curve_fit(gauss, range(len(projY)), projY, p0=p0)
y= j0+coeff[1] y = j0 + coeff[1]
# nu plot # nu plot
x_fit = np.linspace(0, len(projY), 100) x_fit = np.linspace(0, len(projY), 100)
y_fit = gauss(x_fit, *coeff) y_fit = gauss(x_fit, *coeff)
plt.subplot(133) plt.subplot(133)
plt.plot(range(len(projY)), projY) plt.plot(range(len(projY)), projY)
plt.plot(x_fit, y_fit) plt.plot(x_fit, y_fit)
plt.ylabel('Intensity in the box') plt.ylabel("Intensity in the box")
plt.xlabel('Y-pixel of the box') plt.xlabel("Y-pixel of the box")
ga, nu = pyzebra.det2pol(ddist,sttC,nuC,x,y) ga, nu = pyzebra.det2pol(ddist, sttC, nuC, x, y)
return ga[0], omP, nu[0], Int # x0,y0,x0_fit,y0_fit,x1,y1,x1_fit,y1_fit,x2,y2,x2_fit,y2_fit return ga[0], omP, nu[0], Int