Skip to content

Commit 5bb642e

Browse files
tothsatothsa
tothsa
authored and
tothsa
committed
developing Horace capabilities
1 parent 4989811 commit 5bb642e

9 files changed

+443
-41
lines changed

‎load_dat.m

+121
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,121 @@
1+
function D = load_dat(fDir,fName,varargin)
2+
% load TOF spectrum from HDF5 file
3+
4+
bin = varargin(1:4);
5+
6+
dir0 = pwd;
7+
cd(fDir);
8+
9+
% read data dimensions
10+
par = h5info(fName,'/specnd/datcnt');
11+
nQ = par.Dataspace.Size(1);
12+
nE = par.Dataspace.Size(2);
13+
% load grid position information
14+
gridPos = h5read(fName,'/specnd/gridpos');
15+
% add last position afte the end of the file
16+
gridPos = [gridPos nQ*nE+1];
17+
gridVal = h5read(fName,'/specnd/gridval');
18+
19+
EN = h5read(fName,'/specnd/ebin');
20+
21+
22+
% load the saved Q grid for fast data loading
23+
grid0 = h5read(fName,'/specnd/qbin');
24+
grid = cell(1,3);
25+
nGridEdge = zeros(1,3);
26+
27+
for ii = 1:3
28+
grid{ii} = grid0(ii,~isnan(grid0(ii,:)));
29+
nGridEdge(ii) = numel(grid{ii});
30+
end
31+
grid{4} = EN;
32+
nGridEdge(4) = numel(grid{4});
33+
34+
% determine the necessary data cubes to load
35+
idxGrid = [ones(4,1) nGridEdge'-1];
36+
for ii = 1:4
37+
idx = find(grid{ii}<=bin{ii}(1),1,'last');
38+
if ~isempty(idx)
39+
idxGrid(ii,1) = idx;
40+
end
41+
idx = find(grid{ii}>=bin{ii}(end),1,'first')-1;
42+
if ~isempty(idx)
43+
idxGrid(ii,2) = idx;
44+
end
45+
end
46+
47+
% load the necessary data
48+
vGrid = cell(1,4);
49+
for ii = 1:4
50+
vGrid{ii} = idxGrid(ii,1):idxGrid(ii,2);
51+
end
52+
53+
% create the grid list of data blocks to load
54+
gGrid = cell(1,4);
55+
[gGrid{:}] = ndgrid(vGrid{:});
56+
57+
listGrid = sub2ind(nGridEdge-1,gGrid{:});
58+
% data start/end positions in the hdf file
59+
selGrid = ismember(gridVal,listGrid);
60+
startPos = gridPos(selGrid);
61+
endPos = gridPos([false selGrid])-1;
62+
% remove empty blocks
63+
rmIdx = startPos>endPos;
64+
startPos(rmIdx) = [];
65+
endPos(rmIdx) = [];
66+
67+
% load data
68+
dat0 = zeros(sum(endPos-startPos+1),1);
69+
err0 = zeros(sum(endPos-startPos+1),1);
70+
EN0 = zeros(sum(endPos-startPos+1),1);
71+
axi0 = zeros(3,sum(endPos-startPos+1));
72+
73+
runIdx = 1;
74+
75+
fprintf('Loading pixels...\n');
76+
77+
78+
sw_status(0,1)
79+
nPos = numel(startPos);
80+
81+
perc = 0;
82+
for ii = 1:nPos
83+
84+
[idxQ1,idxE1] = ind2sub([nQ nE],startPos(ii));
85+
[idxQ2,~] = ind2sub([nQ nE],endPos(ii));
86+
nDat = idxQ2-idxQ1+1;
87+
dat0(runIdx+(1:nDat)-1) = double(h5read(fName,'/specnd/datcnt',[idxQ1 idxE1],[nDat 1]));
88+
err0(runIdx+(1:nDat)-1) = double(h5read(fName,'/specnd/errmon',[idxQ1 idxE1],[nDat 1]));
89+
axi0(:,runIdx+(1:nDat)-1) = double(h5read(fName,'/specnd/axis', [1 idxQ1 idxE1],[3 nDat 1]));
90+
EN0(runIdx+(1:nDat)-1) = EN(idxE1);
91+
runIdx = runIdx + nDat;
92+
93+
if floor(ii/nPos*100)>perc
94+
perc = floor(ii/nPos*100);
95+
sw_status(ii/nPos*100);
96+
end
97+
end
98+
% load all data
99+
% dat0 = reshape(double(h5read(fName,'/specnd/datcnt',[1 1],[nQ nE])),[nQ*nE 1]);
100+
% err0 = reshape(double(h5read(fName,'/specnd/errmon',[1 1],[nQ nE])),[nQ*nE 1]);
101+
% axi0 = reshape(double(h5read(fName,'/specnd/axis', [1 1 1],[3 nQ nE])),[3 nQ*nE]);
102+
% EN0 = reshape(double(repmat(EN(1:end-1),[nQ 1])),[nQ*nE 1]);
103+
104+
sw_status(100,2);
105+
106+
% remove NaN data values
107+
nanidx = isnan(dat0);
108+
dat0(nanidx) = [];
109+
EN0(nanidx) = [];
110+
err0(nanidx) = [];
111+
axi0(:,nanidx) = [];
112+
113+
fprintf('Pixels are loaded into memory!\n');
114+
fprintf('Binning pixels...\n');
115+
% bin data
116+
D = specnd([axi0' EN0],dat0,err0);
117+
D.bin(varargin{:});
118+
fprintf('Binning ready!\n');
119+
cd(dir0);
120+
121+
end

‎make_horace.m

+79
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,79 @@
1+
%% directories
2+
3+
par_file = '~/Desktop/YBCF_data_80meV/ARCS_2x1_grouping.par';
4+
sqw_file = '~/Desktop/YBCF_data_80meV/ybcf_80meV.sqw';
5+
tmpdir = '~/Desktop/YBCF_data_80meV/horace_tmp/';
6+
indir = '~/Desktop/YBCF_data_80meV/';
7+
8+
%% crystal parameters
9+
10+
%direct geometry
11+
emode = 1;
12+
alatt = [3.875 3.875 7.679];
13+
angdeg = [90,90,90];
14+
% vector || incident beam
15+
u = [1 0 0];
16+
% vector perpendicular to the incident beam,
17+
% pointing towards the large angle detectors on Merlin in the horizontal plane
18+
v = [0 1 0];
19+
%offset angles in case of crystal misorientation (see the Horace manual for details)
20+
omega = 0.0;
21+
%dpsi = 0.0;
22+
%gl = 0.0;
23+
%gs = 0.0;
24+
25+
% Corrected values
26+
dpsi = -2.64;
27+
gl = -0.87;
28+
gs = -1.31;
29+
30+
31+
%% Ei = 80 meV, normalized
32+
33+
fName = dir([indir '*.nxspe']);
34+
fName = {fName.name};
35+
36+
% convert data to sqw
37+
nfiles = numel(fName);
38+
psi = zeros(nfiles,1);
39+
efix = zeros(nfiles,1);
40+
spe_file = cell(nfiles,1);
41+
42+
for idx = 1:nfiles
43+
spe_file{idx} = [indir fName{idx}];
44+
efix(idx) = h5read(spe_file{idx},'/data/NXSPE_info/fixed_energy');
45+
%psi(idx) = h5read(spe_file{idx},'/data/NXSPE_info/psi');
46+
47+
[~,tName] = fileparts(fName{idx});
48+
temp = sscanf(strrep(strrep(tName(13:end),'p','.'),'_',' '),'%f %f %f');
49+
T(idx) = temp(2);
50+
psi(idx) = temp(3);
51+
52+
53+
end
54+
fprintf('Number of files ready to convert: %d.\n',nfiles);
55+
56+
t1 = clock;
57+
gen_sqw(spe_file,par_file,sqw_file,efix,emode,alatt,angdeg,u,v,psi,omega,dpsi,gl,gs);
58+
t2 = clock;
59+
etime(t2,t1)
60+
61+
%% bin
62+
63+
proj.uoffset = [0 0 0];
64+
proj.type = 'ppr';
65+
% (H,K,0) scattering plane
66+
proj.u = [1 0 0];
67+
proj.v = [0 1 0];
68+
69+
% symmetrize all data
70+
Q1 = [-2 0.04 2];
71+
Q2 = [-2 0.04 2];
72+
Q3 = [-6 6];
73+
E0 = [20 25];
74+
75+
t1 = clock;
76+
dat1s = cut_sqw(sqw_file,proj,Q1,Q2,Q3,E0);
77+
t2 = clock;
78+
etime(t2,t1)
79+
% 16s s

‎mfiles/@specnd/plot2.m

+3-2
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@ function plot2(D,varargin)
66
% bins with single element gets a +/- 0.5 width
77
for ii = 1:2
88
if numel(cBin{ii}) == 1
9-
eBin{ii} = cBin{ii}+[-0.5 0.5];
9+
eBin{ii} = cBin{ii}+[-0.5 0.5]';
1010
else
1111
cc = cBin{ii};
1212
eBin{ii} = [(3*cc(1)-cc(2))/2; (cc(1:end-1)+cc(2:end))/2; (3*cc(end)-cc(end-1))/2];
@@ -22,7 +22,8 @@ function plot2(D,varargin)
2222
dat(end+1,:)=0;
2323
dat(:,end+1) = 0;
2424

25-
surf(xx,yy,dat);
25+
sHandle = surf(xx,yy,dat);
26+
set(sHandle,'edgealpha',0)
2627
axis([eBin{1}([1 end])' eBin{2}([1 end])']);
2728
view(2)
2829
xlabel(D.raw.axis.label{1})

‎mfiles/@specnd/specnd.m

+3-3
Original file line numberDiff line numberDiff line change
@@ -82,8 +82,8 @@
8282
D.raw.datcnt.value = varargin{2};
8383
D.raw.errmon.value = 0*D.raw.datcnt.value;
8484
case 3
85-
D.raw.datcnt.value = varargin{1};
86-
D.raw.errmon.value = varargin{2};
85+
D.raw.datcnt.value = varargin{2};
86+
D.raw.errmon.value = varargin{3};
8787
end
8888

8989
% check for event mode
@@ -118,7 +118,7 @@
118118
D.raw.datcnt.value = D.raw.datcnt.value(selector{:});
119119

120120
end
121-
121+
axDim = numel(D.raw.axis.value);
122122
end
123123

124124
if numel(D.raw.axis.name) < axDim

‎mfiles/@specnd/squeeze.m

+20
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
function squeeze(D)
2+
% remove data dimensions with a single coordinate value
3+
4+
if ishistmode(D)
5+
% find singleton indices
6+
singlIdx = find(size(D) == 1);
7+
dimVect = 1:ndims(D);
8+
dimVect(singlIdx) = [];
9+
D.raw.datcnt.value = squeeze(permute(D.raw.datcnt.value,[dimVect singlIdx]));
10+
D.raw.errmon.value = squeeze(permute(D.raw.errmon.value,[dimVect singlIdx]));
11+
12+
D.raw.axis.value = D.raw.axis.value(dimVect);
13+
D.raw.axis.name = D.raw.axis.name(dimVect);
14+
D.raw.axis.label = D.raw.axis.label(dimVect);
15+
16+
else
17+
% TODO
18+
end
19+
20+
end

‎num3compress.m

-14
This file was deleted.

‎num3decompress.m

-22
This file was deleted.

0 commit comments

Comments
 (0)