Commit 1ee2eb71 authored by Miguel Angel Gonzalez's avatar Miguel Angel Gonzalez
Browse files

Issue 31: Some changes to the way or reading/calibrating D2B files depending...

Issue 31: Some changes to the way or reading/calibrating D2B files depending on the mode, raw vs no raw
parent cb090600
......@@ -2,7 +2,7 @@ PRO rdid_d2b, INST,numor,nvers,text,exper,scan,cnt,nd,WOUT,vparm,param,par1,par2
WT ,XT ,YT ,ZT ,OT ,DATE ,PP ,PTXT ,XX ,YY ,ZZ ,NN ,PV ,EE, raw
;+
;** ********
;** read D2B ascii files, should be un_used -> see rdid
;** read D2B ascii files
;-
common calibration, pathcal, cal_d19 , cal_d2b , ang_d2b , cal_d1a , ang_d1a , cal_in13 $
......@@ -13,152 +13,179 @@ PRO rdid_d2b, INST,numor,nvers,text,exper,scan,cnt,nd,WOUT,vparm,param,par1,par2
common c_rdid , dzap, pzap, pzip ,pzup
; *** data will be (number of detectors, height, number of points) ***
; *** WOUT will be 2D (nt*nj,nh), ZZ contains angles ***
; *** for calibration ... ***
; June 2019:
; The way of reading/handling the data has changed following 2018 discussions and new
; measurements on how to calibrate properly the data. The main difference is that now
; the calibration is done with Mantid routines, which compute an efficiency for each
; pixel (including "non active" zones), and determine the true pixel height (not used
; here) and the tube center (not used here). So now, when working in 'raw' mode, the
; raw counts are read, multiplied by the efficiencies and the scan reformed directly.
; On the other hand, before 2018, the 1st step of the calibration consisted in determining
; the top and bottom pixels to be discarded, then 'expanding' (using CONGRID) the
; remaining pixels to 128 (to keep a constant number of pixels per tube), and then
; computing the efficiencies of the 'renormalized pixels'. This behaviour can be
; recovered with the 'no raw' option, so typically one should:
;
; 1. Use 'raw' to read D2B files when using a calibration file generated by Mantid,
; and then discard the required number of pixels on the top and bottom of the detector.
;
; 2. Use 'no raw' when using a calibration file generated by Lamp, using only an
; active detector zone and expanding the remaining pixels to 128.
;
; A modification introduced on 20/6/2019 is that the option to remove some selected
; tubes (given by using the command dud, dets=[] before reading the data) is applied
; now in both modes (before it was inside the IF block corresponding to the 'no raw'
; mode).
; *** data will be (number of detectors, height, number of points) ***
; *** WOUT will be 2D (nt*nj,nh), ZZ contains angles ***
; *** for calibration ... ***
;debug = 0 --> no print, >0 --> print messages
debug = 0
; *** set physical height of detector (30cm) ***
; *** set physical height of detector (30cm) ***
h=30.0
D2TH= param[36]
nj = ROUND(n_elements(WOUT)/(nd)) ; number of points in scan
nt = 128 ; number of detector tubes
nh = nd/nt ; number of height pixels for 128 detector tubes
D2TH = param[36]
nj = ROUND(n_elements(WOUT)/(nd)) ; number of points in scan
nt = 128 ; number of detector tubes
nh = nd/nt ; number of height pixels for 128 detector tubes
; *** rewrite WOUT in matrix form (nh, nt, nj) ***
WOUT= reform(WOUT, nh, nt ,nj ,/overwrite)
; NN = WOUT(0,0) ; set monitor
WOUT = reform(WOUT, nh, nt ,nj ,/overwrite)
YT ='height'
XT ='2*Theta'
YT = 'height'
XT = '2*Theta'
; *** definition of new calibration arrays (default: equal to old)
newang_d2b=ang_d2b
; *** definition of new calibration arrays (default: equal to old)
newang_d2b = ang_d2b
print, 'Linear interpolation of the tubes for ',numor
; *** now we have risers and falls for the detector ***
; *** interpolate data onto 128 (nh) channels ***
; *** reverse the data in the second tube ***
IF raw NE '1' THEN BEGIN
print, 'Linear interpolation of the tubes for ', numor
IF raw NE '1' THEN BEGIN ;first and last useful pixels given in calibration file
rise1=zon_d2b[0,*] & fall1=zon_d2b[1,*] & rise2=zon_d2b[2,*] & fall2=zon_d2b[3,*]
ENDIF ELSE BEGIN
rise1=intarr(nt/2) & fall1=rise1+127 & rise2=rise1+128 & fall2=rise1+255 ; default values - commented!!!
ENDIF ELSE BEGIN ;raw mode: use full tube
rise1=intarr(nt/2) & fall1=rise1+127 & rise2=rise1+128 & fall2=rise1+255
ENDELSE
; *** now we have risers and falls for the detector ***
; *** interpolate data onto 128 (nh) channels ***
; *** reverse the data in the second tube ***
FOR j=0,nj-1 DO BEGIN ; loop thro scan points nj
FOR i=0,nt/2-1 DO BEGIN ;
WONE2 = [WOUT[*,i*2,j],WOUT[*,i*2+1,j]] ; put pairs of tubes together again
WOUT[*,i*2, j]=CONGRID(WONE2[rise1[i]:fall1[i]],nh,/INTERP)
WOUT[*,i*2+1,j]=REVERSE(CONGRID(WONE2[rise2[i]:fall2[i]],nh,/INTERP))
WONE2 = [WOUT[*,i*2,j], WOUT[*,i*2+1,j]] ; put pairs of tubes together again
WOUT[*,i*2, j] = CONGRID(WONE2[rise1[i]:fall1[i]], nh, /INTERP)
WOUT[*,i*2+1,j] = REVERSE(CONGRID(WONE2[rise2[i]:fall2[i]], nh, /INTERP))
ENDFOR
ENDFOR
; *** if calibration available, calibrate by MULTIPLYING ***
IF n_elements(cal_d2b) eq nh*nt THEN BEGIN ;if calibration file exists
if debug ge 1 then print, 'CALIBRATING'
WOUT=float(WOUT)
FOR i=0,nj-1 DO WOUT[*,*,i]=WOUT[*,*,i]*cal_d2b ;newcal has dimension nh*nt i.e. resolution per pixel
fct =1. ;factor set to 1 to preserve angles in ang_d2b
OT =OT+' /'+inf_d2b[0] ; set sub-title
; *** fct*newang_d2b always nd values spaced by ~ -1.25 ***
ENDIF ELSE BEGIN ;calibration does not exist, generate angles
if n_elements(newang_d2b) ne nt then newang_d2b=127.0-findgen(nt)
fct = -1.25
; *** if calibration available, calibrate by MULTIPLYING ***
; *** fct*newang_d2b always nd values spaced by ~ -1.25 ***
IF n_elements(cal_d2b) eq nh*nt THEN BEGIN
if debug ge 1 then print, 'CALIBRATING'
WOUT = float(WOUT)
FOR i=0,nj-1 DO WOUT[*,*,i] = WOUT[*,*,i]*cal_d2b ;newcal has dimension nh*nt i.e. resolution per pixel
fct = 1. ;factor set to 1 to preserve angles in ang_d2b
OT = OT + ' /' + inf_d2b[0] ; set sub-title
ENDIF ELSE BEGIN
;calibration does not exist, generate angles
if n_elements(newang_d2b) ne nt then newang_d2b = 127.0-findgen(nt)
fct = -1.25
ENDELSE
; *** find and remove dud detectors
; *** find detectors with efficiency too different from 1, work with detector number=index+1
; *** find and remove dud detectors
zap = dzap[uniq(dzap, sort(dzap))]
idz = where(zap ge 1 and zap le 128) ; *** test if zap contains real detectors **
if idz[0] ne -1 then zap=zap[idz] else zap = 0 ; *** in case, eliminate spurious values ***
SD = SIZE(zap)
IF SD[0] GT 0 THEN BEGIN
print, 'Dud detectors:'
WOUT[*, zap-1, *] = -999 ; *** set all dud values and angles to -999 ***
newang_d2b[zap-1] = -999
GOOD = where(WOUT ne -999) ; *** good values are those which do not contain -999 ***
WOUT = WOUT[GOOD]
newang_d2b = newang_d2b[where(newang_d2b ne -999)]
; *** redefine the number of detectors
nt = (nt-n_elements(zap)) & print, n_elements(zap),' detectors removed.'
; *** reform WOUT to proper format
WOUT = reform(WOUT, nh, nt, nj ,/overwrite)
ENDIF
IF raw NE '1' THEN BEGIN
if n_elements(pzip) ne 1 then pzip=0
if pzip gt 0 then seuil=pzip else seuil=3.0
if debug ge 4 then print,'cutoff=',seuil
ezap=0
if n_elements(pzip) ne 1 then pzip = 0
if pzip gt 0 then seuil = pzip else seuil = 3.0
if debug ge 4 then print, 'cutoff=', seuil
ezap = 0
; *** find detectors with efficiency too different from 1, work with detector number=index+1
for i=0,nt-1 do begin
cal_1d=cal_d2b[*,i]
maxeff=max(cal_1d,min=mineff) ; find limiting efficiencies over height
cal_1d = cal_d2b[*,i]
maxeff = max(cal_1d, min=mineff) ; find limiting efficiencies over height
if (maxeff gt seuil) or (mineff lt 1/seuil) then ezap=[ezap,i+1] ; compare with seuil
endfor
; *** combine dzap and ezap
zap=[dzap,ezap]
zap=zap[uniq(zap,sort(zap))]
; *** test if zap contains real detectors **
idz=where(zap ge 1 and zap le 128)
; *** in case, eliminate spurious values ***
if idz[0] ne -1 then zap=zap[idz] else zap=0
; *** is there something left? ***
SD=SIZE(zap)
; *** if zap has something, remove corresponding detector tubes from WOUT ***
zap = ezap[uniq(ezap,sort(ezap))]
idz = where(zap ge 1 and zap le 128) ; *** test if zap contains real detectors **
if idz[0] ne -1 then zap = zap[idz] else zap = 0 ; *** in case, eliminate spurious values ***
SD = SIZE(zap)
IF SD[0] GT 0 THEN BEGIN
print, 'Dud detectors:'
; *** set all dud values and angles to -999 ***
WOUT[*,zap-1,*]=-999
newang_d2b[zap-1]=-999
; *** good values are those which do not contain -999 ***
GOOD=where(WOUT ne -999)
; if GOOD(0) eq -1 then return, GOOD(0)
; *** cut out dud values from WOUT and angles, not calibration as already done ***
WOUT=WOUT[GOOD]
newang_d2b=ang_d2b[where(newang_d2b ne -999)]
print, 'Too high/low efficiencies:'
WOUT[*,zap-1,*] = -999 ; *** set all dud values and angles to -999 ***
newang_d2b[zap-1] = -999
GOOD = where(WOUT ne -999) ; *** good values are those which do not contain -999 ***
WOUT = WOUT[GOOD]
newang_d2b = newang_d2b[where(newang_d2b ne -999)]
; *** redefine the number of detectors
nt=(nt-n_elements(zap)) & print,n_elements(zap),' detectors removed.'
nt = (nt-n_elements(zap)) & print, n_elements(zap), ' detectors removed.'
; *** reform WOUT to proper format
WOUT=reform(WOUT, nh, nt, nj ,/overwrite)
WOUT = reform(WOUT, nh, nt, nj ,/overwrite)
ENDIF
ENDIF
; *** set tubes positions for first detector position
; newang_d2b=nt-1-newang_d2b ;turn angles around
if nj gt 1 then XX=fct*newang_d2b+ZZ[0] else XX=1.25*findgen(nh) ; angle not picked up for one point scan
if nj gt 1 then XX = fct*newang_d2b+ZZ[0] else XX = 1.25*findgen(nh) ; angle not picked up for one point scan
; *** reform WOUT to be 2D (nt*nj,nh), transpose first to put nt&nj together ***
if nj gt 1 then WOUT= transpose(WOUT,[1,2,0]) else WOUT= transpose(WOUT,[1,0]) ; WOUT looses 3rd dimension in one point scan
if debug ge 1 then print, 'Concatenating detector positions '
WOUT= reform(WOUT,nt*nj,nh,/overwrite)
; *** do the following if there are at least 2 points in scan ***
; *** create a XX array by adding up all 2th for all dets **
; *** in the same order as in WOUT
WOUT = reform(WOUT,nt*nj, nh, /overwrite)
; *** do the following if there are at least 2 points in scan ***
; *** create a XX array by adding up all 2th for all dets **
; *** in the same order as in WOUT
IF nj gt 1 then BEGIN ; i.e. not single point scan
FOR i=1,nj-1 do XX=[XX , fct*newang_d2b+ZZ[i]]
; *** sort the XX and WOUT arrays, preserve 2D of latter (idx is the index array) ***
idx =sort( XX) & XX =XX[idx] & WOUT=WOUT[idx,*]
Paolo=1 & Did=0
if Paolo then begin
; *** XR contains the 2th's rounded to the nearest step size ***
XR =round(XX/D2TH)*D2TH
; *** GRID contains the fractional indices of XR into XX ***
GRID=(XR-XX)/D2TH+findgen(n_elements(XX))
; *** for each height do a 1D interpolation using WOUT_1D and reassign WOUT***
FOR i=0,nh-1 DO WOUT[*,i]=interpolate(WOUT[*,i],GRID,CUBIC=-0.5)
; *** reassign XX
XX=XR
endif else if Did then begin
FX=round(XX[0]/D2TH)*D2TH & LX=round(XX[n_elements(XX)-1]/D2TH)*D2TH
Un=fix((LX-FX)/D2TH+1)
UX=findgen(Un)*D2TH+FX ; Regular grid without holes and doublons !!!
WTMP=fltarr(Un,nh)
FOR i=0,nh-1 DO WTMP[*,i]=interpol(WOUT[*,i],XX,UX)
XX=UX & WOUT=WTMP
endif
; *** YY is the height, EE the error
YY = findgen(nh)*h/(nh-1) - h/2.0
ZZ = param[46]
WOUT=WOUT>0.
EE = SQRT(WOUT)
FOR i=1,nj-1 do XX = [XX , fct*newang_d2b+ZZ[i]]
; *** sort the XX and WOUT arrays, preserve 2D of latter (idx is the index array) ***
idx = sort(XX) & XX =XX[idx] & WOUT = WOUT[idx,*]
Paolo = 1 & Did=0
if Paolo then begin
XR = round(XX/D2TH)*D2TH ; XR contains the 2th's rounded to the nearest step size ***
GRID = (XR-XX)/D2TH + findgen(n_elements(XX)) ; GRID contains the fractional indices of XR into XX ***
FOR i = 0,nh-1 DO WOUT[*,i] = interpolate(WOUT[*,i], GRID, CUBIC=-0.5) ; for each height do a 1D interpolation using WOUT_1D and reassign WOUT***
XX = XR ; reassign XX
endif else if Did then begin
FX = round(XX[0]/D2TH)*D2TH & LX = round(XX[n_elements(XX)-1]/D2TH)*D2TH
Un = fix((LX-FX)/D2TH+1)
UX = findgen(Un)*D2TH + FX ; Regular grid without holes and doublons !!!
WTMP = fltarr(Un,nh)
FOR i=0,nh-1 DO WTMP[*,i] = interpol(WOUT[*,i],XX,UX)
XX = UX & WOUT = WTMP
endif
YY = findgen(nh)*h/(nh-1) - h/2.0 ; detector height
ZZ = param[46]
WOUT = WOUT>0.
EE = SQRT(WOUT)
ENDIF
END
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment