Skip to content

Commit 48e3d53

Browse files
newvillenewville
committed
interpolate, if needed, MCAs to first energy array
1 parent 40e5cc0 commit 48e3d53

File tree

1 file changed

+127
-99
lines changed

1 file changed

+127
-99
lines changed

‎larch/xrmmap/xrm_mapfile.py

Lines changed: 127 additions & 99 deletions
Original file line numberDiff line numberDiff line change
@@ -6,11 +6,12 @@
66
import numpy as np
77
from pathlib import Path
88
import scipy.stats as stats
9+
from scipy.interpolate import interp1d
910
import json
1011
import multiprocessing as mp
1112
from functools import partial
1213

13-
from pyshortcuts import fix_varname, fix_filename, bytes2str
14+
from pyshortcuts import fix_varname, fix_filename, bytes2str, debugtimer
1415

1516
from larch.utils import isotime, version_ge, unixpath
1617

@@ -272,7 +273,7 @@ def __init__(self, filename=None, folder=None, create_empty=False,
272273
self.force_no_dtc = False
273274
self.all_mcas = all_mcas
274275
self.detector_list = None
275-
276+
self.mca_energies = None
276277
self.compress_args = {'compression': compression}
277278
if compression != 'lzf':
278279
self.compress_args['compression_opts'] = compression_opts
@@ -370,7 +371,7 @@ def __init__(self, filename=None, folder=None, create_empty=False,
370371
cfile = FastMapConfig()
371372
cfile.Read(Path(self.folder, self.ScanFile))
372373
cfile.config['scan']['filename'] = self.filename
373-
print("Create HDF5 File ")
374+
# print("Create HDF5 File ")
374375
self.h5root = h5py.File(self.filename, 'w')
375376
self.write_access = True
376377
if self.dimension is None and isGSEXRM_MapFolder(self.folder):
@@ -396,7 +397,7 @@ def __init__(self, filename=None, folder=None, create_empty=False,
396397
xrd1dbkgdfile=xrd1dbkgd)
397398
elif (self.filename is not None and
398399
self.status == GSEXRM_FileStatus.err_notfound and create_empty):
399-
print("Create HDF5 File")
400+
# print("Create HDF5 File")
400401
self.h5root = h5py.File(self.filename, 'w')
401402
self.write_access = True
402403
create_xrmmap(self.h5root, root=None, dimension=2, start_time=self.start_time)
@@ -598,6 +599,7 @@ def add_data(self, group, name, data, attrs=None, **kws):
598599
kws.update(self.compress_args)
599600
if name in group:
600601
del group[name]
602+
# print("Add Data ", group, name, type(data), kws)
601603
d = group.create_dataset(name, data=data, **kws)
602604
if isinstance(attrs, dict):
603605
for key, val in attrs.items():
@@ -640,6 +642,7 @@ def add_map_config(self, config, nmca=None):
640642
if os.path.exists(roifile):
641643
roidat, calib, extra = readROIFile(roifile)
642644
self.xrmmap.attrs['N_Detectors'] = self.nmca = len(calib['slope'])
645+
# print("CALIB ", calib)
643646
roi_desc, roi_addr, roi_lim = [], [], []
644647
roi_slices = []
645648

@@ -670,6 +673,7 @@ def add_map_config(self, config, nmca=None):
670673
self.roi_addr = roi_addr
671674
self.roi_slices = roi_slices
672675
self.calib = calib
676+
# print("SET CALIB ", calib)
673677
else:
674678
nmca = self.nmca
675679
roilims = np.array([ [[0, 1]] for i in range(nmca)])
@@ -867,7 +871,6 @@ def read_rowdata(self, irow, offset=None):
867871
self.xrdcalfile = calfile
868872

869873
scan_version = getattr(self, 'scan_version', 1.00)
870-
print(" read row data, scan version ", scan_version, self.xrdcalfile)
871874
xrdcal_dat = bytes2str(self.xrmmap['xrd1d'].attrs.get('caldata','{}'))
872875
if self.xrdcalfile is not None and len(xrdcal_dat) < 10:
873876
xrdcal_dat = read_poni(self.xrdcalfile)
@@ -931,7 +934,7 @@ def read_rowdata(self, irow, offset=None):
931934

932935
def add_rowdata(self, row, callback=None, flush=True):
933936
'''adds a row worth of real data'''
934-
# dt = debugtimer()
937+
dt = debugtimer()
935938
if not self.check_hostid():
936939
raise GSEXRM_Exception(NOT_OWNER % self.filename)
937940
if not self.write_access:
@@ -948,50 +951,65 @@ def add_rowdata(self, row, callback=None, flush=True):
948951
pform = '%s, xrdfile=%s' % (pform, row.xrdfile)
949952
print(pform)
950953

951-
# dt.add(" ran callback, print, version %s" %self.version)
954+
dt.add(" ran callback, print, version %s" %self.version)
952955

953-
if version_ge(self.version, '2.0.0'):
954956

955-
mcasum_raw,mcasum_cor = [],[]
957+
if version_ge(self.version, '2.0.0'):
958+
mcasum_raw, mcasum_cor = [],[]
956959
nrows = 0
957960
map_items = sorted(self.xrmmap.keys())
958-
# dt.add(" get %d map items" % len(map_items))
961+
dt.add(" get %d map items" % len(map_items))
959962
for gname in map_items:
960963
g = self.xrmmap[gname]
961964
if bytes2str(g.attrs.get('type', '')).startswith('scalar detect'):
962965
first_det = list(g.keys())[0]
963966
nrows, npts = g[first_det].shape
964967

965-
# dt.add(" got %d map items" % len(map_items))
968+
dt.add(" got %d map items" % len(map_items))
966969
if thisrow >= nrows:
967970
self.resize_arrays(NINIT*(1+nrows/NINIT), force_shrink=False)
968971

969-
# dt.add(" resized ")
972+
dt.add(" resized ")
970973
sclrgrp = self.xrmmap['scalars']
971974
for ai, aname in enumerate(row.scaler_names):
972975
sclrgrp[aname][thisrow, :npts] = row.sisdata[:npts].transpose()[ai]
973-
# dt.add(" add scaler group")
976+
dt.add(" add scaler group")
974977
if self.has_xrf:
975-
976978
npts = min([len(p) for p in row.posvals])
977979
pos = self.xrmmap['positions/pos']
978980
rowpos = np.array([p[:npts] for p in row.posvals])
979981

980982
tpos = rowpos.transpose()
981983
pos[thisrow, :npts, :] = tpos[:npts, :]
982984
nmca, xnpts, nchan = row.counts.shape
985+
986+
if self.mca_energies is None:
987+
off = self.calib['offset']
988+
slo = self.calib['slope']
989+
xoff = max(off) - min(off)
990+
xslo = max(slo) - min(slo)
991+
if xslo > 1.e-6 or xoff > 1.e-4: # needs recalibration
992+
self.mca_energies = []
993+
enx = np.arange(nchan)
994+
for i in range(nmca):
995+
off = self.calib['offset'][i]
996+
slope = self.calib['slope'][i]
997+
quad = self.calib['quad'][i]
998+
self.mca_energies.append(off + enx*(slope+quad*enx))
999+
1000+
9831001
mca_dets = []
984-
# dt.add(" map xrf 1")
1002+
dt.add(" map xrf 1")
9851003
for gname in map_items:
9861004
g = self.xrmmap[gname]
9871005
if bytes2str(g.attrs.get('type', '')).startswith('mca detect'):
9881006
mca_dets.append(gname)
9891007
nrows, npts, nchan = g['counts'].shape
990-
# dt.add(" map xrf 2")
1008+
dt.add(" map xrf 2")
9911009
_nr, npts, nchan = self.xrmmap['mcasum']['counts'].shape
1010+
9921011
npts = min(npts, xnpts, self.npts)
993-
# dt.add(" map xrf 3")
994-
# print("ADD ROW ", self.all_mcas, mca_dets, self.nmca)
1012+
dt.add(" map xrf 3")
9951013
if self.all_mcas:
9961014
for idet, gname in enumerate(mca_dets):
9971015
grp = self.xrmmap[gname]
@@ -1007,7 +1025,7 @@ def add_rowdata(self, row, callback=None, flush=True):
10071025
dtfactor = np.zeros(npts, dtype=np.float32)
10081026
inpcounts = np.zeros(npts, dtype=np.float32)
10091027
outcounts = np.zeros(npts, dtype=np.float32)
1010-
# dt.add(" map xrf 4a: alloc ")
1028+
dt.add(" map xrf 4a: alloc ")
10111029
# print("ADD ", inpcounts.dtype, row.inpcounts.dtype)
10121030
for idet in range(self.nmca):
10131031
realtime += row.realtime[idet, :npts]
@@ -1016,18 +1034,27 @@ def add_rowdata(self, row, callback=None, flush=True):
10161034
outcounts += row.outcounts[idet, :npts]
10171035
livetime /= (1.0*self.nmca)
10181036
realtime /= (1.0*self.nmca)
1019-
# dt.add(" map xrf 4b: time sums")
1037+
dt.add(" map xrf 4b: time sums")
10201038

10211039
sumgrp = self.xrmmap['mcasum']
1022-
sumgrp['counts'][thisrow, :npts, :nchan] = row.total[:npts, :nchan]
1023-
# dt.add(" map xrf 4b: set counts")
1024-
# print("add realtime ", sumgrp['realtime'].shape, self.xrmmap['roimap/det_raw'].shape, thisrow)
1040+
ensum = sumgrp['energy']
1041+
if self.mca_energies is None:
1042+
sumgrp['counts'][thisrow, :npts, :nchan] = row.total[:npts, :nchan]
1043+
else:
1044+
sumx = row.counts[0, :npts, :nchan]*0.0
1045+
for i in range(self.nmca):
1046+
en = self.mca_energies[i]
1047+
sumx += interp1d(en, row.counts[i], kind='linear', fill_value=0,
1048+
copy=False, bounds_error=False)(ensum)
1049+
sumgrp['counts'][thisrow, :npts, :nchan] = sumx
1050+
1051+
dt.add(" map xrf 4b: set counts")
10251052
sumgrp['realtime'][thisrow, :npts] = realtime
10261053
sumgrp['livetime'][thisrow, :npts] = livetime
10271054
sumgrp['dtfactor'][thisrow, :npts] = row.total_dtfactor[:npts]
10281055
sumgrp['inpcounts'][thisrow, :npts] = inpcounts
10291056
sumgrp['outcounts'][thisrow, :npts] = outcounts
1030-
# dt.add(" map xrf 4c: set time data ")
1057+
dt.add(" map xrf 4c: set time data ")
10311058

10321059
if version_ge(self.version, '2.1.0'): # version 2.1
10331060
det_raw = self.xrmmap['roimap/det_raw']
@@ -1059,7 +1086,7 @@ def add_rowdata(self, row, callback=None, flush=True):
10591086
detcor.extend(icor)
10601087
sumraw.append(np.array(iraw).sum(axis=0))
10611088
sumcor.append(np.array(icor).sum(axis=0))
1062-
# dt.add(" map xrf 5a: got simple ROIS")
1089+
dt.add(" map xrf 5a: got simple ROIS")
10631090
det_raw[thisrow, :npts, :] = np.array(detraw).transpose()
10641091
det_cor[thisrow, :npts, :] = np.array(detcor).transpose()
10651092
sum_raw[thisrow, :npts, :] = np.array(sumraw).transpose()
@@ -1084,80 +1111,81 @@ def add_rowdata(self, row, callback=None, flush=True):
10841111
sumcor += mcacor
10851112
roigrp['mcasum'][roiname]['raw'][thisrow,] = sumraw
10861113
roigrp['mcasum'][roiname]['cor'][thisrow,] = sumcor
1087-
# dt.add(" map xrf 6")
1114+
dt.add(" map xrf 6")
10881115
else: # version 1.0.1
1089-
if self.has_xrf:
1090-
nmca, xnpts, nchan = row.counts.shape
1091-
xrm_dets = []
1092-
1093-
nrows = 0
1094-
map_items = sorted(self.xrmmap.keys())
1095-
for gname in map_items:
1096-
g = self.xrmmap[gname]
1097-
if bytes2str(g.attrs.get('type', '')).startswith('mca detect'):
1098-
xrm_dets.append(g)
1099-
nrows, npts, nchan = g['counts'].shape
1100-
1101-
if thisrow >= nrows:
1102-
self.resize_arrays(NINIT*(1+nrows/NINIT), force_shrink=False)
1103-
1104-
_nr, npts, nchan = xrm_dets[0]['counts'].shape
1105-
npts = min(npts, xnpts, self.npts)
1106-
for idet, grp in enumerate(xrm_dets):
1107-
grp['dtfactor'][thisrow, :npts] = row.dtfactor[idet, :npts]
1108-
grp['realtime'][thisrow, :npts] = row.realtime[idet, :npts]
1109-
grp['livetime'][thisrow, :npts] = row.livetime[idet, :npts]
1110-
grp['inpcounts'][thisrow, :npts] = row.inpcounts[idet, :npts]
1111-
grp['outcounts'][thisrow, :npts] = row.outcounts[idet, :npts]
1112-
grp['counts'][thisrow, :npts, :] = row.counts[idet, :npts, :]
1113-
1114-
# here, we add the total dead-time-corrected data to detsum.
1115-
self.xrmmap['detsum']['counts'][thisrow, :npts, :nchan] = row.total[:npts, :nchan]
1116-
1117-
pos = self.xrmmap['positions/pos']
1118-
rowpos = np.array([p[:npts] for p in row.posvals])
1119-
1120-
tpos = rowpos.transpose()
1121-
1122-
pos[thisrow, :npts, :] = tpos[:npts, :]
1123-
1124-
# now add roi map data
1125-
roimap = self.xrmmap['roimap']
1126-
det_raw = roimap['det_raw']
1127-
det_cor = roimap['det_cor']
1128-
sum_raw = roimap['sum_raw']
1129-
sum_cor = roimap['sum_cor']
1130-
1131-
detraw = list(row.sisdata[:npts].transpose())
1132-
1133-
detcor = detraw[:]
1134-
sumraw = detraw[:]
1135-
sumcor = detraw[:]
1136-
1137-
if self.roi_slices is None:
1138-
lims = self.xrmmap['config/rois/limits'][()]
1139-
nrois, nmca, nx = lims.shape
1140-
1141-
self.roi_slices = []
1142-
for iroi in range(nrois):
1143-
x = [slice(lims[iroi, i, 0],
1144-
lims[iroi, i, 1]) for i in range(nmca)]
1145-
self.roi_slices.append(x)
1146-
1147-
for slices in self.roi_slices:
1148-
iraw = [row.counts[i, :npts, slices[i]].sum(axis=1)
1149-
for i in range(nmca)]
1150-
icor = [row.counts[i, :npts, slices[i]].sum(axis=1)*row.dtfactor[i, :npts]
1151-
for i in range(nmca)]
1152-
detraw.extend(iraw)
1153-
detcor.extend(icor)
1154-
sumraw.append(np.array(iraw).sum(axis=0))
1155-
sumcor.append(np.array(icor).sum(axis=0))
1156-
1157-
det_raw[thisrow, :npts, :] = np.array(detraw).transpose()
1158-
det_cor[thisrow, :npts, :] = np.array(detcor).transpose()
1159-
sum_raw[thisrow, :npts, :] = np.array(sumraw).transpose()
1160-
sum_cor[thisrow, :npts, :] = np.array(sumcor).transpose()
1116+
print("version1?")
1117+
# if self.has_xrf:
1118+
# nmca, xnpts, nchan = row.counts.shape
1119+
# xrm_dets = []
1120+
#
1121+
# nrows = 0
1122+
# map_items = sorted(self.xrmmap.keys())
1123+
# for gname in map_items:
1124+
# g = self.xrmmap[gname]
1125+
# if bytes2str(g.attrs.get('type', '')).startswith('mca detect'):
1126+
# xrm_dets.append(g)
1127+
# nrows, npts, nchan = g['counts'].shape
1128+
#
1129+
# if thisrow >= nrows:
1130+
# self.resize_arrays(NINIT*(1+nrows/NINIT), force_shrink=False)
1131+
#
1132+
# _nr, npts, nchan = xrm_dets[0]['counts'].shape
1133+
# npts = min(npts, xnpts, self.npts)
1134+
# for idet, grp in enumerate(xrm_dets):
1135+
# grp['dtfactor'][thisrow, :npts] = row.dtfactor[idet, :npts]
1136+
# grp['realtime'][thisrow, :npts] = row.realtime[idet, :npts]
1137+
# grp['livetime'][thisrow, :npts] = row.livetime[idet, :npts]
1138+
# grp['inpcounts'][thisrow, :npts] = row.inpcounts[idet, :npts]
1139+
# grp['outcounts'][thisrow, :npts] = row.outcounts[idet, :npts]
1140+
# grp['counts'][thisrow, :npts, :] = row.counts[idet, :npts, :]
1141+
#
1142+
# # here, we add the total dead-time-corrected data to detsum.
1143+
# self.xrmmap['detsum']['counts'][thisrow, :npts, :nchan] = row.total[:npts, :nchan]
1144+
#
1145+
# pos = self.xrmmap['positions/pos']
1146+
# rowpos = np.array([p[:npts] for p in row.posvals])
1147+
#
1148+
# tpos = rowpos.transpose()
1149+
#
1150+
# pos[thisrow, :npts, :] = tpos[:npts, :]
1151+
#
1152+
# # now add roi map data
1153+
# roimap = self.xrmmap['roimap']
1154+
# det_raw = roimap['det_raw']
1155+
# det_cor = roimap['det_cor']
1156+
# sum_raw = roimap['sum_raw']
1157+
# sum_cor = roimap['sum_cor']
1158+
#
1159+
# detraw = list(row.sisdata[:npts].transpose())
1160+
#
1161+
# detcor = detraw[:]
1162+
# sumraw = detraw[:]
1163+
# sumcor = detraw[:]
1164+
#
1165+
# if self.roi_slices is None:
1166+
# lims = self.xrmmap['config/rois/limits'][()]
1167+
# nrois, nmca, nx = lims.shape
1168+
#
1169+
# self.roi_slices = []
1170+
# for iroi in range(nrois):
1171+
# x = [slice(lims[iroi, i, 0],
1172+
# lims[iroi, i, 1]) for i in range(nmca)]
1173+
# self.roi_slices.append(x)
1174+
#
1175+
# for slices in self.roi_slices:
1176+
# iraw = [row.counts[i, :npts, slices[i]].sum(axis=1)
1177+
# for i in range(nmca)]
1178+
# icor = [row.counts[i, :npts, slices[i]].sum(axis=1)*row.dtfactor[i, :npts]
1179+
# for i in range(nmca)]
1180+
# detraw.extend(iraw)
1181+
# detcor.extend(icor)
1182+
# sumraw.append(np.array(iraw).sum(axis=0))
1183+
# sumcor.append(np.array(icor).sum(axis=0))
1184+
#
1185+
# det_raw[thisrow, :npts, :] = np.array(detraw).transpose()
1186+
# det_cor[thisrow, :npts, :] = np.array(detcor).transpose()
1187+
# sum_raw[thisrow, :npts, :] = np.array(sumraw).transpose()
1188+
# sum_cor[thisrow, :npts, :] = np.array(sumcor).transpose()
11611189

11621190
if self.has_xrd1d and row.xrdq is not None:
11631191
if thisrow < 2:
@@ -1194,7 +1222,7 @@ def add_rowdata(self, row, callback=None, flush=True):
11941222
self.last_row = thisrow
11951223
self.xrmmap.attrs['Last_Row'] = thisrow
11961224
#self.h5root.flush()
1197-
# dt.add("flushed h5 file")
1225+
dt.add("flushed h5 file")
11981226
# dt.show()
11991227

12001228

0 commit comments

Comments
 (0)