Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

V_DetectorCorrections.ipf

Last change on this file was 1242, checked in by srkline, 2 years ago
updating the IgorVersion? pragma to v7.0 for all files to be consistent.
File size: 68.4 KB

Line
1	#pragma rtGlobals=3 // Use modern global access method and strict wave access.
2	#pragma version=1.0
3	#pragma IgorVersion = 7.00
4
5
6
7	//
8	// functions to apply corrections to the detector panels
9	//
10	// these are meant to be called by the procedures that convert "raw" data to
11	// "adjusted" or corrected data sets
12	//
13
14	////////////////
15	// Constants for detector efficiency and shadowing
16	//
17	// V_TubeEfficiencyShadowCorr()
18	//
19	// JAN 2020
20	///////////////
21	Constant kTube_ri = 0.372 // inner radius of tube [cm]
22	Constant kTube_cc = 0.84 // center to center spacing [cm]
23	Constant kTube_ss = 0.025 // stainless steel shell thickness [cm]
24
25	Constant kSig_2b_He = 0.146 // abs xs for 2 bar He(3) [cm-1 A-1] (multiply this by wavelength)
26	Constant kSig_8b_He = 0.593 // abs xs for 8 bar He(3) [cm-1 A-1] (multiply this by wavelength)
27	Constant kSig_Al = 0.00967 // abs xs for Al [cm-1 A-1] (multiply this by wavelength)
28	Constant kSig_ss = 0.146 // abs xs for 304 SS [cm-1 A-1] (multiply this by wavelength)
29
30
31
32
33	//
34	// detector dead time
35	//
36	// input is the data array (N tubes x M pixels)
37	// input of N x 1 array of dead time values
38	//
39	// output is the corrected counts in data, overwriting the input data
40	//
41	// Note that the equation in Roe (eqn 2.15, p. 63) looks different, but it is really the
42	// same old equation, just written in a more complex form.
43	//
44	// (DONE)
45	// x- verify the direction of the tubes and indexing
46	// x- decide on the appropriate functional form for the tubes
47	// x- need count time as input
48	// x- be sure I'm working in the right data folder (all waves are passed in)
49	// x- clean up when done
50	// x- calculate + return the error contribution?
51	// x- verify the error propagation
52	Function V_DeadTimeCorrectionTubes(dataW,data_errW,dtW,ctTime)
53	Wave dataW,data_errW,dtW
54	Variable ctTime
55
56	// do I count on the orientation as an input, or do I just figure it out on my own?
57	String orientation
58	Variable dimX,dimY
59	dimX = DimSize(dataW,0)
60	dimY = DimSize(dataw,1)
61	if(dimX > dimY)
62	orientation = "horizontal"
63	else
64	orientation = "vertical"
65	endif
66
67	// sum the counts in each tube and divide by time for total cr per tube
68	Variable npt
69
70	if(cmpstr(orientation,"vertical")==0)
71	// this is data dimensioned as (Ntubes,Npix)
72
73	MatrixOp/O sumTubes = sumRows(dataW) // n x 1 result
74	sumTubes /= ctTime //now count rate per tube
75
76	dataW[][] = dataW[p][q]/(1-sumTubes[p]*dtW[p]) //correct the data
77	data_errW[][] = data_errW[p][q]/(1-sumTubes[p]*dtW[p]) // propagate the error wave
78
79	elseif(cmpstr(orientation,"horizontal")==0)
80	// this is data (horizontal) dimensioned as (Npix,Ntubes)
81
82	MatrixOp/O sumTubes = sumCols(dataW) // 1 x m result
83	sumTubes /= ctTime
84
85	dataW[][] = dataW[p][q]/(1-sumTubes[q]*dtW[q])
86	data_errW[][] = data_errW[p][q]/(1-sumTubes[q]*dtW[q])
87
88	else
89	DoAlert 0,"Orientation not correctly passed in DeadTimeCorrectionTubes(). No correction done."
90	endif
91
92	return(0)
93	end
94
95	// test function
96	Function V_testDTCor()
97
98	String detStr = ""
99	String fname = "RAW"
100	Variable ctTime
101
102	detStr = "FR"
103	Wave w = V_getDetectorDataW(fname,detStr)
104	Wave w_err = V_getDetectorDataErrW(fname,detStr)
105	Wave w_dt = V_getDetector_deadtime(fname,detStr)
106
107	ctTime = V_getCount_time(fname)
108
109	// ctTime = 10
110	V_DeadTimeCorrectionTubes(w,w_err,w_dt,ctTime)
111
112	End
113
114
115	//
116	// Non-linear data correction
117	//
118	// DOES NOT modify the data, only calculates the spatial relationship
119	//
120	// input is the data array (N tubes x M pixels)
121	// input of N x M array of quadratic coefficients
122	//
123	// output is wave of corrected real space distance corresponding to each pixel of the data
124	// its distance from the nominal beam center of (0,0)
125	//
126	//
127	// (DONE)
128	// x- UNITS!!!! currently this is mm, which certainly doesn't match anything else!!!
129	//
130	// x- verify the direction of the tubes and indexing
131	// x- be sure I'm working in the right data folder (it is passed in, and the full path is used)
132	// x- clean up when done
133	// x- calculate + return the error contribution? (there is none for this operation)
134	// x- do I want this to return a wave? (no, default names are generated)
135	// x- do I need to write a separate function that returns the distance wave for later calculations?
136	// x- do I want to make the distance array 3D to keep the x and y dims together? Calculate them all right now?
137	// x- what else do I need to pass to the function? (fname=folder? detStr?)
138	// y- (yes,see below) need a separate block or function to handle "B" detector which will be ? different
139	//
140	//
141	Function V_NonLinearCorrection(fname,dataW,coefW,tube_width,detStr,destPath)
142	String fname //can also be a folder such as "RAW"
143	Wave dataW,coefW
144	Variable tube_width
145	String detStr,destPath
146
147
148	// do I count on the orientation as an input, or do I just figure it out on my own?
149	String orientation
150	Variable dimX,dimY
151	dimX = DimSize(dataW,0)
152	dimY = DimSize(dataW,1)
153	if(dimX > dimY)
154	orientation = "horizontal"
155	else
156	orientation = "vertical"
157	endif
158
159	// make a wave of the same dimensions, in the same data folder for the distance
160	// ?? or a 3D wave?
161	Make/O/D/N=(dimX,dimY) $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
162	Make/O/D/N=(dimX,dimY) $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
163	Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
164	Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
165
166	// then per tube, do the quadratic calculation to get the real space distance along the tube
167	// the distance perpendicular to the tube is n*(8.4mm) per tube index
168
169	// DONE
170	// -- GAP was hard-wired, but in 2018 proper values for all 4 gaps were measured
171	// and added to the file header for each detector panel. there is now a read from the
172	// header to get the gap value
173	Variable offset,gap
174
175	gap = V_getDet_panel_gap(fname,detStr)
176
177	// DONE:
178	// -- in case of error, V_getDet_panel_gap() will return -999999
179	// -- it should only apply to data pre-2018 when the field did not exist in the file
180	// -- any VSANS data from 2018+ should read gap from the file and bypass the if()
181
182	if(gap < -100) //-999999 returned if field is missing from file
183
184	if(cmpstr(detStr,"FL") == 0 \|\| cmpstr(detStr,"FR") == 0)
185	gap = 3.5 //mm (measured, JB 1/4/18)
186	endif
187	if(cmpstr(detStr,"FT") == 0 \|\| cmpstr(detStr,"FB") == 0)
188	gap = 3.3 //mm (measured, JB 2/1/18)
189	endif
190	if(cmpstr(detStr,"ML") == 0 \|\| cmpstr(detStr,"MR") == 0)
191	gap = 5.9 //mm (measured, JB 1/4/18)
192	endif
193	if(cmpstr(detStr,"MT") == 0 \|\| cmpstr(detStr,"MB") == 0)
194	gap = 18.3 //mm (measured, JB 2/1/18)
195	endif
196
197	endif
198
199
200	if(cmpstr(orientation,"vertical")==0)
201	// this is data dimensioned as (Ntubes,Npix)
202
203	// adjust the x postion based on the beam center being nominally (0,0) in units of cm, not pixels
204	if(cmpstr(fname,"VCALC")== 0 )
205	offset = VCALC_getPanelTranslation(detStr)
206	offset *= 10 // convert to units of mm
207	// if(cmpstr("L",detStr[1]) == 0)
208	// offset *= -1 //negative value for L
209	// endif
210	else
211	//normal case
212	offset = V_getDet_LateralOffset(fname,detStr)
213	offset *= 10 //convert cm to mm
214	endif
215
216	// calculation is in mm, not cm
217	// offset will be a negative value for the L panel, and positive for the R panel
218	if(kBCTR_CM)
219	if(cmpstr("L",detStr[1]) == 0)
220	// data_realDistX[][] = offset - (dimX - p)*tube_width // TODO should this be dimX-1-p = 47-p?
221	data_realDistX[][] = offset - (dimX - p - 1/2)*tube_width - gap/2 // TODO should this be dimX-1-p = 47-p?
222	else
223	// right
224	// data_realDistX[][] = tube_width*(p+1) + offset + gap //add to the Right det,
225	data_realDistX[][] = tube_width*(p+1/2) + offset + gap/2 //add to the Right det
226	endif
227	else
228	data_realDistX[][] = tube_width*(p)
229	endif
230	data_realDistY[][] = coefW[0][p] + coefW[1][p]q + coefW[2][p]q*q
231
232
233	elseif(cmpstr(orientation,"horizontal")==0)
234	// this is data (horizontal) dimensioned as (Npix,Ntubes)
235	data_realDistY[][] = tube_width*q
236
237	if(cmpstr(fname,"VCALC")== 0 )
238	offset = VCALC_getPanelTranslation(detStr)
239	offset *= 10 // convert to units of mm
240	// if(cmpstr("B",detStr[1]) == 0)
241	// offset *= -1 // negative value for Bottom det
242	// endif
243	else
244	//normal case
245	offset = V_getDet_VerticalOffset(fname,detStr)
246	offset *= 10 //convert cm to mm
247	endif
248
249	if(kBCTR_CM)
250	if(cmpstr("T",detStr[1]) == 0)
251	// data_realDistY[][] = tube_width*(q+1) + offset + gap
252	data_realDistY[][] = tube_width*(q+1/2) + offset + gap/2
253	else
254	// bottom
255	// data_realDistY[][] = offset - (dimY - q)*tube_width // TODO should this be dimY-1-q = 47-q?
256	data_realDistY[][] = offset - (dimY - q - 1/2)*tube_width - gap/2 // TODO should this be dimY-1-q = 47-q?
257	endif
258	else
259	data_realDistY[][] = tube_width*(q)
260	endif
261	data_realDistX[][] = coefW[0][q] + coefW[1][q]p + coefW[2][q]p*p
262
263	else
264	DoAlert 0,"Orientation not correctly passed in NonLinearCorrection(). No correction done."
265	return(0)
266	endif
267
268	return(0)
269	end
270
271
272
273
274	// TODO:
275	// -- the cal_x and y coefficients are totally fake
276	// -- the wave assignment may not be correct.. so beware
277	//
278	//
279	Function V_NonLinearCorrection_B(folder,dataW,cal_x,cal_y,detStr,destPath)
280	String folder
281	Wave dataW,cal_x,cal_y
282	String detStr,destPath
283
284	if(cmpstr(detStr,"B") != 0)
285	return(0)
286	endif
287
288	Print "*Cal_X and Cal_Y for Back are using file values *"
289
290	// cal_x[0] = VCALC_getPixSizeX(detStr)*10 // pixel size in mm VCALC_getPixSizeX(detStr) is [cm]
291	// cal_x[1] = 1
292	// cal_x[2] = 10000
293	// cal_y[0] = VCALC_getPixSizeY(detStr)*10 // pixel size in mm VCALC_getPixSizeX(detStr) is [cm]
294	// cal_y[1] = 1
295	// cal_y[2] = 10000
296
297
298	// do I count on the orientation as an input, or do I just figure it out on my own?
299	Variable dimX,dimY
300
301	// Wave dataW = V_getDetectorDataW(folder,detStr)
302
303	dimX = DimSize(dataW,0)
304	dimY = DimSize(dataW,1)
305
306	// make a wave of the same dimensions, in the same data folder for the distance
307	// ?? or a 3D wave?
308	Make/O/D/N=(dimX,dimY) $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
309	Make/O/D/N=(dimX,dimY) $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
310	Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
311	Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
312
313	// Wave cal_x = V_getDet_cal_x(folder,detStr)
314	// Wave cal_y = V_getDet_cal_y(folder,detStr)
315
316	data_realDistX[][] = cal_x[0]p10 // cal_x and cal_y are in [cm], need mm
317	data_realDistY[][] = cal_y[0]q10
318
319	return(0)
320	end
321
322
323	//
324	//
325	// TODO
326	// -- VERIFY the calculations
327	// -- verify where this needs to be done (if the beam center is changed)
328	// -- then the q-calculation needs to be re-done
329	// -- the position along the tube length is referenced to tube[0], for no particular reason
330	// It may be better to take an average? but [0] is an ASSUMPTION
331	// -- distance along tube is simple interpolation, or do I use the coefficients to
332	// calculate the actual value
333	//
334	// -- distance in the lateral direction is based on tube width, which is a fixed parameter
335	//
336	//
337	Function V_ConvertBeamCtrPix_to_mm(folder,detStr,destPath)
338	String folder,detStr,destPath
339
340	Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
341	Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
342
343	String orientation
344	Variable dimX,dimY,xCtr,yCtr
345	dimX = DimSize(data_realDistX,0)
346	dimY = DimSize(data_realDistX,1)
347	if(dimX > dimY)
348	orientation = "horizontal"
349	else
350	orientation = "vertical"
351	endif
352
353	xCtr = V_getDet_beam_center_x(folder,detStr)
354	yCtr = V_getDet_beam_center_y(folder,detStr)
355
356	Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
357	Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
358	WAVE x_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
359	WAVE y_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
360
361	Variable tube_width = V_getDet_tubeWidth(folder,detStr) //this is also in mm
362
363	//
364	strswitch(detStr) // string switch
365	case "FL":
366	case "ML":
367	// for Left/Right
368	// for left
369	x_mm[0] = data_realDistX[dimX-1][0] + (xCtr-dimX-1)*tube_width
370	y_mm[0] = data_realDistY[0][yCtr]
371
372	break
373	case "FR":
374	case "MR":
375	// for Left/Right
376	// for right
377	x_mm[0] = data_realDistX[0][0] + xCtr*tube_width
378	y_mm[0] = data_realDistY[0][yCtr]
379
380	break
381	case "FT":
382	case "MT":
383	// for Top
384	x_mm[0] = data_realDistX[xCtr][0]
385	y_mm[0] = data_realDistY[0][0] + yCtr*tube_width
386
387	break
388	case "FB":
389	case "MB":
390	// for Bottom
391	x_mm[0] = data_realDistX[xCtr][0]
392	y_mm[0] = data_realDistY[0][dimY-1] + (yCtr-dimY-1)*tube_width
393
394	break
395	default: // optional default expression executed
396	Print "No case matched in V_Convert_FittedPix_2_cm"
397	return(1)
398	endswitch
399
400	return(0)
401	end
402
403	//
404	//
405	// (DONE)
406	// x- VERIFY the calculations
407	// x- verify where this needs to be done (if the beam center is changed)
408	// x- then the q-calculation needs to be re-done
409	// x- the position along the tube length is referenced to tube[0], for no particular reason
410	// It may be better to take an average? but [0] is an ASSUMPTION
411	// x- distance along tube is simple interpolation
412	//
413	// x- distance in the lateral direction is based on tube width, which is a fixed parameter
414	//
415	// the value in pixels is written to the local data folder, NOT to disk (it is recalculated as needed)
416	//
417	Function V_ConvertBeamCtr_to_pix(folder,detStr,destPath)
418	String folder,detStr,destPath
419
420	Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
421	Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
422
423	String orientation
424	Variable dimX,dimY,xCtr,yCtr
425	dimX = DimSize(data_realDistX,0)
426	dimY = DimSize(data_realDistX,1)
427	if(dimX > dimY)
428	orientation = "horizontal"
429	else
430	orientation = "vertical"
431	endif
432
433	xCtr = V_getDet_beam_center_x(folder,detStr) //these are in cm
434	yCtr = V_getDet_beam_center_y(folder,detStr)
435
436	Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_pix")
437	Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_pix")
438	WAVE x_pix = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_pix")
439	WAVE y_pix = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_pix")
440
441	Variable tube_width = V_getDet_tubeWidth(folder,detStr)
442
443	variable edge,delta
444	Variable gap
445
446	// the gap is split equally between the panel pairs
447	// DONE -- replace hard-wired values with V_getDet_panel_gap(fname,detStr) once it is added to the file
448
449	gap = V_getDet_panel_gap(folder,detStr)
450
451	// DONE:
452	// -- check in case of error, value should be read from header
453	// -- it should only apply to data pre-2018 when the field did not exist in the file
454	// -- any VSANS data from 2018+ should read gap from the file.
455
456	if(gap < -100) //-999999 returned if field is missing from file
457	if(cmpstr(detStr,"FL") == 0 \|\| cmpstr(detStr,"FR") == 0)
458	gap = 3.5 //mm (measured, JB 1/4/18)
459	endif
460	if(cmpstr(detStr,"FT") == 0 \|\| cmpstr(detStr,"FB") == 0)
461	gap = 3.3 //mm (measured, JB 2/1/18)
462	endif
463	if(cmpstr(detStr,"ML") == 0 \|\| cmpstr(detStr,"MR") == 0)
464	gap = 5.9 //mm (measured, JB 1/4/18)
465	endif
466	if(cmpstr(detStr,"MT") == 0 \|\| cmpstr(detStr,"MB") == 0)
467	gap = 18.3 //mm (measured, JB 2/1/18)
468	endif
469	endif
470
471	//
472	if(cmpstr(orientation,"vertical")==0)
473	// this is data dimensioned as (Ntubes,Npix)
474
475	if(kBCTR_CM)
476	if(cmpstr("L",detStr[1]) == 0)
477	Make/O/D/N=(dimX) tmpTube
478	tmpTube = data_RealDistX[p][0]
479	FindLevel/P/Q tmpTube xCtr*10
480	if(V_Flag)
481	edge = data_realDistX[47][0] //tube 47
482	delta = abs(xCtr*10 - edge)
483	x_pix[0] = dimX-1 + delta/tube_width
484	else
485	// beam center is on the panel, report the pixel value
486	x_pix[0] = V_LevelX
487	endif
488
489	else
490	// R panel
491	Make/O/D/N=(dimX) tmpTube
492	tmpTube = data_RealDistX[p][0]
493	FindLevel/P/Q tmpTube xCtr*10
494	if(V_Flag)
495	//level not found
496	edge = data_realDistX[0][0]
497	delta = abs(xCtr*10 - edge + gap) // how far past the edge of the panel
498	x_pix[0] = -delta/tube_width //since the left edge of the R panel is pixel 0
499	else
500	// beam center is on the panel, report the pixel value
501	x_pix[0] = V_LevelX
502	endif
503
504	endif
505
506	endif
507
508	// the y-center will be on the panel in this direction
509	Make/O/D/N=(dimY) tmpTube
510	tmpTube = data_RealDistY[0][p]
511	FindLevel /P/Q tmpTube, yCtr*10
512
513	y_pix[0] = V_levelX
514	KillWaves/Z tmpTube
515	// Print x_pix[0],y_pix[0]
516
517	else
518	// this is data (horizontal) dimensioned as (Npix,Ntubes)
519
520	if(kBCTR_CM)
521	if(cmpstr("T",detStr[1]) == 0)
522	Make/O/D/N=(dimY) tmpTube
523	tmpTube = data_RealDistY[p][0]
524	FindLevel/P/Q tmpTube yCtr*10
525	if(V_Flag)
526	edge = data_realDistY[0][0] //tube 0
527	delta = abs(yCtr*10 - edge + gap)
528	y_pix[0] = -delta/tube_width //since the bottom edge of the T panel is pixel 0
529	else
530	y_pix[0] = V_LevelX
531	endif
532
533	else
534	// FM(B) panel
535	Make/O/D/N=(dimY) tmpTube
536	tmpTube = data_RealDistY[p][0]
537	FindLevel/P/Q tmpTube yCtr*10
538	if(V_Flag)
539	edge = data_realDistY[0][47] //y tube 47
540	delta = abs(yCtr*10 - edge)
541	y_pix[0] = dimY-1 + delta/tube_width //since the top edge of the B panels is pixel 47
542	else
543	y_pix[0] = V_LevelX
544	endif
545
546	endif
547	endif
548
549	// the x-center will be on the panel in this direction
550	Make/O/D/N=(dimX) tmpTube
551	tmpTube = data_RealDistX[p][0]
552	FindLevel /P/Q tmpTube, xCtr*10
553
554	x_pix[0] = V_levelX
555	KillWaves/Z tmpTube
556
557	endif
558
559	return(0)
560	end
561
562	// converts from [cm] beam center to pixels
563	//
564	// the value in pixels is written to the local data folder, NOT to disk (it is recalculated as needed)
565	//
566	Function V_ConvertBeamCtr_to_pixB(folder,detStr,destPath)
567	String folder,detStr,destPath
568
569	Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
570	Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
571
572	Variable dimX,dimY,xCtr,yCtr
573	dimX = DimSize(data_realDistX,0)
574	dimY = DimSize(data_realDistX,1)
575
576	xCtr = V_getDet_beam_center_x(folder,detStr) //these are in cm, *10 to get mm
577	yCtr = V_getDet_beam_center_y(folder,detStr)
578
579	Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_pix")
580	Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_pix")
581	WAVE x_pix = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_pix")
582	WAVE y_pix = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_pix")
583
584
585	// simple wave lookup
586	// can't use x_pix[0] = data_RealDistX(xCtr)[0] since the data has no x-scale and (xCtr) is interpreted
587	// as a point value
588
589	//
590	//xCtr, yCtr are in cm, *10 to get mm to compare to distance array
591
592	Make/O/D/N=(dimX) tmpTube
593	tmpTube = data_RealDistX[p][0]
594	FindLevel /P/Q tmpTube, xCtr*10
595
596	x_pix[0] = V_levelX
597	KillWaves/Z tmpTube
598
599
600	Make/O/D/N=(dimY) tmpTube
601	tmpTube = data_RealDistY[0][p]
602	FindLevel /P/Q tmpTube, yCtr*10
603
604	y_pix[0] = V_levelX
605	KillWaves/Z tmpTube
606
607	print "pixel ctr B = ",x_pix[0],y_pix[0]
608
609	return(0)
610	end
611
612	//
613	//
614	// TODO
615	// -- VERIFY the calculations
616	// -- verify where this needs to be done (if the beam center is changed)
617	// -- then the q-calculation needs to be re-done
618	//
619	// -- not much is known about the "B" detector, so this
620	// all hinges on the non-linear corrections being done correctly for that detector
621	//
622	// Variable detCtrX, detCtrY
623	// // get the pixel center of the detector (not the beam center)
624	// detCtrX = trunc( DimSize(dataW,0)/2 ) //
625	// detCtrY = trunc( DimSize(dataW,1)/2 )
626	//
627	//
628	Function V_ConvertBeamCtrPix_to_mmB(folder,detStr,destPath)
629	String folder,detStr,destPath
630
631
632	// DoAlert 0,"Error - Beam center is being interpreted as pixels, but needs to be in cm. V_ConvertBeamCtrPix_to_mmB()"
633
634	Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
635	Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
636
637	Variable xCtr,yCtr
638	xCtr = V_getDet_beam_center_x(folder,detStr)
639	yCtr = V_getDet_beam_center_y(folder,detStr)
640
641	Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
642	Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
643	WAVE x_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
644	WAVE y_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
645
646	x_mm[0] = data_realDistX[xCtr][0]
647	y_mm[0] = data_realDistY[0][yCtr]
648
649	return(0)
650	end
651
652
653
654
655
656
657	/////
658	//
659	// non-linear corrections to the tube pixels
660	// - returns the distance in mm (although this may change)
661	//
662	// c0,c1,c2,pix
663	// c0-c2 are the fit coefficients
664	// pix is the test pixel
665	//
666	// returns the distance in mm (relative to ctr pixel)
667	// ctr is the center pixel, as defined when fitting to quadratic was done
668	//
669	Function V_TubePixel_to_mm(c0,c1,c2,pix)
670	Variable c0,c1,c2,pix
671
672	Variable dist
673	dist = c0 + c1pix + c2pix*pix
674
675	return(dist)
676	End
677	//
678	////
679
680
681	//
682	// TESTING ONLY
683	Proc V_MakeFakeCalibrationWaves()
684	// make these in the RAW data folder, before converting to a work folder
685	// - then they will be "found" by get()
686	// -- only for the tube, not the Back det
687
688	// DoAlert 0, "re-do this and do a better job of filling the fake calibration data"
689
690	DoAlert 0, "Calibration waves are read in from the data file"
691
692	// V_fMakeFakeCalibrationWaves()
693	End
694
695
696
697	//
698	// TESTING ONLY
699	//
700	// orientation does not matter, there are 48 tubes in each bank
701	// so dimension (3,48) for everything.
702	//
703	// -- but the orientation does indicate TB vs LR, which has implications for
704	// the (fictional) dimension of the pixel along the tube axis, at least as far
705	// as for making the fake coefficients.
706	//
707	Function V_fMakeFakeCalibrationWaves()
708
709	Variable ii,pixSize
710	String detStr,fname="RAW",orientation
711
712	for(ii=0;ii<ItemsInList(ksDetectorListNoB);ii+=1)
713	detStr = StringFromList(ii, ksDetectorListNoB, ";")
714	// Wave w = V_getDetectorDataW(fname,detStr)
715	Make/O/D/N=(3,48) $("root:Packages:NIST:VSANS:RAW:entry:instrument:detector_"+detStr+":spatial_calibration")
716	Wave calib = $("root:Packages:NIST:VSANS:RAW:entry:instrument:detector_"+detStr+":spatial_calibration")
717	// !!!! this overwrites what is there
718
719	orientation = V_getDet_tubeOrientation(fname,detStr)
720	if(cmpstr(orientation,"vertical")==0)
721	// this is vertical tube data dimensioned as (Ntubes,Npix)
722	pixSize = 8.4 //V_getDet_y_pixel_size(fname,detStr)
723
724	elseif(cmpstr(orientation,"horizontal")==0)
725	// this is data (horizontal) dimensioned as (Npix,Ntubes)
726	pixSize = 4 //V_getDet_x_pixel_size(fname,detStr)
727
728	else
729	DoAlert 0,"Orientation not correctly passed in NonLinearCorrection(). No correction done."
730	endif
731
732	calib[0][] = -(128/2)pixSize //approx (n/2)pixSixe
733	calib[1][] = pixSize
734	calib[2][] = 2e-4
735
736	endfor
737
738	// now fake calibration for "B"
739	Wave cal_x = V_getDet_cal_x("RAW","B")
740	Wave cal_y = V_getDet_cal_y("RAW","B")
741
742	cal_x = .34 // mm, ignore the other 2 values
743	cal_y = .34 // mm
744	return(0)
745	End
746
747	//
748	// (DONE)
749	// x- MUST VERIFY the definition of SDD and how (if) setback is written to the data files
750	// x- currently I'm assuming that the SDD is the "nominal" value which is correct for the
751	// L/R panels, but is not correct for the T/B panels (must add in the setback)
752	//
753	//
754	//
755	// data_realDistX, Y must be previously generated from running NonLinearCorrection()
756	//
757	// call with:
758	// fname as the WORK folder, "RAW"
759	// detStr = detector String, "FL"
760	// destPath = path to destination WORK folder ("root:Packages:NIST:VSANS:"+folder)
761	//
762	Function V_Detector_CalcQVals(fname,detStr,destPath)
763	String fname,detStr,destPath
764
765	String orientation
766	Variable xCtr,yCtr,lambda,sdd
767
768	// get all of the geometry information
769	orientation = V_getDet_tubeOrientation(fname,detStr)
770
771
772	sdd = V_getDet_ActualDistance(fname,detStr) //sdd derived, including setback [cm]
773
774	// this is the ctr in pixels --xx-- (now it is in cm!)
775	// xCtr = V_getDet_beam_center_x(fname,detStr)
776	// yCtr = V_getDet_beam_center_y(fname,detStr)
777	// this is ctr in mm
778	xCtr = V_getDet_beam_center_x_mm(fname,detStr)
779	yCtr = V_getDet_beam_center_y_mm(fname,detStr)
780	lambda = V_getWavelength(fname)
781	Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
782	Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
783
784	// make the new waves
785	Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qTot_"+detStr)
786	Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qx_"+detStr)
787	Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qy_"+detStr)
788	Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qz_"+detStr)
789	Wave qTot = $(destPath + ":entry:instrument:detector_"+detStr+":qTot_"+detStr)
790	Wave qx = $(destPath + ":entry:instrument:detector_"+detStr+":qx_"+detStr)
791	Wave qy = $(destPath + ":entry:instrument:detector_"+detStr+":qy_"+detStr)
792	Wave qz = $(destPath + ":entry:instrument:detector_"+detStr+":qz_"+detStr)
793
794	// calculate all of the q-values
795	// sdd is passed in [cm]
796
797
798	// after adding in the 680x1656 back detector, load time was 7.8s, without multithreading
799	// with multithreading, 1.9s
800	// qTot = V_CalcQval(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
801	// qx = V_CalcQX(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
802	// qy = V_CalcQY(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
803	// qz = V_CalcQZ(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
804
805	MultiThread qTot = V_CalcQval(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
806	MultiThread qx = V_CalcQX(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
807	MultiThread qy = V_CalcQY(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
808	MultiThread qz = V_CalcQZ(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
809
810	return(0)
811	End
812
813
814	//function to calculate the overall q-value, given all of the necesary trig inputs
815	//
816	// (DONE)
817	// x- verify the calculation (accuracy - in all input conditions)
818	// x- verify the units of everything here, it's currently all jumbled and wrong... and repeated
819	// x- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
820	// to each pixel
821	//
822	//sdd is in [cm]
823	// distX and distY are in [mm]
824	//wavelength is in Angstroms
825	//
826	//returned magnitude of Q is in 1/Angstroms
827	//
828	ThreadSafe Function V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
829	Variable xaxval,yaxval,xctr,yctr,sdd,lam
830	Wave distX,distY
831
832	Variable dx,dy,qval,two_theta,dist
833
834
835	dx = (distX[xaxval][yaxval] - xctr) //delta x in mm
836	dy = (distY[xaxval][yaxval] - yctr) //delta y in mm
837	dist = sqrt(dx^2 + dy^2)
838
839	dist /= 10 // convert mm to cm
840
841	two_theta = atan(dist/sdd)
842
843	qval = 4Pi/lamsin(two_theta/2)
844
845	return qval
846	End
847
848	//calculates just the q-value in the x-direction on the detector
849	// (DONE)
850	// x- verify the calculation (accuracy - in all input conditions)
851	// x- verify the units of everything here, it's currently all jumbled and wrong... and repeated
852	// x- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
853	// to each pixel
854	//
855	//
856	// this properly accounts for qz
857	//
858	ThreadSafe Function V_CalcQX(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
859	Variable xaxval,yaxval,xctr,yctr,sdd,lam
860	Wave distX,distY
861
862	Variable qx,qval,phi,dx,dy,dist,two_theta
863
864	qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
865
866
867	dx = (distX[xaxval][yaxval] - xctr) //delta x in mm
868	dy = (distY[xaxval][yaxval] - yctr) //delta y in mm
869	phi = V_FindPhi(dx,dy)
870
871	//get scattering angle to project onto flat detector => Qr = qval*cos(theta)
872	dist = sqrt(dx^2 + dy^2)
873	dist /= 10 // convert mm to cm
874
875	two_theta = atan(dist/sdd)
876
877	qx = qvalcos(two_theta/2)cos(phi)
878
879	return qx
880	End
881
882	//calculates just the q-value in the y-direction on the detector
883	// (DONE)
884	// x- verify the calculation (accuracy - in all input conditions)
885	// x- verify the units of everything here, it's currently all jumbled and wrong... and repeated
886	// x- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
887	// to each pixel
888	//
889	//
890	// this properly accounts for qz
891	//
892	ThreadSafe Function V_CalcQY(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
893	Variable xaxval,yaxval,xctr,yctr,sdd,lam
894	Wave distX,distY
895
896	Variable qy,qval,phi,dx,dy,dist,two_theta
897
898	qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
899
900
901	dx = (distX[xaxval][yaxval] - xctr) //delta x in mm
902	dy = (distY[xaxval][yaxval] - yctr) //delta y in mm
903	phi = V_FindPhi(dx,dy)
904
905	//get scattering angle to project onto flat detector => Qr = qval*cos(theta)
906	dist = sqrt(dx^2 + dy^2)
907	dist /= 10 // convert mm to cm
908
909	two_theta = atan(dist/sdd)
910
911	qy = qvalcos(two_theta/2)sin(phi)
912
913	return qy
914	End
915
916	//calculates just the q-value in the z-direction on the detector
917	// (DONE)
918	// x- verify the calculation (accuracy - in all input conditions)
919	// x- verify the units of everything here, it's currently all jumbled and wrong... and repeated
920	// x- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
921	// to each pixel
922	//
923	// not actually used for any calculations, but here for completeness if anyone asks, or for 2D data export
924	//
925	// this properly accounts for qz, because it is qz
926	//
927	ThreadSafe Function V_CalcQZ(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
928	Variable xaxval,yaxval,xctr,yctr,sdd,lam
929	Wave distX,distY
930
931	Variable qz,qval,phi,dx,dy,dist,two_theta
932
933	qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
934
935
936	dx = (distX[xaxval][yaxval] - xctr) //delta x in mm
937	dy = (distY[xaxval][yaxval] - yctr) //delta y in mm
938
939	//get scattering angle to project onto flat detector => Qr = qval*cos(theta)
940	dist = sqrt(dx^2 + dy^2)
941	dist /= 10 // convert mm to cm
942
943	two_theta = atan(dist/sdd)
944
945	qz = qval*sin(two_theta/2)
946
947	return qz
948	End
949
950
951	//
952	// (DONE)
953	// x- VERIFY calculations
954	// x- This is the actual solid angle per pixel, not a ratio vs. some "unit SA"
955	// Do I just correct for the different area vs. the "nominal" central area?
956	// x- decide how to implement - YES - directly change the data values (as was done in the past)
957	// or (NOT done this way...use this as a weighting for when the data is binned to I(q). In the second method, 2D data
958	// would need this to be applied before exporting)
959	// x- do I keep a wave note indicating that this correction has been applied to the data
960	// so that it can be "un-applied"? NO
961	// x- do I calculate theta from geometry directly, or get it from Q (Assuming it's present?)
962	// (YES just from geometry, since I need SDD and dx and dy values...)
963	//
964	// "B" is passed, so I need to check for "B" and for panel orientation
965	// -if it is detector "B" (not tubes), then the normal solid angle correction applies
966	// -if it is a tube panel, then I need to know the orientation, to know which angles
967	// and pixel dimensions to use
968	//
969	// *** UPDATED 1/2020 SRK
970	// -using new calculation since the lateral direction of the tubes does not affect the solid angle
971	// projection (see He (2015) and John's memo)
972	//
973	//
974	Function V_SolidAngleCorrection(w,w_err,fname,detStr,destPath)
975	Wave w,w_err
976	String fname,detStr,destPath
977
978	Variable sdd,xCtr,yCtr,lambda
979	String orientation
980
981	// get all of the geometry information
982	orientation = V_getDet_tubeOrientation(fname,detStr)
983	sdd = V_getDet_ActualDistance(fname,detStr)
984
985	// this is ctr in mm
986	xCtr = V_getDet_beam_center_x_mm(fname,detStr)
987	yCtr = V_getDet_beam_center_y_mm(fname,detStr)
988	lambda = V_getWavelength(fname)
989
990	SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr)
991
992	Wave data_realDistX = data_realDistX
993	Wave data_realDistY = data_realDistY
994
995	Duplicate/O w solid_angle,tmp_theta,tmp_dist,tmp_theta_i //in the current df
996
997	//// calculate the scattering angle
998	// dx = (distX - xctr) //delta x in mm
999	// dy = (distY - yctr) //delta y in mm
1000	tmp_dist = sqrt((data_realDistX - xctr)^2 + (data_realDistY - yctr)^2)
1001
1002	tmp_dist /= 10 // convert mm to cm
1003	// sdd is in [cm]
1004
1005	tmp_theta = atan(tmp_dist/sdd) //this is two_theta, the (total) scattering angle
1006
1007	Variable ii,jj,numx,numy,dx,dy
1008	numx = DimSize(tmp_theta,0)
1009	numy = DimSize(tmp_theta,1)
1010
1011	if(cmpstr(detStr,"B")==0)
1012	//detector B is a grid, straightforward cos^3 solid angle
1013	for(ii=0 ;ii<numx;ii+=1)
1014	for(jj=0;jj<numy;jj+=1)
1015
1016	if(ii==0) //do a forward difference if ii==0
1017	dx = (data_realDistX[ii+1][jj] - data_realDistX[ii][jj]) //delta x for the pixel
1018	else
1019	dx = (data_realDistX[ii][jj] - data_realDistX[ii-1][jj]) //delta x for the pixel
1020	endif
1021
1022
1023	if(jj==0)
1024	dy = (data_realDistY[ii][jj+1] - data_realDistY[ii][jj]) //delta y for the pixel
1025	else
1026	dy = (data_realDistY[ii][jj] - data_realDistY[ii][jj-1]) //delta y for the pixel
1027	endif
1028
1029	dx /= 10
1030	dy /= 10 // convert mm to cm (since sdd is in cm)
1031	solid_angle[ii][jj] = dx*dy //this is in cm^2
1032	endfor
1033	endfor
1034
1035	// to cover up any issues w/negative dx or dy
1036	solid_angle = abs(solid_angle)
1037
1038	// solid_angle correction
1039	// == dxdycos^3/sdd^2
1040	solid_angle *= (cos(tmp_theta))^3
1041	solid_angle /= sdd^2
1042
1043	// Here it is! Apply the correction to the intensity (I divide -- to get the counts per solid angle!!)
1044	w /= solid_angle
1045
1046	// correctly apply the correction to the error wave (assume a perfect value?)
1047	w_err /= solid_angle //
1048
1049	else
1050	//
1051	//different calculation for the tubes, different calculation based on XY orientation
1052	//
1053	if(cmpstr(orientation,"vertical")==0)
1054	// L/R panels, tube axis is y-direction
1055	// this is now a different tmp_dist
1056	// convert everything to cm first!
1057	// sdd is in [cm], everything else is in [mm]
1058	tmp_dist = (data_realDistY/10 - yctr/10)/sqrt((data_realDistX/10 - xctr/10)^2 + sdd^2)
1059	tmp_theta_i = atan(tmp_dist) //this is theta_y
1060
1061	else
1062	// horizontal orientation (T/B panels)
1063	// this is now a different tmp_dist
1064	// convert everything to cm first!
1065	// sdd is in [cm], everything else is in [mm]
1066	tmp_dist = (data_realDistX/10 - xctr/10)/sqrt((data_realDistY/10 - yctr/10)^2 + sdd^2)
1067	tmp_theta_i = atan(tmp_dist) //this is theta_x
1068
1069	endif
1070
1071	for(ii=0 ;ii<numx;ii+=1)
1072	for(jj=0;jj<numy;jj+=1)
1073
1074	if(ii==0) //do a forward difference if ii==0
1075	dx = (data_realDistX[ii+1][jj] - data_realDistX[ii][jj]) //delta x for the pixel
1076	else
1077	dx = (data_realDistX[ii][jj] - data_realDistX[ii-1][jj]) //delta x for the pixel
1078	endif
1079
1080
1081	if(jj==0)
1082	dy = (data_realDistY[ii][jj+1] - data_realDistY[ii][jj]) //delta y for the pixel
1083	else
1084	dy = (data_realDistY[ii][jj] - data_realDistY[ii][jj-1]) //delta y for the pixel
1085	endif
1086
1087	dx /= 10
1088	dy /= 10 // convert mm to cm (since sdd is in cm)
1089	solid_angle[ii][jj] = dx*dy //this is in cm^2
1090	endfor
1091	endfor
1092
1093	// to cover up any issues w/negative dx or dy
1094	solid_angle = abs(solid_angle)
1095
1096	// solid_angle correction
1097	// == dxdycos(th)^2*cos(th_i)/sdd^2 using either the theta_x or theta_y value
1098	solid_angle = (cos(tmp_theta))^2cos(tmp_theta_i)
1099	solid_angle /= sdd^2
1100
1101	// Here it is! Apply the correction to the intensity (I divide -- to get the counts per solid angle!!)
1102	w /= solid_angle
1103
1104	//
1105	// correctly apply the correction to the error wave (assume a perfect value?)
1106	w_err /= solid_angle //
1107
1108	endif
1109
1110
1111	// DONE x- clean up after I'm satisfied computations are correct
1112	KillWaves/Z tmp_theta,tmp_dist,tmp_theta_i
1113
1114	return(0)
1115	end
1116
1117	// this is the incorrect solid angle correction that does not take into
1118	// account the tube geometry. It is correct for the high-res detector (and the 30m Ordela)
1119	//
1120	// -- only for testing to prove that the cos(th)^2 *cos(th_i) is correct
1121	//
1122	Function V_SolidAngleCorrection_COS3(w,w_err,fname,detStr,destPath)
1123	Wave w,w_err
1124	String fname,detStr,destPath
1125
1126	Variable sdd,xCtr,yCtr,lambda
1127	String orientation
1128
1129	// get all of the geometry information
1130	orientation = V_getDet_tubeOrientation(fname,detStr)
1131	sdd = V_getDet_ActualDistance(fname,detStr)
1132
1133	// this is ctr in mm
1134	xCtr = V_getDet_beam_center_x_mm(fname,detStr)
1135	yCtr = V_getDet_beam_center_y_mm(fname,detStr)
1136	lambda = V_getWavelength(fname)
1137
1138	SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr)
1139
1140	Wave data_realDistX = data_realDistX
1141	Wave data_realDistY = data_realDistY
1142
1143	Duplicate/O w solid_angle,tmp_theta,tmp_dist //in the current df
1144
1145	//// calculate the scattering angle
1146	// dx = (distX - xctr) //delta x in mm
1147	// dy = (distY - yctr) //delta y in mm
1148	tmp_dist = sqrt((data_realDistX - xctr)^2 + (data_realDistY - yctr)^2)
1149
1150	tmp_dist /= 10 // convert mm to cm
1151	// sdd is in [cm]
1152
1153	tmp_theta = atan(tmp_dist/sdd) //this is two_theta, the (total) scattering angle
1154
1155	Variable ii,jj,numx,numy,dx,dy
1156	numx = DimSize(tmp_theta,0)
1157	numy = DimSize(tmp_theta,1)
1158
1159	for(ii=0 ;ii<numx;ii+=1)
1160	for(jj=0;jj<numy;jj+=1)
1161
1162	if(ii==0) //do a forward difference if ii==0
1163	dx = (data_realDistX[ii+1][jj] - data_realDistX[ii][jj]) //delta x for the pixel
1164	else
1165	dx = (data_realDistX[ii][jj] - data_realDistX[ii-1][jj]) //delta x for the pixel
1166	endif
1167
1168
1169	if(jj==0)
1170	dy = (data_realDistY[ii][jj+1] - data_realDistY[ii][jj]) //delta y for the pixel
1171	else
1172	dy = (data_realDistY[ii][jj] - data_realDistY[ii][jj-1]) //delta y for the pixel
1173	endif
1174
1175	dx /= 10
1176	dy /= 10 // convert mm to cm (since sdd is in cm)
1177	solid_angle[ii][jj] = dx*dy //this is in cm^2
1178	endfor
1179	endfor
1180
1181	// to cover up any issues w/negative dx or dy
1182	solid_angle = abs(solid_angle)
1183
1184	// solid_angle correction
1185	// == dxdycos^3/sdd^2
1186	solid_angle *= (cos(tmp_theta))^3
1187	solid_angle /= sdd^2
1188
1189	// Here it is! Apply the correction to the intensity (I divide -- to get the counts per solid angle!!)
1190	w /= solid_angle
1191
1192	// correctly apply the correction to the error wave (assume a perfect value?)
1193	w_err /= solid_angle //
1194
1195
1196	// DONE x- clean up after I'm satisfied computations are correct
1197	KillWaves/Z tmp_theta,tmp_dist,tmp_theta_i
1198
1199	return(0)
1200	end
1201
1202
1203
1204	//
1205	// Large angle transmission correction
1206	//
1207	// DIVIDE the intensity by this correction to get the right answer
1208	//
1209	//
1210	// Apply the large angle transmssion correction as the data is converted to WORK
1211	// so that whether the data is saved as 2D or 1D, the correction has properly been done.
1212	//
1213	// This is, however, a SAMPLE dependent calculation, not purely instrument geometry.
1214	//
1215	Function V_LargeAngleTransmissionCorr(w,w_err,fname,detStr,destPath)
1216	Wave w,w_err
1217	String fname,detStr,destPath
1218
1219	Variable sdd,xCtr,yCtr,trans,trans_err,uval
1220
1221	// get all of the geometry information
1222	// orientation = V_getDet_tubeOrientation(fname,detStr)
1223	sdd = V_getDet_ActualDistance(fname,detStr)
1224
1225	// this is ctr in mm
1226	xCtr = V_getDet_beam_center_x_mm(fname,detStr)
1227	yCtr = V_getDet_beam_center_y_mm(fname,detStr)
1228	trans = V_getSampleTransmission(fname)
1229	trans_err = V_getSampleTransError(fname)
1230
1231	SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr)
1232
1233	Wave data_realDistX = data_realDistX
1234	Wave data_realDistY = data_realDistY
1235
1236	Duplicate/O w lat_corr,tmp_theta,tmp_dist,lat_err,tmp_err //in the current df
1237
1238	//// calculate the scattering angle
1239	// dx = (distX - xctr) //delta x in mm
1240	// dy = (distY - yctr) //delta y in mm
1241	tmp_dist = sqrt((data_realDistX - xctr)^2 + (data_realDistY - yctr)^2)
1242
1243	tmp_dist /= 10 // convert mm to cm
1244	// sdd is in [cm]
1245
1246	tmp_theta = atan(tmp_dist/sdd) //this is two_theta, the scattering angle
1247
1248	Variable ii,jj,numx,numy,dx,dy,cos_th,arg,tmp
1249	numx = DimSize(tmp_theta,0)
1250	numy = DimSize(tmp_theta,1)
1251
1252
1253	//optical thickness
1254	uval = -ln(trans) //use natural logarithm
1255
1256	for(ii=0 ;ii<numx;ii+=1)
1257	for(jj=0;jj<numy;jj+=1)
1258
1259	cos_th = cos(tmp_theta[ii][jj])
1260	arg = (1-cos_th)/cos_th
1261
1262	// a Taylor series around uval*arg=0 only needs about 4 terms for very good accuracy
1263	// correction= 1 - 0.5uvalarg + (uvalarg)^2/6 - (uvalarg)^3/24 + (uval*arg)^4/120
1264	// OR
1265	if((uval<0.01) \|\| (cos_th>0.99))
1266	//small arg, approx correction
1267	lat_corr[ii][jj] = 1-0.5uvalarg
1268	else
1269	//large arg, exact correction
1270	lat_corr[ii][jj] = (1-exp(-uvalarg))/(uvalarg)
1271	endif
1272
1273	// (DONE)
1274	// x- properly calculate and apply the 2D error propagation
1275	if(trans == 1)
1276	lat_err[ii][jj] = 0 //no correction, no error
1277	else
1278	//sigT, calculated from the Taylor expansion
1279	tmp = (1/trans)(arg/2-arg^2/3uval+arg^3/8uval^2-arg^4/30uval^3)
1280	tmp *= tmp
1281	tmp *= trans_err^2
1282	tmp = sqrt(tmp) //sigT
1283
1284	lat_err[ii][jj] = tmp
1285	endif
1286
1287
1288	endfor
1289	endfor
1290
1291
1292
1293	// Here it is! Apply the correction to the intensity (divide -- to get the proper correction)
1294	w /= lat_corr
1295
1296	// relative errors add in quadrature to the current 2D error
1297	tmp_err = (w_err/lat_corr)^2 + (lat_err/lat_corr)^2ww/lat_corr^2
1298	tmp_err = sqrt(tmp_err)
1299
1300	w_err = tmp_err
1301
1302
1303	// DONE x- clean up after I'm satisfied computations are correct
1304	KillWaves/Z tmp_theta,tmp_dist,tmp_err,lat_err
1305
1306	return(0)
1307	end
1308
1309
1310
1311	//
1312	//test procedure, not called anymore
1313	Proc V_AbsoluteScaling(type,c0,c1,c2,c3,c4,c5,I_err)
1314	String type
1315	Variable c0=1,c1=0.1,c2=0.95,c3=0.1,c4=1,c5=32.0,I_err=0.32
1316	Prompt type,"WORK data type",popup,"CAL;COR;SAM"
1317	Prompt c0, "Sample Transmission"
1318	Prompt c1, "Sample Thickness (cm)"
1319	Prompt c2, "Standard Transmission"
1320	Prompt c3, "Standard Thickness (cm)"
1321	Prompt c4, "I(0) from standard fit (normalized to 1E8 monitor cts)"
1322	Prompt c5, "Standard Cross-Section (cm-1)"
1323	Prompt I_err, "error in I(q=0) (one std dev)"
1324
1325	Variable err
1326	//call the function to do the math
1327	//data from "type" will be scaled and deposited in ABS
1328	err = V_Absolute_Scale(type,c0,c1,c2,c3,c4,c5,I_err)
1329
1330	if(err)
1331	Abort "Error in V_Absolute_Scale()"
1332	endif
1333
1334	//contents are always dumped to ABS
1335	type = "ABS"
1336
1337	//need to update the display with "data" from the correct dataFolder
1338	//reset the current display type to "type"
1339	SVAR gCurDispType = root:Packages:NIST:VSANS:Globals:gCurDispType
1340	gCurDispType = Type
1341
1342	V_UpdateDisplayInformation(Type)
1343
1344	End
1345
1346	//
1347	//
1348	// kappa comes in as s_izero, so be sure to use 1/kappa_err
1349	//
1350	//convert the "type" data to absolute scale using the given standard information
1351	//s_ is the standard
1352	//w_ is the "work" file
1353	//both are work files and should already be normalized to 10^8 monitor counts
1354	Function V_Absolute_Scale(type,absStr)
1355	String type,absStr
1356
1357
1358	Variable w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err
1359
1360	Variable defmon = 1e8,w_moncount,s1,s2,s3,s4
1361	Variable scale,trans_err
1362	Variable err,ii
1363	String detStr
1364
1365	// be sure that the starting data exists
1366	err = V_WorkDataExists(type)
1367	if(err==1)
1368	return(err)
1369	endif
1370
1371	//copy from current dir (type) to ABS
1372	V_CopyHDFToWorkFolder(type,"ABS")
1373
1374	// TODO: -- which monitor to use? Here, I think it should already be normalized to 10^8
1375	//
1376	// w_moncount = V_getMonitorCount(type) //monitor count in "type"
1377
1378	w_moncount = V_getBeamMonNormData(type)
1379
1380	if(w_moncount == 0)
1381	//zero monitor counts will give divide by zero ---
1382	DoAlert 0,"Total monitor count in data file is zero. No rescaling of data"
1383	Return(1) //report error
1384	Endif
1385
1386	w_trans = V_getSampleTransmission(type) //sample transmission
1387	w_thick = V_getSampleThickness(type) //sample thickness
1388	trans_err = V_getSampleTransError(type)
1389
1390
1391	//get the parames from the list
1392	s_trans = NumberByKey("TSTAND", absStr, "=", ";") //parse the list of values
1393	s_thick = NumberByKey("DSTAND", absStr, "=", ";")
1394	s_izero = NumberByKey("IZERO", absStr, "=", ";")
1395	s_cross = NumberByKey("XSECT", absStr, "=", ";")
1396	kappa_err = NumberByKey("SDEV", absStr, "=", ";")
1397
1398
1399	//calculate scale factor
1400	s1 = defmon/w_moncount // monitor count (s1 should be 1)
1401	s2 = s_thick/w_thick
1402	s3 = s_trans/w_trans
1403	s4 = s_cross/s_izero
1404	scale = s1s2s3*s4
1405
1406
1407	// kappa comes in as s_izero, so be sure to use 1/kappa_err
1408
1409	// and now loop through all of the detectors
1410	//do the actual absolute scaling here, modifying the data in ABS
1411	for(ii=0;ii<ItemsInList(ksDetectorListNoB);ii+=1)
1412	detStr = StringFromList(ii, ksDetectorListNoB, ";")
1413	Wave data = V_getDetectorDataW("ABS",detStr)
1414	Wave data_err = V_getDetectorDataErrW("ABS",detStr)
1415
1416	data *= scale
1417	data_err = sqrt(scale^2data_err^2 + scale^2data^2*(kappa_err^2/s_izero^2 +trans_err^2/w_trans^2))
1418	endfor
1419
1420	// do the back detector separately, if it is set to be used
1421	NVAR gIgnoreDetB = root:Packages:NIST:VSANS:Globals:gIgnoreDetB
1422	if(gIgnoreDetB == 0)
1423	detStr = "B"
1424	Wave data = V_getDetectorDataW("ABS",detStr)
1425	Wave data_err = V_getDetectorDataErrW("ABS",detStr)
1426
1427	//get the parames from the list
1428	s_trans = NumberByKey("TSTAND_B", absStr, "=", ";") //parse the list of values
1429	s_thick = NumberByKey("DSTAND_B", absStr, "=", ";")
1430	s_izero = NumberByKey("IZERO_B", absStr, "=", ";")
1431	s_cross = NumberByKey("XSECT_B", absStr, "=", ";")
1432	kappa_err = NumberByKey("SDEV_B", absStr, "=", ";")
1433
1434	//calculate scale factor
1435	s1 = defmon/w_moncount // monitor count (s1 should be 1)
1436	s2 = s_thick/w_thick
1437	s3 = s_trans/w_trans
1438	s4 = s_cross/s_izero
1439	scale = s1s2s3*s4
1440
1441	data *= scale
1442	data_err = sqrt(scale^2data_err^2 + scale^2data^2*(kappa_err^2/s_izero^2 +trans_err^2/w_trans^2))
1443	endif
1444
1445	//********* 15APR02
1446	// DO NOT correct for atenuators here - the COR step already does this, putting all of the data on equal
1447	// footing (zero atten) before doing the subtraction.
1448
1449	Return (0) //no error
1450	End
1451
1452
1453	//
1454	// TODO:
1455	// -- DoAlert 0,"This has not yet been updated for VSANS"
1456	//
1457	//
1458	// match the attenuation of the RAW data to the "type" data
1459	// so that they can be properly added
1460	//
1461	// are the attenuator numbers the same? if so exit
1462	//
1463	// if not, find the attenuator number for type
1464	// - find both attenuation factors
1465	//
1466	// rescale the raw data to match the ratio of the two attenuation factors
1467	// -- adjust the detector count (rw)
1468	// -- the linear data
1469	//
1470	//
1471	Function V_Adjust_RAW_Attenuation(type)
1472	String type
1473
1474	DoAlert 0,"This has not yet been updated for VSANS"
1475
1476	WAVE rw=$("root:Packages:NIST:RAW:realsread")
1477	WAVE linear_data=$("root:Packages:NIST:RAW:linear_data")
1478	WAVE data=$("root:Packages:NIST:RAW:data")
1479	WAVE data_err=$("root:Packages:NIST:RAW:linear_data_error")
1480	WAVE/T tw = $("root:Packages:NIST:RAW:textRead")
1481
1482	WAVE dest_reals=$("root:Packages:NIST:"+type+":realsread")
1483
1484	Variable dest_atten,raw_atten,tol
1485	Variable lambda,raw_atten_err,raw_AttenFactor,dest_attenFactor,dest_atten_err
1486	String fileStr
1487
1488	dest_atten = dest_reals[3]
1489	raw_atten = rw[3]
1490
1491	tol = 0.1 // within 0.1 atten units is OK
1492	if(abs(dest_atten - raw_atten) < tol )
1493	return(0)
1494	endif
1495
1496	fileStr = tw[3]
1497	lambda = rw[26]
1498	// TODO access correct values
1499	raw_AttenFactor = 1//AttenuationFactor(fileStr,lambda,raw_atten,raw_atten_err)
1500	dest_AttenFactor = 1//AttenuationFactor(fileStr,lambda,dest_atten,dest_atten_err)
1501
1502	rw[2] *= dest_AttenFactor/raw_AttenFactor
1503	linear_data *= dest_AttenFactor/raw_AttenFactor
1504
1505	// to keep "data" and linear_data in sync
1506	data = linear_data
1507
1508	return(0)
1509	End
1510
1511	//
1512	// testing procedure, called from a menu selection
1513	//
1514	Proc V_DIV_a_Workfile(type)
1515	String type
1516	Prompt type,"WORK data type",popup,"SAM;EMP;BGD;ADJ;"
1517
1518	//macro will take whatever is in SELECTED folder and DIVide it by the current
1519	//contents of the DIV folder - the function will check for existence
1520	//before proceeding
1521
1522	Abort "This has not yet been updated for VSANS"
1523
1524	Variable err
1525	err = V_DIVCorrection(type) //returns err = 1 if data doesn't exist in specified folders
1526
1527	if(err)
1528	Abort "error in V_DIVCorrection()"
1529	endif
1530
1531	//contents are NOT always dumped to CAL, but are in the new type folder
1532
1533	String newTitle = "WORK_"+type
1534	DoWindow/F VSANS_Data
1535	DoWindow/T VSANS_Data, newTitle
1536	KillStrings/Z newTitle
1537
1538	//need to update the display with "data" from the correct dataFolder
1539	//reset the current displaytype to "type"
1540	String/G root:Packages:NIST:VSANS:Globals:gCurDispType=Type
1541
1542	V_UpdateDisplayInformation(type)
1543
1544	End
1545
1546
1547	//
1548	// TODO:
1549	// x- DoAlert 0,"This has not yet been updated for VSANS"
1550	// -- how is the error propagation handled? Be sure it is calculated correctly when DIV is generated
1551	// and is applied correctly here...
1552	//
1553	//function will divide the contents of "workType" folder with the contents of
1554	//the DIV folder + detStr
1555	// all data is linear scale for the calculation
1556	//
1557	Function V_DIVCorrection(data,data_err,detStr,workType)
1558	Wave data,data_err
1559	String detStr,workType
1560
1561	//check for existence of data in type and DIV
1562	// if the desired data doesn't exist, let the user know, and abort
1563	String destPath=""
1564
1565	NVAR gIgnoreDetB = root:Packages:NIST:VSANS:Globals:gIgnoreDetB
1566	if(cmpstr(detStr,"B")==0 && gIgnoreDetB)
1567	return(0)
1568	endif
1569
1570
1571	if(WaveExists(data) == 0)
1572	Print "The data wave does not exist in V_DIVCorrection()"
1573	Return(1) //error condition
1574	Endif
1575
1576	//check for DIV
1577	// if the DIV workfile doesn't exist, let the user know,and abort
1578	// !! be sure to check first, before trying to access the wave
1579
1580	// WAVE/Z div_data = $("root:Packages:NIST:VSANS:DIV:entry:instrument:detector_"+detStr+":data")
1581	if(WaveExists($("root:Packages:NIST:VSANS:DIV:entry:instrument:detector_"+detStr+":data")) == 0)
1582	Print "The DIV wave does not exist in V_DIVCorrection()"
1583	Return(1) //error condition
1584	Endif
1585	if(WaveExists($("root:Packages:NIST:VSANS:DIV:entry:instrument:detector_"+detStr+":linear_data_error")) == 0)
1586	Print "The DIV error wave does not exist in V_DIVCorrection()"
1587	Return(1) //error condition
1588	Endif
1589	//files exist, proceed
1590
1591	WAVE/Z div_data_err = V_getDetectorDataErrW("DIV",detStr)
1592	WAVE/Z div_data = V_getDetectorDataW("DIV",detStr)
1593
1594
1595
1596	// do the error propagation first, since data is changed by the correction
1597	data_err = sqrt(data_err^2/div_data^2 + div_data_err^2 * data^2/div_data^4 )
1598
1599	// then the correction
1600	data /= div_data
1601
1602
1603	Return(0)
1604	End
1605
1606
1607	//////////////////////////
1608	// detector corrections to stitch the back detector into one proper image
1609	//
1610	//
1611	//
1612
1613
1614	//
1615	// to register the image on the back detector panel
1616	//
1617	// middle portion (552 pix in Y) is held fixed
1618	// top portion of image is shifted right and down
1619	// bottom portion of image is shifted right and up
1620	//
1621	// remainder of image is filled with Zero (NaN causes problems converting to WORK)
1622	//
1623	// currently, data is not added together and averaged, but it could be
1624	//
1625	Function V_ShiftBackDetImage(w,adjW)
1626	Wave w,adjW
1627
1628	NVAR gHighResBinning = root:Packages:NIST:VSANS:Globals:gHighResBinning
1629
1630	// this is necessary for some old data with the 150x150 back (dummy) panel
1631	// the proper back detector has an x-dimension of 680 pixels. Don't do the shift
1632	// if the dimensions are incorrect.
1633	if(DimSize(w,0) < 680)
1634	adjW=w
1635	return(0)
1636	endif
1637
1638	adjW=0
1639
1640	Variable topX,bottomX
1641	Variable topY,bottomY
1642	Variable totalY,ccdX,ccdY
1643
1644	// topX = 7
1645	// topY = 105
1646
1647	// bottomX = 5
1648	// bottomY = 35
1649
1650	// TODOHIGHRES
1651	// the detector pix dimensions are hard-wired, be sure the are correct
1652	switch(gHighResBinning)
1653	case 1:
1654	topX = kShift_topX_bin1
1655	topY = kShift_topY_bin1
1656	bottomX = kShift_bottomX_bin1
1657	bottomY = kShift_bottomY_bin1
1658
1659	totalY = 6624 // total YDim
1660	ccdY = 2208 // = YDim/3
1661	ccdX = 2720 // = xDim
1662	break
1663	case 4:
1664	topX = kShift_topX_bin4
1665	topY = kShift_topY_bin4
1666	bottomX = kShift_bottomX_bin4
1667	bottomY = kShift_bottomY_bin4
1668
1669	totalY = 1656 // total YDim
1670	ccdY = 552 // = YDim/3
1671	ccdX = 680 // = xDim
1672
1673
1674	break
1675	default:
1676	Abort "No binning case matches in V_ShiftBackDetImage"
1677
1678	endswitch
1679
1680	// middle
1681	adjW[][ccdY,ccdY+ccdY] = w[p][q]
1682
1683	//top
1684	adjW[0+topX,ccdX-1][ccdY+ccdY,totalY-1-topY] = w[p-topX][q+topY]
1685
1686	//bottom
1687	adjW[0+bottomX,ccdX-1][0+bottomY,ccdY-1] = w[p-bottomX][q-bottomY]
1688
1689
1690	return(0)
1691	End
1692
1693
1694	Proc pV_MedianFilterBack(folder)
1695	String folder="RAW"
1696
1697	V_MedianFilterBack(folder)
1698	end
1699
1700	Function V_MedianFilterBack(folder)
1701	String folder
1702
1703	Wave w = V_getDetectorDataW(folder,"B")
1704
1705	NVAR gHighResBinning = root:Packages:NIST:VSANS:Globals:gHighResBinning
1706	switch(gHighResBinning)
1707	case 1:
1708	MatrixFilter /N=11 /P=1 median w // /P=n flag sets the number of passes (default is 1 pass)
1709
1710	Print "* median noise filter 11x11 applied to the back detector (1 pass) *"
1711	break
1712	case 4:
1713	MatrixFilter /N=3 /P=1 median w // /P=n flag sets the number of passes (default is 1 pass)
1714
1715	Print "* median noise filter 3x3 applied to the back detector (1 pass) *"
1716	break
1717	default:
1718	Abort "No binning case matches in V_MedianFilterBack"
1719	endswitch
1720
1721	return(0)
1722	End
1723
1724
1725	Proc pV_SubtractReadNoiseBack(folder,ReadNoise)
1726	String folder="RAW"
1727	Variable readNoise=3160
1728
1729	V_SubtractReadNoiseBack(folder,readNoise)
1730	end
1731
1732	Function V_SubtractReadNoiseBack(folder,readNoise)
1733	String folder
1734	Variable readNoise
1735
1736	Wave w = V_getDetectorDataW(folder,"B")
1737	w -= readNoise // a constant value
1738
1739	// MatrixFilter /N=3 median w
1740	// Print "* median noise filter applied to the back detector*"
1741
1742	return(0)
1743	End
1744
1745
1746	Proc pV_MedianAndReadNoiseBack(folder,ReadNoise)
1747	String folder="RAW"
1748	Variable readNoise=3160
1749
1750	V_MedianAndReadNoiseBack(folder,readNoise)
1751	end
1752
1753	Function V_MedianAndReadNoiseBack(folder,readNoise)
1754	String folder
1755	Variable readNoise
1756
1757	Wave w = V_getDetectorDataW(folder,"B")
1758	w -= readNoise // a constant value
1759
1760	NVAR gHighResBinning = root:Packages:NIST:VSANS:Globals:gHighResBinning
1761	switch(gHighResBinning)
1762	case 1:
1763	MatrixFilter /N=11 /P=1 median w // /P=n flag sets the number of passes (default is 1 pass)
1764
1765	Print "* median noise filter 11x11 applied to the back detector (1 pass) *"
1766	break
1767	case 4:
1768	MatrixFilter /N=3 /P=1 median w // /P=n flag sets the number of passes (default is 1 pass)
1769
1770	Print "* median noise filter 3x3 applied to the back detector (1 pass) *"
1771	break
1772	default:
1773	Abort "No binning case matches in V_MedianAndReadNoiseBack"
1774	endswitch
1775
1776	return(0)
1777	End
1778
1779
1780
1781	////////////////
1782	// Detector efficiency and shadowing
1783	///////////////
1784
1785	//
1786	// Tube efficiency + shadowing
1787	//
1788	//
1789	// -- check for the existence of the proper tables (correct wavelength)
1790	// -- generate tables if needed (one-time calculation)
1791	//
1792	// interpolate the table for the correction - to avoid repeated integration
1793	//
1794	// store the tables in: root:Packages:NIST:VSANS:Globals:Efficiency:
1795	//
1796	Function V_TubeEfficiencyShadowCorr(w,w_err,fname,detStr,destPath)
1797	Wave w,w_err
1798	String fname,detStr,destPath
1799
1800	Variable sdd,xCtr,yCtr,lambda
1801	String orientation
1802
1803	// if the panel is "B", exit - since it is not tubes, and this should not be called
1804	if(cmpstr(detStr,"B")==0)
1805	return(1)
1806	endif
1807
1808	// get all of the geometry information
1809	orientation = V_getDet_tubeOrientation(fname,detStr)
1810	sdd = V_getDet_ActualDistance(fname,detStr)
1811
1812	// this is ctr in mm
1813	xCtr = V_getDet_beam_center_x_mm(fname,detStr)
1814	yCtr = V_getDet_beam_center_y_mm(fname,detStr)
1815	lambda = V_getWavelength(fname)
1816
1817	SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr)
1818
1819	Wave data_realDistX = data_realDistX
1820	Wave data_realDistY = data_realDistY
1821
1822	Duplicate/O w tmp_theta_x,tmp_theta_y,tmp_dist,tmp_corr //in the current df
1823
1824	//// calculate the scattering angles theta_x and theta_y
1825
1826	// flip the definitions of x and y for the T/B panels so that x is always lateral WRT the tubes
1827	// and y is always along the length of the tubes
1828
1829	if(cmpstr(orientation,"vertical")==0)
1830	// L/R panels, tube axis is y-direction
1831	// this is now a different tmp_dist
1832	// convert everything to cm first!
1833	// sdd is in [cm], everything else is in [mm]
1834	tmp_dist = (data_realDistY/10 - yctr/10)/sqrt((data_realDistX/10 - xctr/10)^2 + sdd^2)
1835	tmp_theta_y = atan(tmp_dist) //this is theta_y
1836	tmp_theta_x = atan( (data_realDistX/10 - xctr/10)/sdd )
1837
1838	else
1839	// horizontal orientation (T/B panels)
1840	// this is now a different tmp_dist
1841	// convert everything to cm first!
1842	// sdd is in [cm], everything else is in [mm]
1843	tmp_dist = (data_realDistX/10 - xctr/10)/sqrt((data_realDistY/10 - yctr/10)^2 + sdd^2)
1844	tmp_theta_y = atan(tmp_dist) //this is theta_y, along tube direction
1845	tmp_theta_x = atan( (data_realDistY/10 - yctr/10)/sdd ) // this is laterally across tubes
1846	endif
1847
1848
1849	// identify if the 2D efficiency wave has been generated for the data wavelength
1850	//
1851	// if so, declare
1852	// if not, generate
1853
1854	if(WaveExists($"root:Packages:NIST:VSANS:Globals:Efficiency:eff") == 0)
1855	// generate the proper efficiency wave, at lambda
1856	NewDataFolder/O root:Packages:NIST:VSANS:Globals:Efficiency
1857	Print "recalculating efficiency table ..."
1858	V_TubeShadowEfficiencyTables_oneLam(lambda)
1859	// declare the wave
1860	Wave/Z effW = root:Packages:NIST:VSANS:Globals:Efficiency:eff
1861	else
1862	Wave/Z effW = root:Packages:NIST:VSANS:Globals:Efficiency:eff
1863	//is the efficiency at the correct wavelength?
1864	string str=note(effW)
1865	// Print "Note = ",str
1866
1867	if(V_CloseEnough(lambda,NumberByKey("LAMBDA", str,"="),0.1)) //absolute difference of < 0.1 A
1868	// yes, proceed, no need to do anything
1869	else
1870	// no, regenerate the efficiency and then proceed (wave already declared)
1871	Print "recalculating efficiency table ..."
1872	V_TubeShadowEfficiencyTables_oneLam(lambda)
1873	endif
1874	endif
1875
1876
1877	Variable ii,jj,numx,numy,xAngle,yAngle
1878	numx = DimSize(w,0)
1879	numy = DimSize(w,1)
1880
1881	// loop over all of the pixels of the panel and find the interpolated correction (save as a wave)
1882	//
1883	for(ii=0 ;ii<numx;ii+=1)
1884	for(jj=0;jj<numy;jj+=1)
1885
1886	// from the angles, find the (x,y) point to interpolate to get the efficiency
1887
1888	xAngle = tmp_theta_x[ii][jj]
1889	yAngle = tmp_theta_y[ii][jj]
1890
1891	xAngle = abs(xAngle)
1892	yAngle = abs(yAngle)
1893
1894	// the x and y scaling of the eff wave (2D) was set when it was generated (in radians)
1895	// simply reading the scaled xy value does not interpolate!!
1896	// tmp_corr[ii][jj] = effW(xAngle)(yAngle) // NO, returns "stepped" values
1897	tmp_corr[ii][jj] = Interp2D(effW,xAngle,yAngle)
1898
1899	endfor
1900	endfor
1901	//
1902	//
1903	// apply the correction and calculate the error
1904	//
1905	// Here it is! Apply the correction to the intensity (divide -- to get the proper correction)
1906	w /= tmp_corr
1907	//
1908	// relative errors add in quadrature to the current 2D error
1909	// assume that this numerical calculation of efficiency is exact
1910	//
1911	// tmp_err = (w_err/tmp_corr)^2 + (lat_err/lat_corr)^2ww/lat_corr^2
1912	// tmp_err = sqrt(tmp_err)
1913	//
1914	// w_err = tmp_err
1915	//
1916
1917	// TODO
1918	// - clean up after I'm satisfied computations are correct
1919	// KillWaves/Z tmp_theta_x,tmp_theta_y,tmp_dist,tmp_err,tmp_corr
1920
1921	return(0)
1922	end
1923
1924
1925
1926	// the actual integration of the efficiency for an individual pixel
1927	Function V_Efficiency_Integral(pWave,in_u)
1928	Wave pWave
1929	Variable in_u
1930
1931	Variable lambda,th_x,th_y,u_p,integrand,T_sh,max_x,d_ss,d_He
1932
1933	lambda = pWave[0]
1934	th_x = pWave[1]
1935	th_y = pWave[2]
1936
1937	u_p = in_u + kTube_cc * cos(th_x)
1938
1939	// calculate shadow if th_x > 23.727 deg. th_x is input in radians
1940	max_x = 23.727 / 360 * 2*pi
1941	if(th_x < max_x)
1942	T_sh = 1
1943	else
1944
1945	// get d_ss
1946	if(abs(u_p) < kTube_ri)
1947	d_ss = sqrt( (kTube_ri + kTube_ss)^2 - in_u^2 ) - sqrt( kTube_ri^2 - in_u^2 )
1948	elseif (abs(u_p) < (kTube_ri + kTube_ss))
1949	d_ss = sqrt( (kTube_ri + kTube_ss)^2 - in_u^2 )
1950	else
1951	d_ss = 0
1952	endif
1953
1954	// get d_He
1955	if(abs(u_p) < kTube_ri)
1956	d_He = 2 * sqrt( kTube_ri^2 - in_u^2 )
1957	else
1958	d_He = 0
1959	endif
1960
1961	//calculate T_sh
1962	T_sh = exp(-2kSig_sslambdad_ss/cos(th_y)) exp(-kSig_8b_Helambdad_He/cos(th_y))
1963
1964	endif
1965
1966
1967	// calculate the integrand
1968
1969	//note that the in_u value is used here to find d_ss and d_he (not u_p)
1970	// get d_ss
1971	if(abs(in_u) < kTube_ri)
1972	d_ss = sqrt( (kTube_ri + kTube_ss)^2 - in_u^2 ) - sqrt( kTube_ri^2 - in_u^2 )
1973	elseif (abs(in_u) < (kTube_ri + kTube_ss))
1974	d_ss = sqrt( (kTube_ri + kTube_ss)^2 - in_u^2 )
1975	else
1976	d_ss = 0
1977	endif
1978
1979	// get d_He
1980	if(abs(in_u) < kTube_ri)
1981	d_He = 2 * sqrt( kTube_ri^2 - in_u^2 )
1982	else
1983	d_He = 0
1984	endif
1985
1986	integrand = T_shexp(-kSig_sslambdad_ss/cos(th_y))( 1-exp(-kSig_8b_Helambdad_He/cos(th_y)) )
1987
1988	return(integrand)
1989	end
1990
1991	//
1992	// Tube efficiency + shadowing
1993	//
1994	// function to generate the table for interpolation
1995	//
1996	// table is generated for a specific wavelength and normalized to eff(lam,0,0)
1997	//
1998	// below 24 deg (theta_x), there is no shadowing, so the table rows are all identical
1999	//
2000	// Only one table is stored, and the wavelength of that table is stored in the wave note
2001	// -- detector correction checks the note, and recalculates the table if needed
2002	// (calculation takes approx 5 seconds)
2003	//
2004	Function V_TubeShadowEfficiencyTables_oneLam(lambda)
2005	Variable lambda
2006
2007	// storage location for tables
2008	SetDataFolder root:Packages:NIST:VSANS:Globals:Efficiency
2009
2010	//make waves that will be filed with the scattering angles and the result of the calculation
2011	//
2012
2013	//// fill arrays with the scattering angles theta_x and theta_y
2014	// 0 < x < 50
2015	// 0 < y < 50
2016
2017	// *** the definitions of x and y for the T/B panels is flipped so that x is always lateral WRT the tubes
2018	// and y is always along the length of the tubes
2019
2020	Variable ii,jj,numx,numy,dx,dy,cos_th,arg,tmp,normVal
2021	numx = 25
2022	numy = numx
2023
2024	Make/O/D/N=(numx,numy) eff
2025	Make/O/D/N=(numx) theta_x, theta_y,eff_with_shadow,lam_cos
2026
2027	SetScale x 0,(numx2)/3602*pi,"", eff
2028	SetScale y 0,(numy2)/3602*pi,"", eff
2029
2030	Note/K eff // clear the note
2031	Note eff "LAMBDA="+num2str(lambda)
2032
2033	// theta_x = p*2
2034	theta_y = p *2 // value range from 0->45, changes if you change numx
2035
2036	//convert degrees to radians
2037	// theta_x = theta_x/3602pi
2038	theta_y = theta_y/3602pi
2039
2040	// Make/O/D/N=12 lam_wave
2041	// lam_wave = {0.5,0.7,1,1.5,2,3,4,6,8,10,15,20}
2042
2043	// Make/O/D/N=(12*numx) eff_withX_to_interp,lam_cos_theta_y
2044	// eff_withX_to_interp=0
2045	// lam_cos_theta_y=0
2046
2047	Make/O/D/N=3 pWave
2048	pWave[0] = lambda
2049
2050
2051	for(ii=0 ;ii<numx;ii+=1)
2052
2053	for(jj=0;jj<numx;jj+=1)
2054
2055	pWave[1] = indexToScale(eff,ii,0) //set theta x
2056	pWave[2] = indexToScale(eff,jj,1) //set theta y
2057
2058	eff_with_shadow[jj] = Integrate1D(V_Efficiency_Integral,-kTube_ri,kTube_ri,2,0,pWave) // adaptive Gaussian quadrature
2059	eff_with_shadow[jj] /= (2*kTube_ri)
2060
2061	eff[ii][jj] = eff_with_shadow[jj]
2062	endfor
2063
2064	//eff[ii][] = eff_with_shadow[q]
2065	endfor
2066
2067	lam_cos = lambda/cos(theta_y)
2068
2069	Sort lam_cos,eff_with_shadow,lam_cos
2070
2071	//
2072	////// // value for normalization at current wavelength
2073	pWave[0] = lambda
2074	pWave[1] = 0
2075	pWave[2] = 0
2076	////
2077	normVal = Integrate1D(V_Efficiency_Integral,-kTube_ri,kTube_ri,2,0,pWave)
2078	normVal /= (2*kTube_ri)
2079	//
2080	// print normVal
2081	//
2082	eff_with_shadow /= normVal // eff(lam,th_x,th_y) / eff(lam,0,0)
2083
2084	eff /= normVal
2085
2086	// TODO
2087	// - clean up after I'm satisfied computations are correct
2088	// KillWaves/Z tmp_theta,tmp_dist,tmp_err,lat_err
2089
2090	SetDataFolder root:
2091	return(0)
2092	end
2093
2094
2095
2096	//
2097	// Tube efficiency + shadowing
2098	//
2099	//
2100	// TESTING function to generate the tables for interpolation
2101	// and various combinations of the corrections for plotting
2102	//
2103	Function V_TubeShadowEfficiencyTables_withX()
2104
2105
2106	Variable lambda
2107	lambda = 6
2108
2109	Variable theta_val=3 //the single theta_x value that is used
2110
2111	// TODO
2112	// -- better storage location for tables
2113	// bad place for now...
2114	SetDataFolder root:
2115
2116	//make waves that will be filed with the scattering angles and the result of the calculation
2117	//
2118
2119	//// fill arrays with the scattering angles theta_x and theta_y
2120	// 0 < x < 50
2121	// 0 < y < 50
2122
2123	// *** the definitions of x and y for the T/B panels is flipped so that x is always lateral WRT the tubes
2124	// and y is always along the length of the tubes
2125
2126	Variable ii,jj,numx,numy,dx,dy,cos_th,arg,tmp,normVal
2127	numx = 10
2128	numy = 10
2129
2130	// Make/O/D/N=(numx,numy) eff
2131	Make/O/D/N=(numx) theta_x, theta_y,eff_with_shadow,lam_cos
2132
2133	theta_x = p*5
2134	theta_y = p*5 // value range from 0->45, changes if you change numx
2135
2136	//convert degrees to radians
2137	theta_x = theta_x/3602pi
2138	theta_y = theta_y/3602pi
2139
2140	Make/O/D/N=12 lam_wave
2141	lam_wave = {0.5,0.7,1,1.5,2,3,4,6,8,10,15,20}
2142
2143	Make/O/D/N=(12*numx) eff_withX_to_interp,lam_cos_theta_y
2144	eff_withX_to_interp=0
2145	lam_cos_theta_y=0
2146
2147	Make/O/D/N=3 pWave
2148
2149	for(jj=0;jj<12;jj+=1)
2150
2151	pWave[0] = lam_wave[jj]
2152
2153	for(ii=0 ;ii<numx;ii+=1)
2154
2155	pWave[1] = theta_val/3602pi //set theta x to any value
2156	pWave[2] = theta_y[ii]
2157
2158	eff_with_shadow[ii] = Integrate1D(V_Efficiency_Integral,-kTube_ri,kTube_ri,2,0,pWave) // adaptive Gaussian quadrature
2159	eff_with_shadow[ii] /= (2*kTube_ri)
2160
2161	endfor
2162
2163	lam_cos = lam_wave[jj]/cos(theta_y)
2164
2165	// messy indexing for the concatentation
2166	lam_cos_theta_y[jjnumx,(jj+1)numx-1] = lam_cos[p-jj*numx]
2167	eff_withX_to_interp[jjnumx,(jj+1)numx-1] = eff_with_shadow[p-jj*numx]
2168
2169	endfor
2170
2171	Sort lam_cos_theta_y,eff_withX_to_interp,lam_cos_theta_y
2172
2173	//
2174	//////// // value for normalization at what wavelength???
2175	// pWave[0] = 6
2176	// pWave[1] = 0
2177	// pWave[2] = 0
2178	//////
2179	// normVal = Integrate1D(V_Efficiency_Integral,-kTube_ri,kTube_ri,2,0,pWave)
2180	// normVal /= (2*kTube_ri)
2181	////
2182	// print normVal
2183	////
2184	// eff_withX_to_interp /= normVal // eff(lam,th_x,th_y) / eff(lam,0,0)
2185
2186	// TODO
2187	// - clean up after I'm satisfied computations are correct
2188	// KillWaves/Z tmp_theta,tmp_dist,tmp_err,lat_err
2189
2190	return(0)
2191	end
2192
2193
2194
2195	/////////////
2196
2197	//
2198	//
2199	// testing function to calculate the correction for the attenuation
2200	// of the scattered beam by windows downstream of the sample
2201	// (the back window of the sample block, the Si window)
2202	//
2203	// For implementation, this function could be made identical
2204	// to the large angle transmission correction, since the math is
2205	// identical - only the Tw value is different (and should ideally be
2206	// quite close to 1). With Tw near 1, this would be a few percent correction
2207	// at the largest scattering angles.
2208	//
2209	//
2210	Function V_WindowTransmission(tw)
2211	Variable tw
2212
2213	Make/O/D/N=100 theta,method1,method2,arg
2214
2215	theta = p/2
2216	theta = theta/3602pi //convert to radians
2217
2218	// method1 = exp( -ln(tw)/cos(theta) )/tw
2219
2220	Variable tau
2221	tau = -ln(tw)
2222	arg = (1-cos(theta))/cos(theta)
2223
2224	if(tau < 0.01)
2225	method2 = 1 - 0.5tauarg
2226	else
2227	method2 = ( 1 - exp(-tauarg) )/(tauarg)
2228	endif
2229
2230	return(0)
2231	end
2232
2233
2234	//
2235	// Large angle transmission correction for the downstream window
2236	//
2237	// DIVIDE the intensity by this correction to get the right answer
2238	//
2239	// -- this is a duplication of the math for the large angle
2240	// sample tranmission correction. Same situation, but now the
2241	// scattered neutrons are attenuated by whatever windows are
2242	// downstream of the scattering event.
2243	// (the back window of the sample block, the Si window)
2244	//
2245	// For implementation, this function is made identical
2246	// to the large angle transmission correction, since the math is
2247	// identical - only the Tw value is different (and should ideally be
2248	// quite close to 1). With Tw near 1, this would be a few percent correction
2249	// at the largest scattering angles.
2250	//
2251	//
2252	Function V_DownstreamWindowTransmission(w,w_err,fname,detStr,destPath)
2253	Wave w,w_err
2254	String fname,detStr,destPath
2255
2256	Variable sdd,xCtr,yCtr,uval
2257
2258	// get all of the geometry information
2259	// orientation = V_getDet_tubeOrientation(fname,detStr)
2260	sdd = V_getDet_ActualDistance(fname,detStr)
2261
2262	// this is ctr in mm
2263	xCtr = V_getDet_beam_center_x_mm(fname,detStr)
2264	yCtr = V_getDet_beam_center_y_mm(fname,detStr)
2265
2266	// get the value of the overall transmission of the downstream components
2267	// + error if available.
2268	// trans = V_getSampleTransmission(fname)
2269	// trans_err = V_getSampleTransError(fname)
2270	// TODO -- HARD WIRED values, need to set a global or find a place in the header (instrument block?) (reduction?)
2271	// currently globals are forced to one in WorkFolderUtils.ipf as the correction is done
2272	NVAR trans = root:Packages:NIST:VSANS:Globals:gDownstreamWinTrans
2273	NVAR trans_err = root:Packages:NIST:VSANS:Globals:gDownstreamWinTransErr
2274
2275	SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr)
2276
2277	Wave data_realDistX = data_realDistX
2278	Wave data_realDistY = data_realDistY
2279
2280	Duplicate/O w dwt_corr,tmp_theta,tmp_dist,dwt_err,tmp_err //in the current df
2281
2282	//// calculate the scattering angle
2283	// dx = (distX - xctr) //delta x in mm
2284	// dy = (distY - yctr) //delta y in mm
2285	tmp_dist = sqrt((data_realDistX - xctr)^2 + (data_realDistY - yctr)^2)
2286
2287	tmp_dist /= 10 // convert mm to cm
2288	// sdd is in [cm]
2289
2290	tmp_theta = atan(tmp_dist/sdd) //this is two_theta, the scattering angle
2291
2292	Variable ii,jj,numx,numy,dx,dy,cos_th,arg,tmp
2293	numx = DimSize(tmp_theta,0)
2294	numy = DimSize(tmp_theta,1)
2295
2296
2297	//optical thickness
2298	uval = -ln(trans) //use natural logarithm
2299
2300	for(ii=0 ;ii<numx;ii+=1)
2301	for(jj=0;jj<numy;jj+=1)
2302
2303	cos_th = cos(tmp_theta[ii][jj])
2304	arg = (1-cos_th)/cos_th
2305
2306	// a Taylor series around uval*arg=0 only needs about 4 terms for very good accuracy
2307	// correction= 1 - 0.5uvalarg + (uvalarg)^2/6 - (uvalarg)^3/24 + (uval*arg)^4/120
2308	// OR
2309	if((uval<0.01) \|\| (cos_th>0.99))
2310	//small arg, approx correction
2311	dwt_corr[ii][jj] = 1-0.5uvalarg
2312	else
2313	//large arg, exact correction
2314	dwt_corr[ii][jj] = (1-exp(-uvalarg))/(uvalarg)
2315	endif
2316
2317	// (DONE)
2318	// x- properly calculate and apply the 2D error propagation
2319	if(trans == 1)
2320	dwt_err[ii][jj] = 0 //no correction, no error
2321	else
2322	//sigT, calculated from the Taylor expansion
2323	tmp = (1/trans)(arg/2-arg^2/3uval+arg^3/8uval^2-arg^4/30uval^3)
2324	tmp *= tmp
2325	tmp *= trans_err^2
2326	tmp = sqrt(tmp) //sigT
2327
2328	dwt_err[ii][jj] = tmp
2329	endif
2330
2331	endfor
2332	endfor
2333
2334	// Here it is! Apply the correction to the intensity (divide -- to get the proper correction)
2335	w /= dwt_corr
2336
2337	// relative errors add in quadrature to the current 2D error
2338	tmp_err = (w_err/dwt_corr)^2 + (dwt_err/dwt_corr)^2ww/dwt_corr^2
2339	tmp_err = sqrt(tmp_err)
2340
2341	w_err = tmp_err
2342
2343	// DONE x- clean up after I'm satisfied computations are correct
2344	KillWaves/Z tmp_theta,tmp_dist,tmp_err,dwt_err
2345
2346	return(0)
2347	end
2348

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format