Deconvolution of the thickness effect¶
This is the third part of this long journey in the thickness of the silicon sensor of the Pilatus detector: After characterizing the precise position of the Pilatus-1M on the goniometer of ID28 we noticed there are some discrepancies in the ring position and peak width as function of the position on the detector. In a second time the thickness of the detector was modeled and allowed us to define a sparse-matrix which represent the bluring of the signal with a point-spread function which actually depends on the position of the detector. This convolution can be revereted using techniques developped for inverse problems.
We will now correct the images of the first notebook called “goniometer” with the sparse matrix calculated in the second one (called “raytracing”) and check if the pick-width is less chaotic.
In [1]:
%pylab nbagg
Populating the interactive namespace from numpy and matplotlib
In [2]:
import os
#Specific for ESRF
#os.environ["http_proxy"] = "http://proxy.esrf.fr:3128"
In [3]:
import fabio, pyFAI, os
from os.path import basename
from pyFAI.gui import jupyter
from pyFAI.calibrant import get_calibrant
from silx.resources import ExternalResources
downloader = ExternalResources("thick", "http://www.silx.org/pub/pyFAI/testimages")
all_files = downloader.getdir("gonio_ID28.tar.bz2")
for afile in all_files:
print(basename(afile))
gonio_ID28
det130_g45_0001p.npt
det130_g0_0001p.cbf
det130_g17_0001p.npt
det130_g0_0001p.npt
det130_g45_0001p.cbf
det130_g17_0001p.cbf
In [4]:
from scipy.sparse import save_npz, load_npz, csr_matrix, csc_matrix, linalg
#Saved in notebook called "raytracing"
csr = load_npz("csr.npz")
In [5]:
mask = numpy.load("mask.npy")
images = [fabio.open(i).data for i in sorted(all_files) if i.endswith("cbf")]
fig, ax = subplots(1, 3, figsize=(9,3))
for i, img in enumerate(images):
jupyter.display(img, ax=ax[i], label=str(i))
Data type cannot be displayed: application/javascript
In [6]:
fig, ax = subplots(1, 3, figsize=(9,3))
msk = numpy.where(mask)
for i, img in enumerate(images):
fl = img.astype("float32")
fl[msk] = 0.0 # set masked values to 0, negative values could induce errors
bl = fl.ravel()
res = linalg.lsmr(csr.T, bl)
print(res[1:])
cor = res[0].reshape(img.shape)
jupyter.display(cor, ax=ax[i], label=str(i))
Data type cannot be displayed: application/javascript
(1, 23, 66.70987399605569, 10.158270983765981, 1.862080889148464, 4.1457039198038652, 44893427.485308647)
(1, 23, 58.916625196254635, 9.5709507555976376, 1.8666609017532882, 4.0146845391231922, 38902815.355951048)
(1, 24, 28.59430045453111, 4.7422569251327822, 1.8989417893747016, 3.9796981410935794, 16639924.368782267)
It turns out the deconvolution is not that straight forwards and creates some negative wobbles near the rings. This phenomenon is well known in inverse methods which provide a damping factor to limit the effect. This damping factor needs to be adjusted manually to avoid this.
In [7]:
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets
In [8]:
img = images[0]
fl = img.astype("float32")
fl[msk] = 0.0 # set masked values to 0
blured = fl.ravel()
fig, ax = subplots(1, 2, figsize=(8,4))
im0 = ax[0].imshow(numpy.arcsinh(fl))
im1 = ax[1].imshow(numpy.arcsinh(fl))
ax[0].set_title("Original")
ax[1].set_title("Deconvolved")
def deconvol(damp):
res = linalg.lsmr(csr.T, blured, damp=damp, x0=blured)
restored = res[0].reshape(mask.shape)
im1.set_data(numpy.arcsinh(restored))
print("Number of negative pixels: %i"%(restored<0).sum())
interactive_plot = widgets.interactive(deconvol, damp=(0, 5.0))
display(interactive_plot)
Data type cannot be displayed: application/javascript
Failed to display Jupyter Widget of type interactive.
If you're reading this message in Jupyter Notebook or JupyterLab, it may mean that the widgets JavaScript is still loading. If this message persists, it likely means that the widgets JavaScript library is either not installed or not enabled. See the Jupyter Widgets Documentation for setup instructions.
If you're reading this message in another notebook frontend (for example, a static rendering on GitHub or NBViewer), it may mean that your frontend doesn't currently support widgets.
In [9]:
#selection of the damping factor which provides no negative signal:
tot_neg = 100
damp = 0
while tot_neg>0:
damp += 0.1
tot_neg = 0
deconvoluted = []
for i, img in enumerate(images):
fl = img.astype("float32")
fl[msk] = 0.0 # set masked values to 0
bl = fl.ravel()
res = linalg.lsmr(csr.T, bl, damp=damp)
neg=(res[0]<0).sum()
print(i, damp, neg, res[1:])
tot_neg += neg
deconvoluted.append(res[0].reshape(img.shape))
print("damp:",damp)
0 0.1 12878 (2, 19, 4321883.593997509, 7.0740130809038249, 1.698393943810682, 3.4183219524917772, 41657214.505458035)
1 0.1 12253 (2, 20, 3763677.5121062007, 3.9064069822047269, 1.7452525730482935, 3.4658210409578594, 36450685.767392009)
2 0.1 17600 (2, 21, 1583946.5629618806, 1.5743847474835946, 1.778888073385638, 3.5335925761415083, 15124376.494351184)
0 0.2 11660 (2, 14, 7888848.985405592, 5.6256288146636457, 1.4657231247188438, 2.723491733119078, 34985942.907106549)
1 0.2 10761 (2, 14, 6940345.515941161, 5.069311731600906, 1.4706687568610501, 2.7315252177713973, 31183825.221334115)
2 0.2 15000 (2, 14, 2845633.3087043264, 3.1371831766796401, 1.4635183619841836, 2.7248874167756072, 12406590.276522264)
0 0.30000000000000004 9562 (2, 11, 10555837.680557573, 4.2084004259290628, 1.3062417828607145, 2.2651195844181076, 28199160.306425255)
1 0.30000000000000004 9076 (2, 11, 9374221.821183499, 3.6531515632123646, 1.3123978620478625, 2.2399738034325511, 25515415.77776191)
2 0.30000000000000004 12266 (2, 11, 3768064.0522358264, 2.1385262209007752, 1.3033266565553487, 2.2309901394162353, 9874730.9157402962)
0 0.4 8059 (2, 9, 12476186.72001587, 3.8802413730896794, 1.1904619320515117, 1.8863968666028703, 22440000.641676854)
1 0.4 7441 (2, 9, 11156347.45337748, 3.4289630350685742, 1.193761729637159, 1.8838616188409905, 20505761.17148133)
2 0.4 9597 (2, 9, 4425091.2015181035, 1.9388546751796458, 1.1824568018372001, 1.8855234891543979, 7802642.8943203576)
0 0.5 6476 (2, 7, 13846834.842369597, 13.263426650611899, 1.0568874444572918, 1.6323276625305716, 17877903.969196949)
1 0.5 5692 (2, 7, 12441366.2042652, 11.348221224974166, 1.0615548249406035, 1.6339004608457397, 16439162.569837855)
2 0.5 7095 (2, 8, 4891168.198022051, 1.0345853896447619, 1.1168383072609167, 1.6321674724241761, 6188789.7184904292)
0 0.6 4923 (2, 7, 14831285.463075664, 2.383871906669496, 1.0568874444572918, 1.4278004347472653, 14361329.60120401)
1 0.6 4022 (2, 7, 13370023.65804597, 2.0472902634859418, 1.0615548249406035, 1.4290684915424074, 13258771.307992153)
2 0.6 4842 (2, 7, 5224588.54627925, 1.0793039647024578, 1.0487777157964333, 1.4265615929562203, 4956588.9562078137)
0 0.7 3443 (2, 6, 15548310.636277616, 4.6723322304274255, 0.9840543618324515, 1.2645032189796421, 11669711.721603643)
1 0.7 2684 (2, 6, 14049027.610942049, 4.1096615813816468, 0.9898624060121443, 1.2669688761269708, 10802702.566435993)
2 0.7 3086 (2, 6, 5466763.484668505, 2.1026340070406078, 0.9751918481711386, 1.264136638465831, 4019140.024491996)
0 0.7999999999999999 2237 (2, 6, 16079695.284954334, 1.3294073199019898, 0.9840543618324515, 1.1337895832243896, 9602257.0319290441)
1 0.7999999999999999 1704 (2, 6, 14553505.51696309, 1.1725314068967438, 0.9898624060121443, 1.1355211262715892, 8905330.036925476)
2 0.7999999999999999 1871 (2, 6, 5645880.542526024, 0.59424705441206371, 0.9751918481711386, 1.1335209261469905, 3302039.1182649173)
0 0.8999999999999999 1348 (2, 5, 16480747.611802146, 5.88848599824325, 0.9054518506346702, 1.039865943648604, 8000075.5414817166)
1 0.8999999999999999 1041 (2, 5, 14934905.21790837, 5.3131121168099753, 0.9129611496133215, 1.0436940315508072, 7429263.867318524)
2 0.8999999999999999 1107 (2, 5, 5780866.319376789, 2.6618523453105669, 0.895637461222218, 1.0391140348035892, 2747949.5313641271)
0 0.9999999999999999 767 (2, 5, 16788870.296727166, 2.4060587913606692, 0.9054518506346702, 1.0326083316265215, 6744334.3023225199)
1 0.9999999999999999 632 (2, 5, 15228285.34073762, 2.1756707956130987, 0.9129611496133215, 1.0357683855630029, 6269229.1360698529)
2 0.9999999999999999 674 (2, 5, 5884458.333562568, 1.0837286606434389, 0.895637461222218, 1.0322579088364725, 2314610.1023986028)
0 1.0999999999999999 433 (2, 5, 17029593.71573765, 1.0471308411680613, 0.9054518506346702, 1.0271348620122089, 5748096.3353474392)
1 1.0999999999999999 352 (2, 5, 15457694.549778676, 0.94848713848455624, 0.9129611496133215, 1.0297855860872132, 5347107.3793465523)
2 1.0999999999999999 384 (2, 5, 5965320.9551829165, 0.47035647210922549, 0.895637461222218, 1.0270157355540273, 1971384.3568484683)
0 1.2 245 (2, 4, 17220587.87557631, 11.855358405370115, 0.818337233068597, 1.0228447483444125, 4948101.2751028389)
1 1.2 211 (2, 4, 15639833.203990966, 10.667667735287811, 0.8280449594647742, 1.0240351125170375, 4605539.0509930383)
2 1.2 237 (2, 5, 6029435.580222552, 0.21609741363975563, 0.895637461222218, 1.0229300886008843, 1696117.4126032754)
0 1.3 143 (2, 4, 17374277.48592824, 6.6438913328505063, 0.818337233068597, 1.0195514395489986, 4298170.1464291876)
1 1.3 111 (2, 4, 15786472.117560865, 5.9864846986379838, 0.8280449594647742, 1.0205829336798344, 4002395.5081131491)
2 1.3 155 (2, 4, 6080999.829532508, 2.9845210580302393, 0.8083501992645462, 1.0196907865808638, 1472708.3768822292)
0 1.4000000000000001 90 (2, 4, 17499542.775204774, 3.8582334846833568, 0.818337233068597, 1.016915565901702, 3764334.2263497142)
1 1.4000000000000001 76 (2, 4, 15906038.829056432, 3.4803379823921365, 0.8280449594647742, 1.0178173471693504, 3506552.3747703591)
2 1.4000000000000001 96 (2, 4, 6123009.27499327, 1.7310892564689488, 0.8083501992645462, 1.0170829998730164, 1289352.8012345743)
0 1.5000000000000002 56 (2, 4, 17602833.310866546, 2.3131510977289076, 0.818337233068597, 1.0147747747162661, 3321358.5302527202)
1 1.5000000000000002 46 (2, 4, 16004662.178194474, 2.0885021356719804, 0.8280449594647742, 1.0155693862838346, 3094813.4956444195)
2 1.5000000000000002 65 (2, 4, 6157636.924453535, 1.0368336672318383, 0.8083501992645462, 1.0149548617756046, 1137303.558753982)
0 1.6000000000000003 38 (2, 4, 17688905.85443919, 1.4271357216035403, 0.818337233068597, 1.0130134143579892, 2950291.7275948594)
1 1.6000000000000003 30 (2, 4, 16086866.980071092, 1.2895178180100904, 0.8280449594647742, 1.0137185412720142, 2749716.2542595086)
2 1.6000000000000003 39 (2, 4, 6186483.87907338, 0.63917261239093048, 0.8083501992645462, 1.0131969317880869, 1010004.7138191476)
0 1.7000000000000004 29 (2, 4, 17761319.91984709, 0.90354319007651873, 0.818337233068597, 1.011547454110197, 2636742.3656661543)
1 1.7000000000000004 18 (2, 4, 16156041.903024914, 0.81693823315851599, 0.8280449594647742, 1.0121771505991219, 2457974.5777046219)
2 1.7000000000000004 28 (2, 4, 6210747.298094724, 0.40439646963781672, 0.8083501992645462, 1.0117289283879789, 902485.71152076731)
0 1.8000000000000005 23 (2, 4, 17822775.505701326, 0.58558193750495946, 0.818337233068597, 1.0103147366772813, 2369661.0411364655)
1 1.8000000000000005 10 (2, 4, 16214759.014770027, 0.5297413075741314, 0.8280449594647742, 1.0108981361916489, 2209372.3870919882)
2 1.8000000000000005 15 (2, 4, 6231334.663883276, 0.2619374331399405, 0.8083501992645462, 1.010490990259741, 810935.24911589024)
0 1.9000000000000006 17 (2, 4, 17875346.688535817, 0.38765774648737256, 0.818337233068597, 1.0092685343141001, 2140473.6901405477)
1 1.9000000000000006 6 (2, 4, 16264995.160215642, 0.35085381371694896, 0.8280449594647742, 1.0098134547256175, 1995972.6073321986)
2 1.9000000000000006 8 (2, 4, 6248942.6985554835, 0.17331921935393643, 0.8083501992645462, 1.0094378095984924, 732399.02764947049)
0 2.0000000000000004 15 (2, 4, 17920645.541291054, 0.26164541750381859, 0.818337233068597, 1.0083731994109995, 1942459.1804616847)
1 2.0000000000000004 5 (2, 4, 16308287.509445166, 0.23689964239881259, 0.8280449594647742, 1.0088816353525585, 1811547.2201532284)
2 2.0000000000000004 2 (2, 4, 6264112.675685391, 0.11693102290937513, 0.8083501992645462, 1.0085346211593638, 664563.27804327966)
0 2.1000000000000005 14 (2, 4, 17959938.7724845, 0.17974708690925484, 0.818337233068597, 1.0076011719777833, 1770297.9325031796)
1 2.1000000000000005 5 (2, 3, 16345844.356045963, 11.527255800288904, 0.7313672012567888, 1.0081646332914587, 1651162.8364928684)
2 2.1000000000000005 1 (2, 4, 6277269.739256002, 0.080301092826323264, 0.8083501992645462, 1.0077544125555722, 605598.10830568522)
0 2.2000000000000006 13 (2, 3, 17994231.832293328, 10.263410253899691, 0.7208770294028594, 1.0068508486240206, 1619741.2455065292)
1 2.2000000000000006 3 (2, 3, 16378625.079487013, 8.8104017698063171, 0.7313672012567888, 1.0074479285894242, 1510876.4487291707)
2 2.2000000000000006 0 (2, 3, 6288751.232052226, 4.3585883233108191, 0.7124343875322071, 1.0071559405479567, 554042.94457297574)
0 2.3000000000000007 11 (2, 3, 18024330.32737255, 7.9301559985591359, 0.7208770294028594, 1.0062728549795792, 1487366.5904062653)
1 2.3000000000000007 1 (2, 3, 16407398.575351965, 6.8091967991917768, 0.7313672012567888, 1.0068212584313183, 1387509.8940314848)
2 2.3000000000000007 0 (2, 3, 6298827.366561279, 3.3670362457128267, 0.7124343875322071, 1.0065496190645882, 508721.879663839)
0 2.400000000000001 10 (2, 3, 18050885.401900202, 6.1907738951009987, 0.7208770294028594, 1.0057648372544621, 1370394.7243476359)
1 2.400000000000001 1 (2, 3, 16432786.456492491, 5.3168717291679908, 0.7313672012567888, 1.0062702064648275, 1278481.0515201602)
2 2.400000000000001 0 (2, 3, 6307716.489964505, 2.6280485363813355, 0.7124343875322071, 1.0060170858006161, 468680.48472360586)
0 2.500000000000001 7 (2, 3, 18074427.647790488, 4.879198689130428, 0.7208770294028594, 1.0053159707485797, 1266551.6874910567)
1 2.500000000000001 1 (2, 3, 16455295.350526728, 4.1912710752911124, 0.7313672012567888, 1.0057831196576557, 1181676.3168828806)
2 2.500000000000001 0 (2, 3, 6315596.484650351, 2.0709433172317415, 0.7124343875322071, 1.0055468583778133, 433138.17947539786)
0 2.600000000000001 6 (2, 3, 18095392.711381465, 3.8796648067317228, 0.7208770294028594, 1.0049174397556742, 1173963.7080113278)
1 2.600000000000001 1 (2, 3, 16475341.299150117, 3.3332508036151904, 0.7313672012567888, 1.0053504985950636, 1095353.3870023305)
2 2.600000000000001 0 (2, 3, 6322613.37317426, 1.6464660276376597, 0.7124343875322071, 1.005129599401849, 401451.99252611451)
0 2.700000000000001 4 (2, 3, 18114140.821997482, 3.1103838672934323, 0.7208770294028594, 1.0045620074969484, 1091076.4586229334)
1 2.700000000000001 0 (2, 3, 16493268.372779518, 2.6727367455912159, 0.7313672012567888, 1.0049645385526969, 1018066.5058819719)
2 2.700000000000001 0 (2, 3, 6328887.878079462, 1.3198298265370041, 0.7124343875322071, 1.0047576540723313, 373088.74073076935)
0 2.800000000000001 4 (2, 3, 18130971.824247207, 2.5128476424079511, 0.7208770294028594, 1.0042436894335907, 1016592.4992789062)
1 2.800000000000001 0 (2, 3, 16509363.003949437, 2.1595821567826863, 0.7313672012567888, 1.0046187803518472, 948608.50462896307)
2 2.800000000000001 0 (2, 3, 6334520.4696034575, 1.0661568781892281, 0.7124343875322071, 1.0044246997593191, 347603.48945320735)
0 2.9000000000000012 3 (2, 3, 18146136.85095344, 2.0447118582107167, 0.7208770294028594, 1.003957502589699, 949422.42121676472)
1 2.9000000000000012 0 (2, 3, 16523865.120650608, 1.7574823256596444, 0.7313672012567888, 1.0043078420813099, 885965.50644858391)
2 2.9000000000000012 0 (2, 3, 6339595.284117616, 0.86744702770234017, 0.7124343875322071, 1.0041254776778017, 324622.74059592187)
0 3.0000000000000013 2 (2, 3, 18159847.46242451, 1.674995373508624, 0.7208770294028594, 1.0036992713812964, 888646.40054799616)
1 3.0000000000000013 0 (2, 3, 16536976.865283042, 1.4398671279835318, 0.7313672012567888, 1.0040272110346629, 829281.26220905629)
2 3.0000000000000013 0 (2, 3, 6344183.19120208, 0.71053358433247538, 0.7124343875322071, 1.0038555854417486, 303831.20969523396)
0 3.1000000000000014 0 (2, 3, 18172282.85776532, 1.3807986212804906, 0.7208770294028594, 1.0034654758863191, 833483.72448448092)
1 3.1000000000000014 0 (2, 3, 16548869.475819526, 1.1870921704503978, 0.7313672012567888, 1.0037730809420273, 777828.86774469551)
2 3.1000000000000014 0 (2, 3, 6348344.213110899, 0.58568651378386238, 0.7124343875322071, 1.0036113152406074, 284961.35003635875)
damp: 3.1000000000000014
In [10]:
# damp = 1
# deconvoluted = []
# for i, img in enumerate(images):
# fl = img.astype("float32")
# fl[msk] = 0.0 # set masked values to 0
# bl = fl.ravel()
# res = linalg.lsmr(csr.T, bl, damp=damp)
# neg=(res[0]<-1).sum()
# print(i, damp, neg, res[1:])
# tot_neg += neg
# deconvoluted.append(res[0].reshape(img.shape))
In [11]:
fig, ax = subplots(1, 3, figsize=(9,3))
for i,img in enumerate(deconvoluted):
jupyter.display(img, ax=ax[i], label=str(i))
Data type cannot be displayed: application/javascript
At this point we are back to the initial state: can we fit the goniometer and check the width of the peaks to validate if it got better.
In [12]:
wavelength=0.6968e-10
calibrant = get_calibrant("LaB6")
calibrant.wavelength = wavelength
print(calibrant)
detector = pyFAI.detector_factory("Pilatus1M")
detector.mask = mask
fimages = []
fimages.append(fabio.cbfimage.CbfImage(data=deconvoluted[0], header={"angle":0}))
fimages.append(fabio.cbfimage.CbfImage(data=deconvoluted[1], header={"angle":17}))
fimages.append(fabio.cbfimage.CbfImage(data=deconvoluted[2], header={"angle":45}))
LaB6 Calibrant at wavelength 6.968e-11
In [13]:
def get_angle(metadata):
"""Takes the basename (like det130_g45_0001.cbf ) and returns the angle of the detector"""
return metadata["angle"]
In [14]:
from pyFAI.goniometer import GeometryTransformation, GoniometerRefinement, Goniometer
from pyFAI.gui import jupyter
In [15]:
gonioref2d = GoniometerRefinement.sload("id28.json", pos_function=get_angle)
In [16]:
for idx, fimg in enumerate(fimages):
sg = gonioref2d.new_geometry(str(idx), image=fimg.data, metadata=fimg.header, calibrant=calibrant)
print(sg.label, "Angle:", sg.get_position())
0 Angle: 0
1 Angle: 17
2 Angle: 45
In [17]:
fig,ax = subplots(1,3, figsize=(9,3))
for lbl, sg in gonioref2d.single_geometries.items():
idx = int(lbl)
print(sg.extract_cp())
jupyter.display(sg=sg, ax=ax[idx])
Data type cannot be displayed: application/javascript
ControlPoints instance containing 6 group of point:
LaB6 Calibrant with 120 reflections at wavelength 6.968e-11
Containing 6 groups of points:
# a ring 0: 361 points
# b ring 1: 361 points
# c ring 2: 326 points
# d ring 3: 131 points
# e ring 4: 48 points
# f ring 5: 2 points
ControlPoints instance containing 13 group of point:
LaB6 Calibrant with 120 reflections at wavelength 6.968e-11
Containing 13 groups of points:
# g ring 0: 183 points
# h ring 1: 181 points
# i ring 2: 179 points
# j ring 3: 129 points
# k ring 4: 110 points
# l ring 5: 98 points
# m ring 6: 82 points
# n ring 7: 77 points
# o ring 8: 72 points
# p ring 9: 68 points
# q ring 10: 37 points
# r ring 11: 16 points
# s ring 12: 1 points
ControlPoints instance containing 26 group of point:
LaB6 Calibrant with 120 reflections at wavelength 6.968e-11
Containing 26 groups of points:
# t ring 7: 37 points
# u ring 8: 55 points
# v ring 9: 67 points
# w ring 10: 64 points
# x ring 11: 62 points
# y ring 12: 61 points
# z ring 13: 57 points
#aa ring 14: 55 points
#ab ring 15: 55 points
#ac ring 16: 53 points
#ad ring 17: 53 points
#ae ring 18: 51 points
#af ring 19: 51 points
#ag ring 20: 49 points
#ah ring 21: 47 points
#ai ring 22: 47 points
#aj ring 23: 45 points
#ak ring 24: 45 points
#al ring 25: 43 points
#am ring 26: 43 points
#an ring 27: 42 points
#ao ring 28: 41 points
#ap ring 29: 41 points
#aq ring 30: 41 points
#ar ring 31: 19 points
#as ring 32: 3 points
In [18]:
gonioref2d.refine2()
Cost function before refinement: 3.32962453706e-08
[ 2.84547605e-01 8.86540301e-02 8.93070277e-02 5.48717450e-03
5.54901577e-03 1.74619285e-02 -2.09889821e-05]
fun: 2.0418327956501834e-08
jac: array([ 1.13396290e-07, -3.58306167e-07, 3.13610199e-06,
-9.22968526e-07, 1.54877071e-07, -3.13002441e-05,
3.56611340e-08])
message: 'Optimization terminated successfully.'
nfev: 62
nit: 6
njev: 6
status: 0
success: True
x: array([ 2.84378750e-01, 8.86587707e-02, 8.93044104e-02,
5.48712194e-03, 5.54964221e-03, 1.74617227e-02,
-2.08467289e-05])
Cost function after refinement: 2.0418327956501834e-08
GonioParam(dist=0.28437875037614968, poni1=0.088658770743275148, poni2=0.08930441042948345, rot1=0.0054871219448940269, rot2=0.005549642205035156, scale1=0.017461722749243933, scale2=-2.0846728913328391e-05)
maxdelta on: dist (0) 0.2845476046521533 --> 0.284378750376
Out[18]:
array([ 2.84378750e-01, 8.86587707e-02, 8.93044104e-02,
5.48712194e-03, 5.54964221e-03, 1.74617227e-02,
-2.08467289e-05])
In [19]:
#Create a MultiGeometry integrator from the refined geometry:
angles = []
dimages = []
for sg in gonioref2d.single_geometries.values():
angles.append(sg.get_position())
dimages.append(sg.image)
multigeo = gonioref2d.get_mg(angles)
multigeo.radial_range=(0, 63)
print(multigeo)
# Integrate the whole set of images in a single run:
res_mg = multigeo.integrate1d(dimages, 10000)
ax = jupyter.plot1d(res_mg, label="multigeo")
for lbl, sg in gonioref2d.single_geometries.items():
ai = gonioref2d.get_ai(sg.get_position())
img = sg.image * ai.dist * ai.dist / ai.pixel1 / ai.pixel2
res = ai.integrate1d(img, 5000, unit="2th_deg", method="splitpixel")
ax.plot(*res, "--", label=lbl)
ax.legend()
MultiGeometry integrator with 3 geometries on (0, 63) radial range (2th_deg) and (-180, 180) azimuthal range (deg)
Data type cannot be displayed: application/javascript
Out[19]:
<matplotlib.legend.Legend at 0x7f82f7107c18>
In [20]:
#Let's focus on the inner most ring on the image taken at 45°:
#ax.set_xlim(21.5, 21.7)
ax.set_xlim(29.0, 29.2)
ax.set_ylim(0, 2e10)
Out[20]:
(0, 20000000000.0)
In [21]:
#Peak FWHM
from scipy.interpolate import interp1d
from scipy.optimize import bisect
def calc_fwhm(integrate_result, calibrant):
"calculate the tth position and FWHM for each peak"
delta = integrate_result.intensity[1:] - integrate_result.intensity[:-1]
maxima = numpy.where(numpy.logical_and(delta[:-1]>0, delta[1:]<0))[0]
minima = numpy.where(numpy.logical_and(delta[:-1]<0, delta[1:]>0))[0]
maxima += 1
minima += 1
tth = []
FWHM = []
for tth_rad in calibrant.get_2th():
tth_deg = tth_rad*integrate_result.unit.scale
if (tth_deg<=integrate_result.radial[0]) or (tth_deg>=integrate_result.radial[-1]):
continue
idx_theo = abs(integrate_result.radial-tth_deg).argmin()
id0_max = abs(maxima-idx_theo).argmin()
id0_min = abs(minima-idx_theo).argmin()
I_max = integrate_result.intensity[maxima[id0_max]]
I_min = integrate_result.intensity[minima[id0_min]]
tth_maxi = integrate_result.radial[maxima[id0_max]]
I_thres = (I_max + I_min)/2.0
if minima[id0_min]>maxima[id0_max]:
if id0_min == 0:
min_lo = integrate_result.radial[0]
else:
min_lo = integrate_result.radial[minima[id0_min-1]]
min_hi = integrate_result.radial[minima[id0_min]]
else:
if id0_min == len(minima) -1:
min_hi = integrate_result.radial[-1]
else:
min_hi = integrate_result.radial[minima[id0_min+1]]
min_lo = integrate_result.radial[minima[id0_min]]
f = interp1d(integrate_result.radial, integrate_result.intensity-I_thres)
tth_lo = bisect(f, min_lo, tth_maxi)
tth_hi = bisect(f, tth_maxi, min_hi)
FWHM.append(tth_hi-tth_lo)
tth.append(tth_deg)
return tth, FWHM
fig, ax = subplots()
ax.plot(*calc_fwhm(res_mg, calibrant), "o", label="multi")
for lbl, sg in gonioref2d.single_geometries.items():
ai = gonioref2d.get_ai(sg.get_position())
img = sg.image * ai.dist * ai.dist / ai.pixel1 / ai.pixel2
res = ai.integrate1d(img, 5000, unit="2th_deg", method="splitpixel")
t,w = calc_fwhm(res, calibrant=calibrant)
ax.plot(t, w,"-o", label=lbl)
ax.set_title("Peak shape as function of the angle")
ax.set_xlabel(res_mg.unit.label)
ax.legend()
Data type cannot be displayed: application/javascript
Out[21]:
<matplotlib.legend.Legend at 0x7f82f7128668>
conclusion:¶
The results do not look enhanced compared to the initial fit: neither the peak position, nor the FWHM looks enhanced. Maybe there is an error somewhere.
In [ ]: