Commit 536bd012 authored by Miguel Angel Gonzalez's avatar Miguel Angel Gonzalez
Browse files

Issue 54. Update Cosmos to handle new NeXus versions

parent 4ed839b3
;TS 2021-04-27 Bug fixes
;TS 2021-04-20 Implementation of recent NEXUS file changes
;TS 2021-03-31 Corrected TOF distance for interchopper distance
;TS 2021-02-25 + 03-26 added reading of pcen, chopsep,chopwin,cr from datafile
;TS 2021-02-25 Corrected Pcen from 106 to 132.5
;TS 2021-02-20 Got rid of some of the many print outputs
;TS 2021-02-17 Needed to correct reading of flipper states from Nexus file
;TS 2020-02-15 added analyzer to san dan options in coherent and incoherent case
;TS 2019-01-11 Added a print output for upper and lower end of the wavelength of each time channel
;TS 2018-06-04 Removed the restriction that om_fwhm is disregarded for polarized runs since the resolution calculation now should be based only on the UP polarization.
......@@ -198,28 +205,19 @@ pro cosmos_anal_backgroundaverage, arr, err, maskbk, str
backgmax = findgen(xsize) * gradmax + gradmax
backgmin = findgen(xsize) * gradmin + gradmin
backgerr=0.5*abs(backgmax-backgmin)
;if t eq 100 or t eq 500 or t eq 600 or t eq 700 or t eq 800 then begin
;print,(t),gradoff(0,t),gradoff(1,t)
;print,arr(maskbk,t)
;endif
backg = findgen(xsize) * gradoff[0,t] + gradoff[1,t] ; 1D array in x
arr[*, t] = arr[*, t] - backg
;if t eq 100 or t eq 500 or t eq 600 or t eq 700 or t eq 800 then begin
;print,arr(maskbk,t)
;print,backg(maskbk)
;endif
; add error in quadrature
err[*, t] = sqrt(err[*, t]^2 +backgerr^2)
;cosmos_logmessage, ' Background average slope to ' + str + ' runs with regression factors y = a x + b :'
; cosmos_logmessage, string(9b) + 'min' + string([9b, 9b]) + 'max'
; cosmos_logmessage, ' a' + string(9b) + cosmos_str_make(min(gradoff[1, *])) + string(9b) + cosmos_str_make(max(gradoff[1, *]))
; cosmos_logmessage, ' b' + string(9b) + cosmos_str_make(min(gradoff[0, *])) + string(9b) + cosmos_str_make(max(gradoff[0, *]))
endif else begin
;only one background strip so must subtract a constant -the eman of that background strip
bk=mean(arr[maskbk,t])
arr[*, t] = arr[*, t] - bk
bkgnderrorsquared = total(err[maskbk, t]^2) / (msize^2)
err[*, t] = sqrt(err[*, t]^2 + bkgnderrorsquared)
;only one background strip so must subtract a constant -the eman of that background strip
bk=mean(arr[maskbk,t])
arr[*, t] = arr[*, t] - bk
bkgnderrorsquared = total(err[maskbk, t]^2) / (msize^2)
err[*, t] = sqrt(err[*, t]^2 + bkgnderrorsquared)
endelse
endfor
......@@ -227,6 +225,7 @@ pro cosmos_anal_backgroundaverage, arr, err, maskbk, str
return
end
; subtract fitted background
pro cosmos_anal_backgroundfit, arr, err, maskbk, str
......@@ -247,10 +246,9 @@ pro cosmos_anal_backgroundfit, arr, err, maskbk, str
backgmin = findgen(xsize) * (fit[1,i]-sigz[1, i]) + fit[0,i]+sigz[0, i] ; 1D array in x
backgerr = 0.5*abs(backgmax-backgmin)
;print,i,fit(1,i),sigz(1,i),fit(0,i),sigz(0,i)
arr[*, i] = arr[*, i] - backg
; add error in quadrature
err[*, i] = sqrt(err[*, i]^2 + backgerr^2)
endfor
......@@ -518,7 +516,7 @@ end
;-
pro cosmos_tofpolcorr_db, array, array_err, wavelength
cosmos_logmessage, 'Direct beam polarization correction done.', /warning
cosmos_logmessage, 'Direct beam polarization correction done.' ;TS 2021-04-27 Removed the output as warning
cosmos_tofpolcorr_factors, wavelength, F1L, F2L, PhiL, P1L, P2L, eF1L, eF2L, ePhiL, eP1L, eP2L, InMat
;* Calculate the corrected counts & error (do first errors because original array needed)
......@@ -534,7 +532,7 @@ end
function cosmos_tofpolcorr_noana, array
cosmos_logmessage, 'Performing polarization correction without analyzer...', /warning
cosmos_logmessage, 'Performing polarization correction without analyzer...' ;TS 2021-04-27 Removed the output as warning
s = size(array)
if s[0] ne 3 and s[0] ne 4 then begin
......@@ -1110,106 +1108,6 @@ function cosmos_anal_factor, array_a, array_b
end
; new grouping procedure - groups integer number of q points into each new bin
; modified 07/05 to group integral numbers of q points subject to constant multiplier of instrument resolution (given by groupnum)
; MG 28.01.2015 --> average over number of grouped points (needed to treat both coherent/incoherent cases)
; MG 6.2.2015 --> routine rewritten because BC and PG arrived to the conclusion that the Q resolution was not OK
;function cosmos_anal_group, result, groupnum
;
; ; display min and max fractional resolutions for dataset
; fractionalerror = result[1, *] / result[0, *]
; minfracerror = min(fractionalerror)
; maxfracerror = max(fractionalerror)
;
; if (groupnum le 0.0) then begin
;
; groupnum = 1.0
; cosmos_logmessage, ' Grouping data to calculated Q resolution of instrument = (' + cosmos_str_make(minfracerror * 100.0) + '%, ' + cosmos_str_make(maxfracerror * 100.0) + '%) (Auto)'
;
; endif else begin
;
; cosmos_logmessage, ' Grouping data to ' + cosmos_str_make(groupnum) + ' times calculated Q resolution of instrument = (' + cosmos_str_make(minfracerror * 100.0) + $
; '%, ' + cosmos_str_make(maxfracerror * 100.0) + '%)'
;
; endelse
;
; i = 0
; newpt = 0
; sortedresult = cosmos_anal_sort(result)
; groupedresult = sortedresult ; preserve unbundled data!
;
; j = 0 ; number of points in new bin
; minj = 0;
; maxj = 0;
;
; for i = 0, n_elements(result[0, *]) - 1 do begin
;
; if (j eq 0) then begin
;
; ; groupedresult[0, *] are min values
; groupedresult[0, newpt] = sortedresult[0, i] - sortedresult[1, i]
; ; groupedresult[1, *] are max values
; groupedresult[1, newpt] = sortedresult[0, i] + sortedresult[1, i]
; groupedresult[2:4, newpt] = sortedresult[2:4, i]
; currentthreshold = sortedresult[0, i] + (sortedresult[1, i] * groupnum)
; j = j + 1
;
; endif else begin
;
; if (sortedresult[0, i] ge currentthreshold) then begin
;
; if (maxj eq 0) then begin
; minj = j
; maxj = j
; endif else begin
; if (j lt minj) then minj = j
; if (j gt maxj) then maxj = j
; endelse
;
; ;average before going to next point
; if j gt 0 then groupedresult[2:4, newpt] = groupedresult[2:4, newpt] / float(j)
;
; newpt = newpt + 1
; groupedresult[0, newpt] = sortedresult[0, i] - sortedresult[1, i]
; groupedresult[1, newpt] = sortedresult[0, i] + sortedresult[1, i]
; groupedresult[2:4, newpt] = sortedresult[2:4, i]
; currentthreshold = sortedresult[0, i] + (sortedresult[1, i] * groupnum)
; j = 1
;
; endif else begin
;
; groupedresult[0, newpt] = groupedresult[0, newpt] < (sortedresult[0, i] - sortedresult[1, i])
; groupedresult[1, newpt] = groupedresult[1, newpt] > (sortedresult[0, i] + sortedresult[1, i])
; groupedresult[2, newpt] = groupedresult[2, newpt] + sortedresult[2, i]
; groupedresult[3, newpt] = sqrt(groupedresult[3, newpt]^2 + sortedresult[3, i]^2)
; groupedresult[4, newpt] = groupedresult[4, newpt] + sortedresult[4, i]
; j = j + 1
;
; endelse
;
; endelse
;
; endfor
;
; ;average for last point
; if j gt 0 then groupedresult[2:4, newpt] = groupedresult[2:4, newpt] / float(j)
;
; groupedresult = groupedresult[*, 0 : newpt]
;
; ; turn min, max back into value, error
; for i = 0, n_elements(groupedresult[0, *]) - 1 do begin
;
; groupedresult[1, i] = (groupedresult[1, i] - groupedresult[0, i]) / 2.
; groupedresult[0, i] = groupedresult[0, i] + groupedresult[1, i]
;
; endfor
;
; cosmos_logmessage, ' Min/max numbers of old points in new bins = ' + cosmos_str_make(minj) + ', ' + cosmos_str_make(maxj)
;
; return, groupedresult
;
;end
function cosmos_anal_group, result, groupnum
; display min and max fractional resolutions for dataset
......@@ -1253,11 +1151,6 @@ function cosmos_anal_group, result, groupnum
;PG: I would use the absolute value, not the square
;weights = abs(sortedresult[2,points] / sortedresult[3,points])
weights[*] = 1 ;PG: actually, maybe it shouldn't be weighted at all
;print, sortedresult[0, i]
;print, old_threshold, new_threshold
;help, q_values, weights
;print, q_values
;print, weights
;average Q (weighted by errors)
groupedresult[0, newpt] = total(weights*q_values) / total(weights) ;PG=Maybe the final q shouldn't be weighted
......@@ -2576,7 +2469,7 @@ pro cosmos_anal_run, result, groupedresult, direct, reflect, instr, theta, runno
angle_bragg = angle_centre + $
0.5*atan((newpeakdir-c_params.pcen)*c_params.pixelwidth/pardir.sdetd) - $
0.5*atan(((float(peakref)+0.5)-c_params.pcen)*c_params.pixelwidth/parref.sdetd)
cosmos_logmessage, 'DAN coherent option NO ANALYZER = !!!!!' + cosmos_str_make(tra) , /WARNING
;cosmos_logmessage, 'DAN coherent option NO ANALYZER = !!!!!' + cosmos_str_make(tra) , /WARNING
endif else begin ;TS 2020-02-15 with analyzer
angle_bragg = angle_centre - $
0.5*atan((newpeakdir-c_params.pcen)*c_params.pixelwidth/pardir.sdetd) + $
......@@ -2592,12 +2485,12 @@ pro cosmos_anal_run, result, groupedresult, direct, reflect, instr, theta, runno
angle_bragg = angle_centre + $
0.5*atan((newpeakdir-c_params.pcen)*c_params.pixelwidth/pardir.sdetd) - $
0.5*atan((newpeakref-c_params.pcen)*c_params.pixelwidth/parref.sdetd)
cosmos_logmessage, 'DAN incoherent option NO ANALYZER = !!!!!' + cosmos_str_make(tra) , /WARNING
cosmos_logmessage, 'DAN incoherent option' ; TS 2021-04-27 Removed the warning
endif else begin ;TS 2020-02-15 with analyzer
angle_bragg = angle_centre - $
0.5*atan((newpeakdir-c_params.pcen)*c_params.pixelwidth/pardir.sdetd) + $
0.5*atan((newpeakref-c_params.pcen)*c_params.pixelwidth/parref.sdetd)
cosmos_logmessage, 'DAN incoherent option WITH ANALYZER = !!!!!' + cosmos_str_make(tra) , /WARNING
cosmos_logmessage, 'DAN incoherent option WITH ANALYZER = !!!!!' + cosmos_str_make(tra) ;TS 2021-04-27 Removed the warning
endelse
endelse
......@@ -2711,7 +2604,7 @@ pro cosmos_anal_run, result, groupedresult, direct, reflect, instr, theta, runno
angle_bragg = (theta * !pi / 180.) - $
0.5*atan((float(peakref)-c_params.pcen)*c_params.pixelwidth/parref.sdetd)+ $
0.5*atan((newpeakref-c_params.pcen)*c_params.pixelwidth/parref.sdetd)
cosmos_logmessage, 'Table option coherent NO ANALYZER = !!!!!' + cosmos_str_make(tra) , /WARNING
;cosmos_logmessage, 'Table option coherent NO ANALYZER = !!!!!' + cosmos_str_make(tra) , /WARNING
endif else begin ;TS 2020-02-15 With analyzer option for coherent binning. Need to mirror the reflection.
angle_bragg = (parref.san * !pi / 180.) + $
0.5*atan((float(peakref)-c_params.pcen)*c_params.pixelwidth/parref.sdetd) - $
......@@ -2751,9 +2644,6 @@ pro cosmos_anal_run, result, groupedresult, direct, reflect, instr, theta, runno
;TS 2019-01-11 Lower and upper bound of wavelength for each time channel.
;lambdaLowerbound = 1e10 * (c_params.planckperkg * ((findgen(tsize)) * parref.channelwidth + parref.delay) / ref_total_tofd)
;lambdaUpperbound = 1e10 * (c_params.planckperkg * ((findgen(tsize)+1) * parref.channelwidth + parref.delay) / ref_total_tofd)
;print,'array', lambda
;print,'lambdaLowerbound', lambdaLowerbound
;print,'lambdaUpperbound', lambdaUpperbound
cosmos_logmessage, ' '
cosmos_logmessage, ' Constants used in calculating lambda :'
......@@ -2831,7 +2721,7 @@ pro cosmos_anal_run, result, groupedresult, direct, reflect, instr, theta, runno
cosmos_tofpolcorr_db, cntdiri, cntdiri_err, lambda
endelse
endif else begin
cosmos_logmessage, 'Flipper state does not correspond to 00. Please use appropriate run for DB00 correction', /WARNING
cosmos_logmessage, 'DB Flipper state does not correspond to 00. State is :' + cosmos_str_make(pardir.f1) + cosmos_str_make(pardir.f2), /WARNING ;TS 2021-02-17 Output to identify problem.
endelse
endif ;endif polarized DB
endif
......@@ -2859,24 +2749,12 @@ pro cosmos_anal_run, result, groupedresult, direct, reflect, instr, theta, runno
;PG on 16/2/15: This had to be moved up, such that the information whether the sample is bent exists for Bob's procedure which follows
ratio = parref.s3w/parref.s2w
;coll = c_params.interslit
;sdd =pardir.sdetd;
;sdr =parref.sdetd;
sdr = parref.slit3 + parref.sdetd ;s3 to detector dist in m
vs = sdr+(ratio*c_params.interslit*0.001)/(1+ratio) ;virtual source position for the collimation [bc]
;print,'vs=',vs
da = 0.68 * sqrt((parref.s2w^2 + parref.s3w^2)/c_params.interslit^2) ;fwhm incoming angular width
s2_fwhm = (0.68*parref.s2w)/c_params.interslit ;fwhm angles of first slit PG:on 13/2/15 corrected for interslit distance, this has to be anticipated to choose the case for coherent COSMOS
detres = (2.2)*0.001 ;FWHM of detector resolution assuming a gaussian (m) same for both figaro and d17
da_det = sqrt((da*vs)^2+detres^2) ;fwhm on detector of da from virtual source in m
print,'da_det=',da_det
print,'det_fwhm=',det_fwhm
print,'detdb_fwhm=',detdb_fwhm
print, 'pixelwidth=', c_params.pixelwidth
; (CARE: should be an instrument parameter, CHECK: same resolution for Figaro/D17? BioRef in future?)
;MAG: Patch to make nodirect option to work, but this has to be corrected
;I guess coherent method makes no sense without DB - need to be able to calculate resol without DB
;To revise when last version of resolution calculation is defined by Bob and Philipp
if c_params.om_fwhm lt 0 then begin
......@@ -2885,7 +2763,6 @@ pro cosmos_anal_run, result, groupedresult, direct, reflect, instr, theta, runno
om_fwhm_real = 0.5*sqrt((det_fwhm)^2-da_det^2)/pardir.sdetd
if sqrt((det_fwhm)^2-da_det^2) ge c_params.pixelwidth then om_fwhm = 0.5*sqrt((det_fwhm)^2-da_det^2)/pardir.sdetd ;PG on 27/7/15: I think the threshhold should be the pixel size (4x)
if sqrt((det_fwhm)^2-da_det^2) lt c_params.pixelwidth then om_fwhm = 0
print, 'firstif=', sqrt((det_fwhm)^2-da_det^2)
endif else begin
om_fwhm_real = 0
om_fwhm = 0
......@@ -2896,7 +2773,6 @@ pro cosmos_anal_run, result, groupedresult, direct, reflect, instr, theta, runno
om_fwhm_real = 0.5*sqrt((det_fwhm)^2-detdb_fwhm^2)/pardir.sdetd
if sqrt((det_fwhm)^2-detdb_fwhm^2) ge c_params.pixelwidth then om_fwhm = 0.5*sqrt((det_fwhm)^2-detdb_fwhm^2)/pardir.sdetd
if sqrt((det_fwhm)^2-detdb_fwhm^2) lt c_params.pixelwidth then om_fwhm = 0
print, 'secondif=', sqrt((det_fwhm)^2-detdb_fwhm^2)
endif else begin
om_fwhm_real = 0
om_fwhm = 0
......@@ -2915,8 +2791,7 @@ pro cosmos_anal_run, result, groupedresult, direct, reflect, instr, theta, runno
endelse
twotheta = 2.0 * angle_bragg ;2theta of the specular
print, 'om_fwhm (deg)=', om_fwhm*180/!pi
print, 'om_fwhm_real (deg)=', om_fwhm_real*180/!pi
;Coherent treatment (only for normal_run)
if normal_run and calc_method eq 'coherent' then begin
......@@ -2942,11 +2817,11 @@ pro cosmos_anal_run, result, groupedresult, direct, reflect, instr, theta, runno
if tra gt 250 then begin ;;TS 2020-02-15 NO analyzer option for coherent binning. Standard.
tth_max = twotheta + float(peakref-fgdref[0]) * b_tth ; 2th foreground max
tth_min = twotheta + float(peakref-fgdref[1]) * b_tth ; 2th forground min
cosmos_logmessage, ' Coherent binning NO ANALYZER = !!!!!' + cosmos_str_make(tra)
;cosmos_logmessage, ' Coherent binning NO ANALYZER = !!!!!' + cosmos_str_make(tra)
endif else begin ;TS 2020-02-15 with analyzer option for coherent binning. Need to mirror the reflection.
tth_max = twotheta + float(fgdref[1]-peakref) * b_tth ; 2th foreground max
tth_min = twotheta + float(fgdref[0]-peakref) * b_tth ; 2th forground min
cosmos_logmessage, ' Coherent binning WITH ANALYZER = !!!!!' + cosmos_str_make(tra)
;cosmos_logmessage, ' Coherent binning WITH ANALYZER = !!!!!' + cosmos_str_make(tra) ;TS 2021-04-27
endelse
endelse
......@@ -2954,19 +2829,15 @@ pro cosmos_anal_run, result, groupedresult, direct, reflect, instr, theta, runno
;TS 2015-06-01 ((detres/(2*parref.sdetd)) changed to ((detres/(parref.sdetd))
;since the comparison to S2 does not need the factor of 2 in the detector resolution
if ((om_fwhm gt 0) AND ((detres/(parref.sdetd)) gt s2_fwhm)) then begin ;bent case
cosmos_logmessage, 'Bent case Lam max. ', /WARNING
;cosmos_logmessage, 'Bent case Lam max. ', /WARNING ;TS 2021-04-27
lam_max_p = (sin((tth_max)*0.5)/sin(twotheta/2.))*lambda_end ; max lam needed for bent case [bob 13/4/15]
lam_min_p = (sin((tth_min)*0.5)/sin(twotheta/2.))*lambda_ini ; min lam needed for bent case [bob 13/4/15]
endif else begin ;divergent case
;print,'divvy case'
lam_max_p = (sin((tth_max-twotheta/2.0)*0.5)/sin(twotheta/4.))*lambda_end ; max lam needed for divergent case [bob 13/4/15]
lam_min_p = (sin((tth_min-twotheta/2.0)*0.5)/sin(twotheta/4.))*lambda_ini ; min lam needed for divergent case [bob 13/4/15]
cosmos_logmessage, 'Divergent case Lam max. ', /WARNING
;cosmos_logmessage, 'Divergent case Lam max. ', /WARNING ;TS 2021-04-27
endelse
;print,lambda_ini,lambda_end,lam_min_p,lam_max_p
aa = floor((lam_min_p-lambda[0])/b_lam) ; wavelength bin of the lowest wavelength needed for lammin
bb = floor((lam_max_p-lambda[0])/b_lam) ; wavelength bin of highest wavelength needed for lammax
......@@ -3013,10 +2884,10 @@ pro cosmos_anal_run, result, groupedresult, direct, reflect, instr, theta, runno
gama = float(peakref-i)*b_tth ;2th shift from specular
if tra gt 250 then begin ;;TS 2020-02-15 NO analyzer option for coherent binning. Standard
gama = float(peakref-i)*b_tth ;2th shift from specular
cosmos_logmessage, ' NO ANALYZER b = !!!!!'
;cosmos_logmessage, ' NO ANALYZER b = !!!!!'
endif else begin ;TS 2020-02-15 with analyzer option for coherent binning. Need to mirror the reflection.
gama = -1.0*float(peakref-i)*b_tth ;2th shift from specular
cosmos_logmessage, ' WITH ANALYZER b = !!!!!'
;cosmos_logmessage, ' WITH ANALYZER b = !!!!!'
endelse
endelse
;TS 2015-06-01 ((detres/(2*parref.sdetd)) changed to ((detres/(parref.sdetd))
......@@ -3025,15 +2896,9 @@ pro cosmos_anal_run, result, groupedresult, direct, reflect, instr, theta, runno
if ((om_fwhm gt 0) AND ((detres/(parref.sdetd)) gt s2_fwhm)) then begin
lam_bot = sin((twotheta_arr[j])*0.5)*(lam-0.5*b_lam)/(sin((twotheta_arr[j]+gama+b_tth/2.)*0.5)) ;bot lam ;specular wavelength at same qz ;PG on 16/02/15: only valid for bent sample and collimated incoming beam
lam_top = sin((twotheta_arr[j])*0.5)*(lam+0.5*b_lam)/(sin((twotheta_arr[j]+gama-b_tth/2.)*0.5)) ;top lam ;specular wavelength at same qz ;PG on 16/02/15: only valid for bent sample and collimated incoming beam
;if j eq 300 then cosmos_logmessage, 'Bent sample coherent lam_bot + lam_top calculation ', /WARNING
;if j eq 400 then print, 'bent'
endif else begin
lam_bot = sin((twotheta_arr[j]+gama)*0.5)*(lam-0.5*b_lam)/(sin((twotheta_arr[j]+3.*gama+b_tth)*0.5)) ;For divergent incoming beam; PG on 28/04/15: b_tth should not be halfed here (see paper)
lam_top = sin((twotheta_arr[j]+gama)*0.5)*(lam+0.5*b_lam)/(sin((twotheta_arr[j]+3.*gama-b_tth)*0.5))
;if j eq 300 then cosmos_logmessage, 'Divergent sample coherent lam_bot + lam_top calculation ', /WARNING
;cosmos_logmessage, 'Divergent beam coherent lambda calculation ', /WARNING
;if j eq 400 then print, 'div'
;if j eq 500 then print,i,peakref,lam,lam_bot,lam_top,gamma*180./!pi,twotheta*180./!pi
;these are lines of constant q with the origin shifted by gamma
endelse
;calculate the limits in lambda on the chosen 2th line in pixels
......@@ -3085,64 +2950,38 @@ pro cosmos_anal_run, result, groupedresult, direct, reflect, instr, theta, runno
chop_res = (c_params.chopsep + (c_params.planckperkg * (parref.openangle+openofflam[lambda_range]) * parref.period / $
(360. * lambda[lambda_range] * 1.e-10))) / (2 * parref.tofd) ;0.5*max chopper lambda resolution
;chop_res = (c_params.chopsep + (c_params.planckperkg * parref.openangle * parref.period / $
; (360. * lambda[lambda_range] * 1.e-10))) / (2 * parref.tofd) ;0.5*max chopper lambda resolution
det_res = c_params.planckperkg * parref.channelwidth / $;calculate max detector lambda resolution
( lambda[lambda_range] * 1.e-10) / (2.0 * parref.tofd)
err_res = 0.98 * (3*chop_res^2 + det_res^2 + 3*chop_res*det_res)/(2*chop_res + det_res) ;FWHM of resulting trapezoid
;PG on 07/08/2017: Added wavelength smearing due to finite beam size at chopper 1 position:
tempratio=(parref.tofd-sdr)*1000/c_params.interslit
tempa=tempratio*abs(parref.s2w-parref.s3w)+parref.s2w
tempb=tempratio*(parref.s2w+parref.s3w)+parref.s2w
tempwidthfwhm=0.49*(tempb^3-tempa^3)/(tempb^2-tempa^2)
;print, 'new' , tempratio , tempa , tempb
;calculate beam width lambda resolution
width_res = tempwidthfwhm/1000.0 * parref.period / $
(2 * !pi * c_params.cr)* c_params.planckperkg / $
( lambda[lambda_range] * 1.e-10) / parref.tofd
err_res = sqrt(err_res^2+width_res^2)
;print, width_res
;err_res = sqrt(err_res^2)
;s3_fwhm=(0.68*parref.s3w)/c_params.interslit ; fwhm angles of second slit PG on 16/02/15: corrected for interslit distance
s3_fwhm=(0.68*parref.s3w)/(parref.slit3*1000.0+parref.sdetd*1000.0) ; fwhm angles of second slit PG on 28/04/15: corrected for interslit distance
;print,'s2 s3= ',s2_fwhm*180./!pi,s3_fwhm*180./!pi
;Resolution calculation (different depending on coherent/incoherent treatment)
if calc_method eq 'coherent' then begin
if om_fwhm gt 0 then begin
;TS 2015-06-01: corrected da for S3-det distance
;;TS 2015-06-01: distinguish between cases of divergence larger or smaller than bendiness of sample
if s2_fwhm ge 2*om_fwhm then err_ray = sqrt( (detres/(parref.slit3+parref.sdetd))^2 + s3_fwhm^2 + om_fwhm^2 )/angle_bragg ;PG on 16/2/15: I don;t think this is correct
if s2_fwhm ge 2*om_fwhm then cosmos_logmessage, 'Bent sample resolution with s2fwhm>2*om_fwhm ', /WARNING ;2018-06-06 TS
if s2_fwhm lt 2*om_fwhm then err_ray = sqrt( (detres/(2*(parref.slit3+parref.sdetd)))^2 + s3_fwhm^2 + s2_fwhm^2 )/angle_bragg ;PG on 16/2/15: I don;t think this is correct
if s2_fwhm lt 2*om_fwhm then cosmos_logmessage, 'Bent sample resolution with s2fwhm<2*om_fwhm ', /WARNING ;2018-06-06 TS
;cosmos_logmessage, 'detres ' + cosmos_str_make(detres), /WARNING ;2018-06-06 TS for debugging
;cosmos_logmessage, 'slit3 ' + cosmos_str_make(parref.slit3), /WARNING ;2018-06-06 TS for debugging
;cosmos_logmessage, 'sdet ' + cosmos_str_make(parref.sdetd), /WARNING ;2018-06-06 TS for debugging
;cosmos_logmessage, 's2w ' + cosmos_str_make(parref.s2w), /WARNING ;2018-06-06 TS for debugging
;cosmos_logmessage, 's3w ' + cosmos_str_make(parref.s3w), /WARNING ;2018-06-06 TS for debugging
;cosmos_logmessage, 'interslit ' + cosmos_str_make(c_params.interslit), /WARNING ;2018-06-06 TS for debugging
;cosmos_logmessage, 'om_fwhm ' + cosmos_str_make(om_fwhm), /WARNING ;2018-06-06 TS for debugging
;cosmos_logmessage, '2x om_fwhm ' + cosmos_str_make(2*om_fwhm), /WARNING ;2018-06-06 TS for debugging
;cosmos_logmessage, 's2_fwhm ' + cosmos_str_make(s2_fwhm), /WARNING ;2018-06-06 TS for debugging
;cosmos_logmessage, 's3_fwhm ' + cosmos_str_make(s3_fwhm), /WARNING ;2018-06-06 TS for debugging
;tester=detres/(2*(parref.slit3+parref.sdetd)) ;2018-06-06 TS for debugging
;cosmos_logmessage, 'detr ' + cosmos_str_make(tester), /WARNING ;2018-06-06 TS for debugging
;cosmos_logmessage, 'da ' + cosmos_str_make(da*180/!pi), /WARNING ;2018-06-06 TS for debugging
;cosmos_logmessage, 'angle_bragg ' + cosmos_str_make(angle_bragg), /WARNING ;2018-06-06 TS for debugging
;cosmos_logmessage, 'n err_ray ' + cosmos_str_make(err_ray), /WARNING ;2018-06-06 TS for debugging
;err_ray = sqrt((detres/(2*parref.sdetd))^2 + s3_fwhm^2+s2_fwhm^2)/angle_bragg ;This is the case for a bent sample
;TS 2015-06-01: corrected da for S3-det distance
;;TS 2015-06-01: distinguish between cases of divergence larger or smaller than bendiness of sample
if s2_fwhm ge 2*om_fwhm then err_ray = sqrt( (detres/(parref.slit3+parref.sdetd))^2 + s3_fwhm^2 + om_fwhm^2 )/angle_bragg ;PG on 16/2/15: I don;t think this is correct
;if s2_fwhm ge 2*om_fwhm then cosmos_logmessage, 'Bent sample resolution with s2fwhm>2*om_fwhm ', /WARNING ;TS 2021-04-27
if s2_fwhm lt 2*om_fwhm then err_ray = sqrt( (detres/(2*(parref.slit3+parref.sdetd)))^2 + s3_fwhm^2 + s2_fwhm^2 )/angle_bragg ;PG on 16/2/15: I don;t think this is correct
endif else begin
if s2_fwhm gt (detres/parref.sdetd) then begin
err_ray = sqrt((detres/(parref.slit3+parref.sdetd))^2 + s3_fwhm^2)/angle_bragg ;This is the case for a flat sample and a divergent incoming beam
cosmos_logmessage, 'Divergent beam resolution ', /WARNING
;cosmos_logmessage, 'Divergent beam resolution ', /WARNING
endif else begin
cosmos_logmessage, 'Constant q foreground COSMOS for angle ' + cosmos_str_make(runno) + ' is not appropriate!', /WARNING ;This is the case were the detector resolution is worse than the incoming divergence for a flat sample
err_ray = sqrt(da^2 + (detres/(parref.slit3+parref.sdetd))^2)/angle_bragg
......@@ -3161,9 +3000,7 @@ pro cosmos_anal_run, result, groupedresult, direct, reflect, instr, theta, runno
cosmos_logmessage, ' '
endif
; see paper - PG on 16/2/2015 decided to not use the sample bendiness in the incoherent case
;err_ray = sqrt(da^2+om_fwhm^2)/angle_bragg
err_ray = sqrt(da^2)/angle_bragg
;print, om_fwhm
if om_fwhm gt 0 then begin
if ((abs(parref.s2w-pardir.s2w) lt 0.04) and (abs(parref.s3w-pardir.s3w) lt 0.04)) then cosmos_logmessage, 'The sample is bent or there is a lot of diffuse scattering for angle ' + cosmos_str_make(runno) + '!' ;PG on 28/04/15
endif
......@@ -3171,7 +3008,6 @@ pro cosmos_anal_run, result, groupedresult, direct, reflect, instr, theta, runno
err_ray_temp = 0.68*sqrt((((total(foregdref)+1)*c_params.pixelwidth)^2 + (parref.s3w/1000)^2)/parref.sdetd^2)/angle_bragg ;PG on 16/2/2015
; PG: Choose the smaller of the two angular resolutions (either defined by incoming divergence or fgd width)
err_ray=min([err_ray,err_ray_temp])
cosmos_logmessage, 'Incoherent resolution ', /WARNING
endelse ;end of resolution calculation (coherent/incoherent)
totf = fgdref[1]-fgdref[0]
......@@ -3192,9 +3028,7 @@ pro cosmos_anal_run, result, groupedresult, direct, reflect, instr, theta, runno
cosmos_logmessage, ' constants.interslit = ' + cosmos_str_make(c_params.interslit)
cosmos_logmessage, ' sample waviness (fwhm) = ' + cosmos_str_make(om_fwhm_real*180./!pi) + 'deg, incoming angular divergence (fwhm) = ' + cosmos_str_make(da*180/!pi) + 'deg'
cosmos_logmessage, ' angle_bragg = ' + cosmos_str_make(angle_bragg) + ' rad = ' + cosmos_str_make(angle_bragg*180./3.1416) + 'deg'
;cosmos_logmessage, ' twotheta of chosen pixel line = ' + cosmos_str_make(twotheta*180./3.1416) + 'deg'
cosmos_logmessage, ' err_ray = ' + cosmos_str_make(err_ray)
cosmos_logmessage, ' TRA = ' + cosmos_str_make(tra) ;;TS 2020-02-15
cosmos_logmessage, ' '
result = fltarr(5, result_size) ;[Q, dQ, I, dI, lambda]
......@@ -3257,7 +3091,6 @@ pro cosmos_anal_run, result, groupedresult, direct, reflect, instr, theta, runno
endif else begin
cosmos_logmessage, 'Reference Q values cannot be used. Different number of points!', /warning
cosmos_logmessage, 'Try to slightly adjust the Min/Max Lambda', /warning ;TS 2016-10-21 Had a problem that the wavelength is sometimes one point longer...
;print, 'lambda length problem'
endelse
endif
......@@ -3599,9 +3432,6 @@ endif else begin
;reflection from downwards
newtheta = (theta0 - deltath) * !radeg ;effective incident angle in deg
endelse
print,'newtheta',(theta0* !radeg),newtheta
print,'poffoff',poffoff
print,'openofflam',openofflam
end
......@@ -3656,11 +3486,7 @@ pro cosmos_gravnew,theta0,lambda,newtheta,delta,dz,poffoff,openofflam,xchopper,x
newtheta = -1.0 * newtheta * !radeg
theta0 = -1.0 * theta0
endif else newtheta = newtheta * !radeg
; print,'newtheta',(theta0* !radeg),newtheta
; print,'poffoff',poffoff
; print,'openofflam',openofflam
end
......@@ -3800,7 +3626,7 @@ pro cosmos_raw_read, files, monitors, counts, params, path, override, reflect=re
ok = file_test(testfile)
if ok then begin
;cosmos_logmessage, 'Parameters coming from ASCII file' ;TS 2021-04-27 Removed the log message
for f = 0, n_elements(files) - 1 do begin
if strcompress(files[f]) gt ' ' then begin
......@@ -3931,7 +3757,7 @@ pro cosmos_raw_read, files, monitors, counts, params, path, override, reflect=re
cosmos_logmessage, 'No valid value for TOF delay in numor ' + file_name, /WARNING
cosmos_logmessage, 'Using default value = ' + cosmos_str_make(c_params.edelay), /WARNING
endelse
endif else if (strcmp(strtrim(c_params.inst,2),'figaro',/FOLD_CASE)) then begin
; poff
......@@ -4031,21 +3857,7 @@ pro cosmos_raw_read, files, monitors, counts, params, path, override, reflect=re
;MAG: Setting arrays sizes from numor
if (strcmp(strtrim(c_params.inst,2),'d17',/FOLD_CASE)) then begin
; History of detector centers
case c_params.data_format of
1: c_params.pcen = 135.79
2: c_params.pcen = 132.5
3: c_params.pcen = 132.5
4: c_params.pcen = 132.5
5: c_params.pcen = 132.5
else: begin
cosmos_logmessage, 'c_params.data_format = ' + string(c_params.data_format) + ' not handled!', /WARNING
cosmos_logmessage, 'Using c_params.pcen = 132.5', /WARNING
c_params.pcen = 132.5
end
endcase
;pixels_x
tmp = long(par1[98]-par1[97]+1)
if (tmp gt 1) then begin
......@@ -4127,6 +3939,7 @@ pro cosmos_raw_read, files, monitors, counts, params, path, override, reflect=re
; 100 : y2
; 101 : nx, detector grouping factor in x (ie nx = 1 => base detector has 286 pixels in x, nx = 2 => base detector has 144 pixels in x)
; 102 : ny, detector grouping factor in y (ny = 1 => 276 pixels)
; 105 : pcen, center of detector
; params2[ ] :
; 2 : san, sample angle
; 15 : sample to detector distance (millimetres)
......@@ -4144,9 +3957,36 @@ pro cosmos_raw_read, files, monitors, counts, params, path, override, reflect=re
if (strcmp(strtrim(c_params.inst,2),'d17',/FOLD_CASE)) then begin
;TS 2021-02-25
; History of detector centers
; Pcen
case c_params.data_format of
1: c_params.pcen = 135.79
2: c_params.pcen = 132.5
3: c_params.pcen = 132.5
4: c_params.pcen = 132.5
5: c_params.pcen = 132.5 ;TS 2021-02-25 corrected the value
else: begin
cosmos_logmessage, 'c_params.data_format = ' + string(c_params.data_format) + ' not handled!', /WARNING
cosmos_logmessage, 'Using c_params.pcen = 132.5', /WARNING ;TS 2021-02-25 corrected the value
c_params.pcen = 132.5 ;TS 2021-02-25 corrected the value
end
endcase
;Pcen which is now written in datafile ;TS 2021-02-25
tmp = par1[105]
if (tmp gt 0.0) then begin
c_params.pcen = tmp
endif else begin
cosmos_logmessage, 'No valid value for Pcen in numor ' + file_name, /WARNING
cosmos_logmessage, 'Using default value = ' + cosmos_str_make(c_params.pcen), /WARNING
endelse
; TS 2021-02-25 added reading of chopper parameter from ascii file depending on the ascii version
;chopsep + chopwin + (cr) chopradius
;until cycle 183, chopper gap was correctly written in par2[14],
;but this is not the case from cycle 191 onwards, so relying on
;hard-coded values only
;hard-coded values
;From cycle 211 the values are written in the ascii file par2[14] ;TS 2021-02-25
if c_params.data_format le 2 then begin
c_params.chopsep = c_params.chopsep_d17[0]
c_params.chopwin = c_params.chopwin_d17[0]
......@@ -4155,12 +3995,47 @@ pro cosmos_raw_read, files, monitors, counts, params, path, override, reflect=re
c_params.chopsep = c_params.chopsep_d17[1]
c_params.chopwin = c_params.chopwin_d17[1]
c_params.cr = c_params.cr_d17[1]
endif else begin
endif else if c_params.data_format eq 4 then begin
c_params.chopsep = c_params.chopsep_d17[2]
c_params.chopwin = c_params.chopwin_d17[2]
c_params.cr = c_params.cr_d17[2]
endif else begin ; TS 2021-02-25 added reading of chopper parameter from ascii file
;read in the chopper separation. New from cycle 211, Ascii version 5
tmp = float(par2[14])
if (tmp gt 0.0) then begin
if (tmp le 10.0) then begin
c_params.chopsep = tmp ;
endif else begin
c_params.chopsep = tmp/1000 ;
endelse
endif else begin
cosmos_logmessage, 'No valid value for chopsep in numor ' + file_name, /WARNING
cosmos_logmessage, 'Using default value = ' + cosmos_str_make(c_params.chopsep_d17[2]), /WARNING
c_params.chopsep = c_params.chopsep_d17[2]