The code.ill.fr has been recreated and upgraded with the latest version this weekend, If you encounter any problem please inform the Helpdesk.

Commit 8c6901c1 authored by Gonzalez, Miguel's avatar Gonzalez, Miguel
Browse files

Cleaning and improving calibration macros for D2B

parent f1e031f6
......@@ -12,36 +12,60 @@
;Vanadium data is required for blocks 1 and 2.
;A highly structured pattern is required for 3.
;This angle correction has only be done once and is not included here.
;
;There are now 2 functions to calculate the efficiencies:
;
;calib_eff_2d_iterative is an updated version of the standard iterative
;procedure and the first choice,
;
;calib_eff_2d_coarse is a more approximate version, but also faster and
;more robust, so it can be used in problematic cases
;
;---------------------------------------
;
;set variables; first and last vana numors and title
a=31524 & b=31547 & n=b-a+1
t='vana:'+string(a)+' -'+string(b)
a = 31524
b = 31547
n = b-a+1
t = 'vana:' + string(a) + ' -' + string(b)
;set current directory (variable c) to read calibration files in local directory
cd,current=c
cd, current=c
;make sure no detectors are excluded
dud,dets=[0]
;set calibartion file with no active zones or efficiencies set
calibration,file=c+'\d2bcal.raw'
dud, dets=[0]
;set calibration file with no active zones and efficiencies = 1
calibration, file = c + '\d2bcal.raw'
;read small angle and high angle data
w1=rdrun(a) & w2=rdrun(b)
;
;run active zones calibration
w3=calib_zone(w1,w2,2,'d2bcal.raw','d2bcal071.zone',t)
;***you have to inspect the zone calib file at this point***
;********to check that all zones are well defined***
;**************all '*'s have to be removed***
w1 = rdrun(a)
w2 = rdrun(b)
;determine active zones in detector
w3 = calib_zone(w1, w2, 2, 'd2bcal.raw', 'd2bcal.zone', t)
;you should inspect the zone calib file and check that all zones are well defined
;good limits are {5:20}, {115:130}, {125:140}, {235:250}
;a * marks any tube with a value outside those limits
;use w3 to see zone boundaries
;set calibration file with zones defined
calibration,file=c+'\d2bcal071.zone'
;
;run efficiency calibration
w4=calib_eff_2d_bms(a,n,20,'d2bcal071.zone','d2bcal071.2d',t)
;set calibration file with zones and efficiencies defined
calibration,file=c+'\d2bcal071.2d'
;
calibration, file = c + '\d2bcal.zone'
;run efficiency calibration (grouping pixels vertically in blocks of 8)
;varying the threshold limits (0.9 and 1.1 here) can improve the result
w4 = calib_eff_2d_iterative(a, n, 20, 8, 0.9, 1.1, 'd2bcal.zone', 'd2bcal.2d', t)
;call for coarse method
;w4 = calib_eff_2d_coarse(a, n, 20, 1, 0.5, 1.5, 'd2bcal.zone', 'd2bcal.2d', t)
;set new calibration file with zones and efficiencies defined
calibration, file= c + '\d2bcal.2d'
;read data to check efficiency correction
w5=rdrun(a) & w6=rdrun(b)
;
;If this calib file is OK, copy it to
;/home/cs/lambda/CALIBRATION
w5 = rdrun(a)
w6 = rdrun(b)
;if this calib file is OK, copy it to /home/cs/lambda/CALIBRATION
;with an adequate name (e.g. including the date or the cycle number)
FUNCTION calib_eff_2d, numor, n, wi, cal_fil_in, cal_fil_out, title
;+
;******* ************
;**
;** User callable
;** The call is e.g. w2=calib_eff_2d(51577, 5, 20, 'd2bcal.raw', 'd2bcal.eff', 'title')
;**
;** 51577 is the first scan, 5 is the number of scans, 20 is the number of a working workspace
;** w2 is an output workspace containing the table of efficiencies
;** d2bcal.raw is the old file, d2bcal.eff is the new one
;** data must be read in with a calibration file that defines the active zones
;** since efficiencies are to be determined for these zones,
;** *** but the efficiencies table in ***
;** *** the input calibration file should contain only 1's ***
;** this macro as for calib_eff_tubes_1d, but efficiencies per vertical block (height hb) rather than per tube
;** in "1d" version integrate over height at the start
;** ---------------------------------------------------
;** METHOD
;** ------
;** Sum data using initial efficiencies to get average diffraction pattern
;** Difference between each scan and the average gives efficiency per tube, height group and scan
;** Average efficiency over different scans gives new effciency per tube, height group
;** Repeat
;** Currently the beamstop is not omitted.
;** For this, identify tubes that scan the beam stop
;** Use tube_score(scan,tube)
;** ---------------------------------------------------
;-
;read the first numor and get array dimensions
datp=1
win=rdrun(numor,DATP=datp,w=wi)
print, n_elements(datp)
npts=(size(win))[1] & print, 'npoints ',npts
nh =(size(win))[2] & print, 'nh ',nh
nt=128 ;number of detector tubes
nsteps=npts/nt ;number of scan steps e.g.25, n is number of different scans
tolerance=1e-4 ;used to stop refinment cycle for efficiencies
max_cycle=30 ;max no. of cycles in efficiency calculation
min_cycle=3 ;min no. of cycles in efficiency calculation
hb=16 ;height of blocks for efficiency calculation
nb=nh/hb ;number of blocks in height
;define matrix of efficiencies per scan (eff) and global (effg & effgold) (new and old)
eff=fltarr(n,nt,nb) & effg=fltarr(nt,nb) & effgold=fltarr(nt,nb)
eff[*,*,*]=1.0 & effg[*,*]=1.0 & effgold[*,*]=2.0
effout =fltarr(nh,nt) & effout[*,*]=1.0
efflimit=1.5 ; limit to exclude bad detectors from sum
tube_score=intarr(nt) & tube_score[*]=0 ; array to describe good or bad tubes
;tube_score(41)=1 & tube_score(45)=1 & tube_score(101)=1 & tube_score(102)=1 & tube_score(121)=1 ; bad tubes to be omitted from total
;effg(41)=2.1 & effg(45)=2.1 & effg(101)=2.1 & effg(102)=2.1 & effg(121)=2.1 ; fix efficiency
;make array to store all calibration workspaces and enter the first scan, read the rest
w_all_3d=fltarr(n,npts,nh)
w_all_3d[0,*,*]=win
print, 'reading data'
for i=1,n-1 do begin
newnumor=numor+i
win=rdrun(newnumor)
w_all_3d[i,*,*]=win
endfor
;define arrays for storing summed/total counts
n_totpts=npts+(n-1)*nsteps ;& print, n_totpts
w_tot=fltarr(n_totpts)
i_tot=intarr(n_totpts)
;while loop to refine efficiencies, stop when they are no longer changing according to the tolerance
icount=0 & deviation=tolerance*2
dev_old=0. & dev_new=0. & del_dev=1.; variables to follow gradient of deviation, stop on sign change
while (icount lt max_cycle) and (abs(del_dev) gt 1.0e-5) do begin
icount=icount+1 ;& print, 'doing efficiencies ', icount
;calculate the total including efficiencies, initialise total and counter arrays
w_tot[*]=0.0 & i_tot[*]=0
for i=0,n-1 do begin ; loop thro scans
for j=0,nt-1 do begin ; loop thro tubes
if (tube_score[j] lt 1) then begin
alllowlim=nsteps*j & allhilim=alllowlim+nsteps-1 ;limits to define a single tube
totlowlim=nsteps*(i+j) & tothilim=totlowlim+nsteps-1 ;limits to define a single tube, including scan offset
for jh=0,nb-1 do begin ; loop thro height blocks
hmin=jh*hb & hmax=(jh+1)*hb-1
i_tot[totlowlim:tothilim]=i_tot[totlowlim:tothilim]+hb ; count how many times adding into w_tot
w_tot[totlowlim:tothilim]=w_tot[totlowlim:tothilim]+total(w_all_3d[i,alllowlim:allhilim,hmin:hmax],3)*effg[j,jh]
endfor
endif else begin
; if (i eq 0) then print,'detector missing', j
endelse
endfor
endfor
;calculate the average
w_tot[*]=w_tot[*]/i_tot[*]
;calculate differences between individual scans and the average i.e. efficiencies per scan i
;averaged over the number of scan steps
for i=0,n-1 do begin ; loop over different scans
for j=0,nt-1 do begin ; loop over detector tubes
if (tube_score[j] lt 1) then begin
alllowlim=nsteps*j & allhilim=alllowlim+nsteps-1 ;limits to define a single tube
totlowlim=nsteps*(i+j) & tothilim=totlowlim+nsteps-1 ;limits to define a single tube, including scan offset
for jh=0,nb-1 do begin ; loop thro height blocks, nb is number of blocks height hb
hmin=jh*hb & hmax=(jh+1)*hb-1 ; & print, 'hmin and hmax', hmin, hmax
eff[i,j,jh]=total(w_tot[totlowlim:tothilim])/total(w_all_3d[i,alllowlim:allhilim,hmin:hmax])
endfor
endif
endfor
endfor
;calulate average/global efficiency for detectors not omitted
for j=0,nt-1 do begin
for jh=0,nb-1 do begin
if (tube_score[j] lt 1) then effg[j,jh]=hb*total(eff[*,j,jh],1)/n
endfor
endfor
;calculate deviation from previous efficiencies
dev_new=total(abs(effg-effgold))/nt/hb & print,'deviation, icount ',dev_new, icount
del_dev=dev_old-dev_new & dev_old=dev_new & print, del_dev
;update previous efficiencies
effgold=effg
endwhile ;*** efficiencies refined in effg
ave = total(effg)/n_elements(effg) & print, 'ave efficiency ', ave
effg=effg/ave
;write effg into effout, the transpose
for j=0,nt-1 do begin
for jh=0,nb-1 do begin
hmin=jh*hb & hmax=(jh+1)*hb-1
effout[hmin:hmax,j]=effg[j,jh]
endfor
endfor
;normalise efficiencies
; ave=total(effout)/n_elements(effout) & effout=effout/ave & print, 'average eff is ',ave
; ****************** update calibration file ******************
;effout in format (nh,nt), write new calibration file
print, 'writing calibration file'
openr,f_in,cal_fil_in,/get_lun
openw,f_out,cal_fil_out,/get_lun
;write new zones info
printf,f_out,title
;copy zones from old calibration file
junk=''
readf,f_in,junk ;skip title
while junk ne 'efficiencies' do begin
readf,f_in,junk
printf,f_out,junk
endwhile
;write efficiencies
printf,f_out,effout
;read down to angles in old file
while junk ne 'angles' do readf,f_in,junk
;then write angles from old file
printf,f_out,'angles'
while junk ne 'calib_end' do begin
readf,f_in,junk
printf,f_out,junk
endwhile
;printf,f_out,'calib_end'
free_lun, f_in & free_lun, f_out
return, effout
end
FUNCTION calib_eff_2d_B, numor, n, wi, bms_minthreshold, bms_maxthreshold, cal_fil_in, cal_fil_out, title
;+
;******* ****************
;**
;** User callable
;** The call is e.g. w2=calib_eff_2d(51577, 5, 20, 0.9, 1.1, 'd2bcal.raw', 'd2bcal.eff', 'title')
;**
;** 51577 is the first scan, 5 is the number of scans, 20 is the number of a working workspace
;** 0.9 is the fraction of the average intensity used as lower threshold to determine tubes masked by beamstop,
;** 1.1 is the fraction of the average intensity used as upper threshold to determine tubes masked by beamstop,
;** w2 is an output workspace containing the table of efficiencies
;** d2bcal.raw is the old file, d2bcal.eff is the new one
;** data must be read in with a calibration file that defines the active zones
;** since efficiencies are to be determined for these zones,
;** *** but the efficiencies table in ***
;** *** the input calibration file should contain only 1's ***
;** this macro as for calib_eff_tubes_1d, but efficiencies per vertical block (height hb) rather than per tube
;** in "1d" version integrate over height at the start
;** ---------------------------------------------------
;** METHOD
;** ------
;** Sum data using initial efficiencies to get average diffraction pattern
;** Difference between each scan and the average gives efficiency per tube, height group and scan
;** Average efficiency over different scans gives new effciency per tube, height group
;** Repeat
;** Currently the beamstop is not omitted.
;** For this, identify tubes that scan the beam stop
;** Use tube_score(scan,tube)
;** ---------------------------------------------------
;-
forward_function rdrun
;read the first numor and get array dimensions
datp=1
win=rdrun(numor,DATP=datp,w=wi)
print, n_elements(datp)
npts=(size(win))[1] & print, 'npoints ',npts
nh =(size(win))[2] & print, 'nh ',nh
nt=128 ;number of detector tubes
nsteps=npts/nt ;number of scan steps e.g.25, n is number of different scans
tolerance=1e-4 ;used to stop refinment cycle for efficiencies
max_cycle=30 ;max no. of cycles in efficiency calculation
min_cycle=3 ;min no. of cycles in efficiency calculation
hb=16 ;height of blocks for efficiency calculation
nb=nh/hb ;number of blocks in height
;define matrix of efficiencies per scan (eff) and global (effg & effgold) (new and old)
eff=fltarr(n,nt,nb) & effg=fltarr(nt,nb) & effgold=fltarr(nt,nb)
eff[*,*,*]=1.0 & effg[*,*]=1.0 & effgold[*,*]=2.0
effout =fltarr(nh,nt) & effout[*,*]=1.0
;make array to store all calibration workspaces and enter the first scan, read the rest
w_all_3d=fltarr(n,npts,nh)
w_all_3d[0,*,*]=win
print, 'reading data'
for i=1,n-1 do begin
newnumor=numor+i
win=rdrun(newnumor)
w_all_3d[i,*,*]=win
endfor
;find average intensity and scale it, used to find beamstop
ave_min = bms_minthreshold* total(win)/npts/nb
ave_max = bms_maxthreshold * total(win)/npts/nb
;define arrays for storing summed/total counts
n_totpts=npts+(n-1)*nsteps ;& print, n_totpts
w_tot=fltarr(n_totpts)
i_tot=intarr(n_totpts)
;while loop to refine efficiencies, stop when they are no longer changing according to the tolerance
icount=0 & deviation=tolerance*2
dev_old=0. & dev_new=0. & del_dev=1.; variables to follow gradient of deviation, stop on sign change
while (icount lt max_cycle) and (abs(del_dev) gt 1.0e-5) do begin
icount=icount+1 & print, 'doing efficiencies ', icount
;calculate the total including efficiencies, initialise total and counter arrays
w_tot[*]=0.0 & i_tot[*]=0
for i=0,n-1 do begin ; loop thro scans
for j=0,nt-1 do begin ; loop thro tubes
alllowlim=nsteps*j & allhilim=alllowlim+nsteps-1 ;limits to define a single tube
totlowlim=nsteps*(i+j) & tothilim=totlowlim+nsteps-1 ;limits to define a single tube, including scan offset
for jh=0,nb-1 do begin ; loop thro height blocks
hmin=jh*hb & hmax=(jh+1)*hb-1
tot_temp=total(w_all_3d[i,alllowlim:allhilim,hmin:hmax],3) ; a vector nsteps long
if (mean(tot_temp) gt ave_min and mean(tot_temp) lt ave_max) then begin
i_tot[totlowlim:tothilim]=i_tot[totlowlim:tothilim]+hb ; count how many times adding into w_tot
w_tot[totlowlim:tothilim]=w_tot[totlowlim:tothilim]+tot_temp*effg[j,jh]
endif else begin
;print, 'scan,tube,height omitted',i,j,jh
endelse
endfor
endfor
endfor
;calculate the average
index=where(i_tot gt 0,count) & w_tot[index]=w_tot[index]/i_tot[index]
;calculate differences between individual scans and the average i.e. efficiencies per scan i
;averaged over the number of scan steps
for i=0,n-1 do begin ; loop over different scans
for j=0,nt-1 do begin ; loop over detector tubes
alllowlim=nsteps*j & allhilim=alllowlim+nsteps-1 ;limits to define a single tube
totlowlim=nsteps*(i+j) & tothilim=totlowlim+nsteps-1 ;limits to define a single tube, including scan offset
for jh=0,nb-1 do begin ; loop thro height blocks, nb is number of blocks height hb
hmin=jh*hb & hmax=(jh+1)*hb-1
totot_temp=total(w_all_3d[i,alllowlim:allhilim,hmin:hmax])
if (totot_temp gt 0.0) then begin
eff[i,j,jh]=total(w_tot[totlowlim:tothilim])/totot_temp
endif else begin
eff[i,j,jh]=-1.0
endelse
endfor
endfor
endfor
;calulate average/global efficiency for detectors not omitted
for j=0,nt-1 do begin
for jh=0,nb-1 do begin
index=where(eff[*,j] gt 0.,count)
effg[j,jh]=hb*total(eff[index,j,jh],1)/count
endfor
endfor
;calculate deviation from previous efficiencies
dev_new=total(abs(effg-effgold))/nt/hb & print,'deviation, icount ',dev_new, icount
del_dev=dev_old-dev_new & dev_old=dev_new & print, del_dev
;update previous efficiencies
effgold=effg
endwhile ;*** efficiencies refined in effg
ave = total(effg)/n_elements(effg) & print, 'ave efficiency ', ave
effg=effg/ave
;write effg into effout, the transpose
for j=0,nt-1 do begin
for jh=0,nb-1 do begin
hmin=jh*hb & hmax=(jh+1)*hb-1
effout[hmin:hmax,j]=effg[j,jh]
endfor
endfor
;normalise efficiencies
; ave=total(effout)/n_elements(effout) & effout=effout/ave & print, 'average eff is ',ave
; ****************** update calibration file ******************
;effout in format (nh,nt), write new calibration file
print, 'writing calibration file'
openr,f_in,cal_fil_in,/get_lun
openw,f_out,cal_fil_out,/get_lun
;write new zones info
printf,f_out,title + ' (method B)'
;copy zones from old calibration file
junk=''
readf,f_in,junk ;skip title
while junk ne 'efficiencies' do begin
readf,f_in,junk
printf,f_out,junk
endwhile
;write efficiencies
printf,f_out,effout
;read down to angles in old file
while junk ne 'angles' do readf,f_in,junk
;then write angles from old file
printf,f_out,'angles'
while junk ne 'calib_end' do begin
readf,f_in,junk
printf,f_out,junk
endwhile
;printf,f_out,'calib_end'
free_lun, f_in & free_lun, f_out
return, effout
end
FUNCTION calib_eff_2d_C, numor, n, wi, cal_fil_in, cal_fil_out, title
;+
;******* ****************
;**
;** User callable
;** The call is e.g. w2=calib_eff_2d(51577, 5, 20, 'd2bcal.raw', 'd2bcal.eff', 'title')
;**
;** 51577 is the first scan, 5 is the number of scans, 20 is the number of a working workspace
;** w2 is an output workspace containing the table of efficiencies
;** d2bcal.raw is the old file, d2bcal.eff is the new one
;** data must be read in with a calibration file that defines the active zones
;** since efficiencies are to be determined for these zones,
;** *** but the efficiencies table in ***
;** *** the input calibration file should contain only 1's ***
;** this macro as for calib_eff_tubes_1d, but efficiencies per vertical block (height hb) rather than per tube
;** in "1d" version integrate over height at the start
;** ---------------------------------------------------
;** METHOD
;** ------
;** Sum data using initial efficiencies to get average diffraction pattern
;** Difference between each scan and the average gives efficiency per tube, height group and scan
;** Average efficiency over different scans gives new effciency per tube, height group
;** Repeat
;** Currently the beamstop is not omitted.
;** For this, identify tubes that scan the beam stop
;** Use tube_score(scan,tube)
;** ---------------------------------------------------
;-
forward_function rdrun
;read the first numor and get array dimensions
datp=1
win=rdrun(numor,DATP=datp,w=wi)
print, n_elements(datp)
npts=(size(win))[1] & print, 'npoints ',npts
nh =(size(win))[2] & print, 'nh ',nh
nt=128 ;number of detector tubes
nsteps=npts/nt ;number of scan steps e.g.25, n is number of different scans
tolerance=1e-4 ;used to stop refinment cycle for efficiencies
max_cycle=30 ;max no. of cycles in efficiency calculation
min_cycle=3 ;min no. of cycles in efficiency calculation
hb=16 ;height of blocks for efficiency calculation
nb=nh/hb ;number of blocks in height
;define matrix of efficiencies per scan (eff) and global (effg & effgold) (new and old)
eff=fltarr(n,nt,nb) & effg=fltarr(nt,nb) & effgold=fltarr(nt,nb)
eff[*,*,*]=1.0 & effg[*,*]=1.0 & effgold[*,*]=2.0
effout =fltarr(nh,nt) & effout[*,*]=1.0
;make array to store all calibration workspaces and enter the first scan, read the rest
w_all_3d=fltarr(n,npts,nh)
w_all_3d[0,*,*]=win
print, 'reading data'
for i=1,n-1 do begin
newnumor=numor+i
win=rdrun(newnumor)
w_all_3d[i,*,*]=win
endfor
w_all_3d_blocks = fltarr(n,npts,nb)
for jh=0,nb-1 do begin ; loop thro height blocks
hmin=jh*hb & hmax=(jh+1)*hb-1
w_all_3d_blocks[*,*,jh] = total(w_all_3d[*,*,hmin:hmax],3)
endfor
median_pixel = median(w_all_3d)
median_block = median(w_all_3d_blocks)
median_tube = median(w_all_3d) * 128
counts_in_blocks = fltarr(nt,nb)
counter_in_blocks = fltarr(nt,nb)
tolerance = 0.5 * median_tube * nsteps
median_block_nsteps = median_block * nsteps
tubes_done = fltarr(nt)
adding = []
for i=0,n-1 do begin ; loop thro scans
for j=0,nt-1 do begin ; loop thro tubes
tube_lowlim=nsteps*j & tube_hilim=tube_lowlim+nsteps-1 ;limits to define a single tube
;only consider full tubes not masked by beamstop
if (total(w_all_3d_blocks[i, tube_lowlim:tube_hilim, *]) gt tolerance) then begin
tubes_done[j] = 1.0
for k=0,nb-1 do begin ; loop thro height blocks
adding = [adding, total(w_all_3d_blocks[i, tube_lowlim:tube_hilim, k])]
counts_in_blocks[j,k] += total(w_all_3d_blocks[i, tube_lowlim:tube_hilim, k])
counter_in_blocks[j,k] += 1
endfor
endif
endfor
endfor
eff_blocks = fltarr(nt,nb)*0.0 + 1.0
check = where(counter_in_blocks le 0.0, missing_blocks, complement=good_blocks)
if missing_blocks gt 0 then print, 'Problem: Efficiency could not be calculated from some detector blocks. Set to 1'
eff_blocks[good_blocks] = median_block_nsteps/(counts_in_blocks[good_blocks]/counter_in_blocks[good_blocks])
;write effg into effout, the transpose
for j=0,nt-1 do begin
for jh=0,nb-1 do begin
hmin=jh*hb & hmax=(jh+1)*hb-1
effout[hmin:hmax,j]=eff_blocks[j,jh]
endfor
endfor
; ****************** update calibration file ******************
;effout in format (nh,nt), write new calibration file
print, 'writing calibration file'
openr,f_in,cal_fil_in,/get_lun
openw,f_out,cal_fil_out,/get_lun
;write new zones info
printf,f_out,title + ' (method C)'
;copy zones from old calibration file
junk=''
readf,f_in,junk ;skip title
while junk ne 'efficiencies' do begin
readf,f_in,junk
printf,f_out,junk
endwhile
;write efficiencies
printf,f_out,effout
;read down to angles in old file
while junk ne 'angles' do readf,f_in,junk
;then write angles from old file
printf,f_out,'angles'
while junk ne 'calib_end' do begin
readf,f_in,junk
printf,f_out,junk
endwhile
;printf,f_out,'calib_end'
free_lun, f_in & free_lun, f_out
return, effout
end
FUNCTION calib_eff_2d_D, numor, n, wi, cal_fil_in, cal_fil_out, title
;+
;******* ****************
;**
;** User callable
;** The call is e.g. w2=calib_eff_2d(51577, 5, 20, 'd2bcal.raw', 'd2bcal.eff', 'title')
;**
;** 51577 is the first scan, 5 is the number of scans, 20 is the number of a working workspace
;** w2 is an output workspace containing the table of efficiencies
;** d2bcal.raw is the old file, d2bcal.eff is the new one
;** data must be read in with a calibration file that defines the active zones
;** since efficiencies are to be determined for these zones,
;** *** but the efficiencies table in ***
;** *** the input calibration file should contain only 1's ***
;** this macro as for calib_eff_tubes_1d, but efficiencies per vertical block (height hb) rather than per tube
;** in "1d" version integrate over height at the start
;** ---------------------------------------------------
;** METHOD
;** ------
;** Sum data using initial efficiencies to get average diffraction pattern
;** Difference between each scan and the average gives efficiency per tube, height group and scan
;** Average efficiency over different scans gives new effciency per tube, height group
;** Repeat
;** Currently the beamstop is not omitted.
;** For this, identify tubes that scan the beam stop
;** Use tube_score(scan,tube)
;** ---------------------------------------------------
;-
forward_function rdrun
;read the first numor and get array dimensions
datp=1
win=rdrun(numor,DATP=datp,w=wi)
print, n_elements(datp)
npts=(size(win))[1] & print, 'npoints ',npts
nh =(size(win))[2] & print, 'nh ',nh
nt=128 ;number of detector tubes
nsteps=npts/nt ;number of scan steps e.g.25, n is number of different scans
tolerance=1e-4 ;used to stop refinment cycle for efficiencies
max_cycle=30 ;max no. of cycles in efficiency calculation
min_cycle=3 ;min no. of cycles in efficiency calculation
hb=1 ;height of blocks for efficiency calculation
nb=nh/hb ;number of blocks in height
;define matrix of efficiencies per scan (eff) and global (effg & effgold) (new and old)
eff=fltarr(n,nt,nb) & effg=fltarr(nt,nb) & effgold=fltarr(nt,nb)
eff[*,*,*]=1.0 & effg[*,*]=1.0 & effgold[*,*]=2.0
effout =fltarr(nh,nt) & effout[*,*]=1.0
;make array to store all calibration workspaces and enter the first scan, read the rest
w_all_3d=fltarr(n,npts,nh)
w_all_3d[0,*,*]=win
print, 'reading data'
for i=1,n-1 do begin
newnumor=numor+i
win=rdrun(newnumor)
w_all_3d[i,*,*]=win
endfor
print, 'data read!'