Context Navigation

V_DetectorCorrections.ipf @ 1095

Last change on this file since 1095 was 1095, checked in by srkline, 4 years ago

Added two model functions for white beam smearing.

Many other small changes for processing of the back detector, shuffling of VSANS menu items, and consistent naming of V_ procedures.

File size: 44.9 KB

Line
1	#pragma rtGlobals=3 // Use modern global access method and strict wave access.
2	#pragma version=1.0
3	#pragma IgorVersion=6.1
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
16
17	//
18	// detector dead time
19	//
20	// input is the data array (N tubes x M pixels)
21	// input of N x 1 array of dead time values
22	//
23	// output is the corrected counts in data, overwriting the input data
24	//
25	// Note that the equation in Roe (eqn 2.15, p. 63) looks different, but it is really the
26	// same old equation, just written in a more complex form.
27	//
28	// (DONE)
29	// x- verify the direction of the tubes and indexing
30	// x- decide on the appropriate functional form for the tubes
31	// x- need count time as input
32	// x- be sure I'm working in the right data folder (all waves are passed in)
33	// x- clean up when done
34	// x- calculate + return the error contribution?
35	// x- verify the error propagation
36	Function V_DeadTimeCorrectionTubes(dataW,data_errW,dtW,ctTime)
37	Wave dataW,data_errW,dtW
38	Variable ctTime
39
40	// do I count on the orientation as an input, or do I just figure it out on my own?
41	String orientation
42	Variable dimX,dimY
43	dimX = DimSize(dataW,0)
44	dimY = DimSize(dataw,1)
45	if(dimX > dimY)
46	orientation = "horizontal"
47	else
48	orientation = "vertical"
49	endif
50
51	// sum the counts in each tube and divide by time for total cr per tube
52	Variable npt
53
54	if(cmpstr(orientation,"vertical")==0)
55	// this is data dimensioned as (Ntubes,Npix)
56
57	MatrixOp/O sumTubes = sumRows(dataW) // n x 1 result
58	sumTubes /= ctTime //now count rate per tube
59
60	dataW[][] = dataW[p][q]/(1-sumTubes[p]*dtW[p]) //correct the data
61	data_errW[][] = data_errW[p][q]/(1-sumTubes[p]*dtW[p]) // propagate the error wave
62
63	elseif(cmpstr(orientation,"horizontal")==0)
64	// this is data (horizontal) dimensioned as (Npix,Ntubes)
65
66	MatrixOp/O sumTubes = sumCols(dataW) // 1 x m result
67	sumTubes /= ctTime
68
69	dataW[][] = dataW[p][q]/(1-sumTubes[q]*dtW[q])
70	data_errW[][] = data_errW[p][q]/(1-sumTubes[q]*dtW[q])
71
72	else
73	DoAlert 0,"Orientation not correctly passed in DeadTimeCorrectionTubes(). No correction done."
74	endif
75
76	return(0)
77	end
78
79	// test function
80	Function V_testDTCor()
81
82	String detStr = ""
83	String fname = "RAW"
84	Variable ctTime
85
86	detStr = "FR"
87	Wave w = V_getDetectorDataW(fname,detStr)
88	Wave w_err = V_getDetectorDataErrW(fname,detStr)
89	Wave w_dt = V_getDetector_deadtime(fname,detStr)
90
91	ctTime = V_getCount_time(fname)
92
93	// ctTime = 10
94	V_DeadTimeCorrectionTubes(w,w_err,w_dt,ctTime)
95
96	End
97
98
99	//
100	// Non-linear data correction
101	//
102	// DOES NOT modify the data, only calculates the spatial relationship
103	//
104	// input is the data array (N tubes x M pixels)
105	// input of N x M array of quadratic coefficients
106	//
107	// output is wave of corrected real space distance corresponding to each pixel of the data
108	// its distance from the nominal beam center of (0,0)
109	//
110	//
111	// (DONE)
112	// x- UNITS!!!! currently this is mm, which certainly doesn't match anything else!!!
113	//
114	// x- verify the direction of the tubes and indexing
115	// x- be sure I'm working in the right data folder (it is passed in, and the full path is used)
116	// x- clean up when done
117	// x- calculate + return the error contribution? (there is none for this operation)
118	// x- do I want this to return a wave? (no, default names are generated)
119	// x- do I need to write a separate function that returns the distance wave for later calculations?
120	// x- do I want to make the distance array 3D to keep the x and y dims together? Calculate them all right now?
121	// x- what else do I need to pass to the function? (fname=folder? detStr?)
122	// y- (yes,see below) need a separate block or function to handle "B" detector which will be ? different
123	//
124	//
125	Function V_NonLinearCorrection(fname,dataW,coefW,tube_width,detStr,destPath)
126	String fname //can also be a folder such as "RAW"
127	Wave dataW,coefW
128	Variable tube_width
129	String detStr,destPath
130
131
132	// do I count on the orientation as an input, or do I just figure it out on my own?
133	String orientation
134	Variable dimX,dimY
135	dimX = DimSize(dataW,0)
136	dimY = DimSize(dataW,1)
137	if(dimX > dimY)
138	orientation = "horizontal"
139	else
140	orientation = "vertical"
141	endif
142
143	// make a wave of the same dimensions, in the same data folder for the distance
144	// ?? or a 3D wave?
145	Make/O/D/N=(dimX,dimY) $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
146	Make/O/D/N=(dimX,dimY) $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
147	Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
148	Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
149
150	// then per tube, do the quadratic calculation to get the real space distance along the tube
151	// the distance perpendicular to the tube is n*(8.4mm) per tube index
152
153	// TODO
154	// -- GAP IS HARD-WIRED as constant values
155	Variable offset,gap
156
157	// kPanelTouchingGap is in mm
158	// the gap is split equally between the panel pairs
159	// (DONE) -- replace all of this with V_getDet_panel_gap(fname,detStr) once it is added to the file
160
161	gap = V_getDet_panel_gap(fname,detStr)
162
163	// TODO:
164	// -- once the gap fields have been verified, this check can be removed
165	// -- it should only apply to data pre-2018 when the field did not exist in the file
166	// -- any VSANS data from 2018+ should read gap from the file.
167
168	if(gap < -100) //-999999 returned if field is missing from file
169
170	if(cmpstr(detStr,"FL") == 0 \|\| cmpstr(detStr,"FR") == 0)
171	gap = 3.5 //mm (measured, JB 1/4/18)
172	endif
173	if(cmpstr(detStr,"FT") == 0 \|\| cmpstr(detStr,"FB") == 0)
174	gap = 3.3 //mm (measured, JB 2/1/18)
175	endif
176	if(cmpstr(detStr,"ML") == 0 \|\| cmpstr(detStr,"MR") == 0)
177	gap = 5.9 //mm (measured, JB 1/4/18)
178	endif
179	if(cmpstr(detStr,"MT") == 0 \|\| cmpstr(detStr,"MB") == 0)
180	gap = 18.3 //mm (measured, JB 2/1/18)
181	endif
182
183	endif
184
185
186	if(cmpstr(orientation,"vertical")==0)
187	// this is data dimensioned as (Ntubes,Npix)
188
189	// adjust the x postion based on the beam center being nominally (0,0) in units of cm, not pixels
190	if(cmpstr(fname,"VCALC")== 0 )
191	offset = VCALC_getPanelTranslation(detStr)
192	offset *= 10 // convert to units of mm
193	// if(cmpstr("L",detStr[1]) == 0)
194	// offset *= -1 //negative value for L
195	// endif
196	else
197	//normal case
198	offset = V_getDet_LateralOffset(fname,detStr)
199	offset *= 10 //convert cm to mm
200	endif
201
202	// calculation is in mm, not cm
203	// offset will be a negative value for the L panel, and positive for the R panel
204	if(kBCTR_CM)
205	if(cmpstr("L",detStr[1]) == 0)
206	// data_realDistX[][] = offset - (dimX - p)*tube_width // TODO should this be dimX-1-p = 47-p?
207	data_realDistX[][] = offset - (dimX - p - 1/2)*tube_width - gap/2 // TODO should this be dimX-1-p = 47-p?
208	else
209	// right
210	// data_realDistX[][] = tube_width*(p+1) + offset + gap //add to the Right det,
211	data_realDistX[][] = tube_width*(p+1/2) + offset + gap/2 //add to the Right det
212	endif
213	else
214	data_realDistX[][] = tube_width*(p)
215	endif
216	data_realDistY[][] = coefW[0][p] + coefW[1][p]q + coefW[2][p]q*q
217
218
219	elseif(cmpstr(orientation,"horizontal")==0)
220	// this is data (horizontal) dimensioned as (Npix,Ntubes)
221	data_realDistY[][] = tube_width*q
222
223	if(cmpstr(fname,"VCALC")== 0 )
224	offset = VCALC_getPanelTranslation(detStr)
225	offset *= 10 // convert to units of mm
226	// if(cmpstr("B",detStr[1]) == 0)
227	// offset *= -1 // negative value for Bottom det
228	// endif
229	else
230	//normal case
231	offset = V_getDet_VerticalOffset(fname,detStr)
232	offset *= 10 //convert cm to mm
233	endif
234
235	if(kBCTR_CM)
236	if(cmpstr("T",detStr[1]) == 0)
237	// data_realDistY[][] = tube_width*(q+1) + offset + gap
238	data_realDistY[][] = tube_width*(q+1/2) + offset + gap/2
239	else
240	// bottom
241	// data_realDistY[][] = offset - (dimY - q)*tube_width // TODO should this be dimY-1-q = 47-q?
242	data_realDistY[][] = offset - (dimY - q - 1/2)*tube_width - gap/2 // TODO should this be dimY-1-q = 47-q?
243	endif
244	else
245	data_realDistY[][] = tube_width*(q)
246	endif
247	data_realDistX[][] = coefW[0][q] + coefW[1][q]p + coefW[2][q]p*p
248
249	else
250	DoAlert 0,"Orientation not correctly passed in NonLinearCorrection(). No correction done."
251	return(0)
252	endif
253
254	return(0)
255	end
256
257
258
259
260	// TODO:
261	// -- the cal_x and y coefficients are totally fake
262	// -- the wave assignment may not be correct.. so beware
263	//
264	//
265	Function V_NonLinearCorrection_B(folder,dataW,cal_x,cal_y,detStr,destPath)
266	String folder
267	Wave dataW,cal_x,cal_y
268	String detStr,destPath
269
270	if(cmpstr(detStr,"B") != 0)
271	return(0)
272	endif
273
274	Print "*Cal_X and Cal_Y for Back are using default values instead of file values *"
275
276	cal_x[0] = VCALC_getPixSizeX(detStr)*10 // pixel size in mm VCALC_getPixSizeX(detStr) is [cm]
277	cal_x[1] = 1
278	cal_x[2] = 10000
279	cal_y[0] = VCALC_getPixSizeY(detStr)*10 // pixel size in mm VCALC_getPixSizeX(detStr) is [cm]
280	cal_y[1] = 1
281	cal_y[2] = 10000
282
283
284	// do I count on the orientation as an input, or do I just figure it out on my own?
285	Variable dimX,dimY
286
287	// Wave dataW = V_getDetectorDataW(folder,detStr)
288
289	dimX = DimSize(dataW,0)
290	dimY = DimSize(dataW,1)
291
292	// make a wave of the same dimensions, in the same data folder for the distance
293	// ?? or a 3D wave?
294	Make/O/D/N=(dimX,dimY) $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
295	Make/O/D/N=(dimX,dimY) $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
296	Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
297	Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
298
299	// Wave cal_x = V_getDet_cal_x(folder,detStr)
300	// Wave cal_y = V_getDet_cal_y(folder,detStr)
301
302	data_realDistX[][] = cal_x[0]*p
303	data_realDistY[][] = cal_y[0]*q
304
305	return(0)
306	end
307
308
309	//
310	//
311	// TODO
312	// -- VERIFY the calculations
313	// -- verify where this needs to be done (if the beam center is changed)
314	// -- then the q-calculation needs to be re-done
315	// -- the position along the tube length is referenced to tube[0], for no particular reason
316	// It may be better to take an average? but [0] is an ASSUMPTION
317	// -- distance along tube is simple interpolation, or do I use the coefficients to
318	// calculate the actual value
319	//
320	// -- distance in the lateral direction is based on tube width, which is a fixed parameter
321	//
322	//
323	Function V_ConvertBeamCtr_to_mm(folder,detStr,destPath)
324	String folder,detStr,destPath
325
326	Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
327	Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
328
329	String orientation
330	Variable dimX,dimY,xCtr,yCtr
331	dimX = DimSize(data_realDistX,0)
332	dimY = DimSize(data_realDistX,1)
333	if(dimX > dimY)
334	orientation = "horizontal"
335	else
336	orientation = "vertical"
337	endif
338
339	xCtr = V_getDet_beam_center_x(folder,detStr)
340	yCtr = V_getDet_beam_center_y(folder,detStr)
341
342	Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
343	Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
344	WAVE x_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
345	WAVE y_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
346
347	Variable tube_width = V_getDet_tubeWidth(folder,detStr)
348
349	//
350	if(cmpstr(orientation,"vertical")==0)
351	// this is data dimensioned as (Ntubes,Npix)
352	// data_realDistX[][] = tube_width*p
353	// data_realDistY[][] = coefW[0][p] + coefW[1][p]q + coefW[2][p]q*q
354	x_mm[0] = tube_width*xCtr
355	y_mm[0] = data_realDistY[0][yCtr]
356	else
357	// this is data (horizontal) dimensioned as (Npix,Ntubes)
358	// data_realDistX[][] = coefW[0][q] + coefW[1][q]p + coefW[2][q]p*p
359	// data_realDistY[][] = tube_width*q
360	x_mm[0] = data_realDistX[xCtr][0]
361	y_mm[0] = tube_width*yCtr
362	endif
363
364	return(0)
365	end
366
367	//
368	//
369	// (DONE)
370	// x- VERIFY the calculations
371	// x- verify where this needs to be done (if the beam center is changed)
372	// x- then the q-calculation needs to be re-done
373	// x- the position along the tube length is referenced to tube[0], for no particular reason
374	// It may be better to take an average? but [0] is an ASSUMPTION
375	// x- distance along tube is simple interpolation
376	//
377	// x- distance in the lateral direction is based on tube width, which is a fixed parameter
378	//
379	// the value in pixels is written to the local data folder, NOT to disk (it is recalculated as needed)
380	//
381	Function V_ConvertBeamCtr_to_pix(folder,detStr,destPath)
382	String folder,detStr,destPath
383
384	Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
385	Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
386
387	String orientation
388	Variable dimX,dimY,xCtr,yCtr
389	dimX = DimSize(data_realDistX,0)
390	dimY = DimSize(data_realDistX,1)
391	if(dimX > dimY)
392	orientation = "horizontal"
393	else
394	orientation = "vertical"
395	endif
396
397	xCtr = V_getDet_beam_center_x(folder,detStr) //these are in cm
398	yCtr = V_getDet_beam_center_y(folder,detStr)
399
400	Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_pix")
401	Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_pix")
402	WAVE x_pix = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_pix")
403	WAVE y_pix = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_pix")
404
405	Variable tube_width = V_getDet_tubeWidth(folder,detStr)
406
407	variable edge,delta
408	Variable gap
409
410	// kPanelTouchingGap is in mm
411	// the gap is split equally between the panel pairs
412	// TODO -- replace all of this with V_getDet_panel_gap(fname,detStr) once it is added to the file
413	// these hard-wired values were determined from 6A and WB beam centers. LR values were exactly the same for
414	// both beam considitions (+/- 0.0 mm). FTB was +/- 0.8 mm, MTB +/- 2 mm
415	if(cmpstr(detStr,"FL") == 0 \|\| cmpstr(detStr,"FR") == 0)
416	gap = 3.8 //mm
417	endif
418	if(cmpstr(detStr,"FT") == 0 \|\| cmpstr(detStr,"FB") == 0)
419	gap = 5 //mm
420	endif
421	if(cmpstr(detStr,"ML") == 0 \|\| cmpstr(detStr,"MR") == 0)
422	gap = 5.9 //mm
423	endif
424	if(cmpstr(detStr,"MT") == 0 \|\| cmpstr(detStr,"MB") == 0)
425	gap = 5 //mm
426	endif
427	// TODO: this is the line to keep, to replace the hard-wired values
428	// gap = V_getDet_panel_gap(fname,detStr)
429
430	//
431	if(cmpstr(orientation,"vertical")==0)
432	// this is data dimensioned as (Ntubes,Npix)
433
434	if(kBCTR_CM)
435	if(cmpstr("L",detStr[1]) == 0)
436	edge = data_realDistX[47][0] //tube 47
437	delta = abs(xCtr*10 - edge)
438	x_pix[0] = dimX-1 + delta/tube_width
439	else
440	// R panel
441	edge = data_realDistX[0][0]
442	delta = abs(xCtr*10 - edge + gap)
443	x_pix[0] = -delta/tube_width //since the left edge of the R panel is pixel 0
444	endif
445	endif
446
447	Make/O/D/N=(dimY) tmpTube
448	tmpTube = data_RealDistY[0][p]
449	FindLevel /P/Q tmpTube, yCtr
450
451	y_pix[0] = V_levelX
452	KillWaves/Z tmpTube
453	else
454	// this is data (horizontal) dimensioned as (Npix,Ntubes)
455
456	if(kBCTR_CM)
457	if(cmpstr("T",detStr[1]) == 0)
458	edge = data_realDistY[0][0] //tube 0
459	delta = abs(yCtr*10 - edge + gap)
460	y_pix[0] = -delta/tube_width //since the bottom edge of the T panel is pixel 0
461	else
462	// FM(B) panel
463	edge = data_realDistY[0][47] //y tube 47
464	delta = abs(yCtr*10 - edge)
465	y_pix[0] = dimY-1 + delta/tube_width //since the top edge of the B panels is pixel 47
466	endif
467	endif
468
469
470	Make/O/D/N=(dimX) tmpTube
471	tmpTube = data_RealDistX[p][0]
472	FindLevel /P/Q tmpTube, xCtr
473
474	x_pix[0] = V_levelX
475	KillWaves/Z tmpTube
476
477
478	endif
479
480	return(0)
481	end
482
483	// converts from [cm] beam center to pixels
484	//
485	// the value in pixels is written to the local data folder, NOT to disk (it is recalculated as needed)
486	//
487	Function V_ConvertBeamCtr_to_pixB(folder,detStr,destPath)
488	String folder,detStr,destPath
489
490	Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
491	Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
492
493	Variable dimX,dimY,xCtr,yCtr
494	dimX = DimSize(data_realDistX,0)
495	dimY = DimSize(data_realDistX,1)
496
497	xCtr = V_getDet_beam_center_x(folder,detStr) //these are in cm, *10 to get mm
498	yCtr = V_getDet_beam_center_y(folder,detStr)
499
500	Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_pix")
501	Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_pix")
502	WAVE x_pix = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_pix")
503	WAVE y_pix = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_pix")
504
505
506	// simple wave lookup
507	// can't use x_pix[0] = data_RealDistX(xCtr)[0] since the data has no x-scale and (xCtr) is interpreted
508	// as a point value
509
510	//
511	//xCtr, yCtr are in cm, *10 to get mm to compare to distance array
512
513	Make/O/D/N=(dimX) tmpTube
514	tmpTube = data_RealDistX[p][0]
515	FindLevel /P/Q tmpTube, xCtr*10
516
517	x_pix[0] = V_levelX
518	KillWaves/Z tmpTube
519
520
521	Make/O/D/N=(dimY) tmpTube
522	tmpTube = data_RealDistY[0][p]
523	FindLevel /P/Q tmpTube, yCtr*10
524
525	y_pix[0] = V_levelX
526	KillWaves/Z tmpTube
527
528	print "pixel ctr B = ",x_pix[0],y_pix[0]
529
530	return(0)
531	end
532
533	//
534	//
535	// TODO
536	// -- VERIFY the calculations
537	// -- verify where this needs to be done (if the beam center is changed)
538	// -- then the q-calculation needs to be re-done
539	//
540	// -- not much is known about the "B" detector, so this
541	// all hinges on the non-linear corrections being done correctly for that detector
542	//
543	// Variable detCtrX, detCtrY
544	// // get the pixel center of the detector (not the beam center)
545	// detCtrX = trunc( DimSize(dataW,0)/2 ) //
546	// detCtrY = trunc( DimSize(dataW,1)/2 )
547	//
548	//
549	Function V_ConvertBeamCtr_to_mmB(folder,detStr,destPath)
550	String folder,detStr,destPath
551
552
553	DoAlert 0,"Error - Beam center is being interpreted as pixels, but needs to be in cm. V_ConvertBeamCtr_to_mmB()"
554
555	Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
556	Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
557
558	Variable xCtr,yCtr
559	xCtr = V_getDet_beam_center_x(folder,detStr)
560	yCtr = V_getDet_beam_center_y(folder,detStr)
561
562	Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
563	Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
564	WAVE x_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
565	WAVE y_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
566
567	x_mm[0] = data_realDistX[xCtr][0]
568	y_mm[0] = data_realDistY[0][yCtr]
569
570	return(0)
571	end
572
573
574
575
576
577
578	/////
579	//
580	// non-linear corrections to the tube pixels
581	// - returns the distance in mm (although this may change)
582	//
583	// c0,c1,c2,pix
584	// c0-c2 are the fit coefficients
585	// pix is the test pixel
586	//
587	// returns the distance in mm (relative to ctr pixel)
588	// ctr is the center pixel, as defined when fitting to quadratic was done
589	//
590	Function V_TubePixel_to_mm(c0,c1,c2,pix)
591	Variable c0,c1,c2,pix
592
593	Variable dist
594	dist = c0 + c1pix + c2pix*pix
595
596	return(dist)
597	End
598	//
599	////
600
601
602	//
603	// TESTING ONLY
604	Proc V_MakeFakeCalibrationWaves()
605	// make these in the RAW data folder, before converting to a work folder
606	// - then they will be "found" by get()
607	// -- only for the tube, not the Back det
608
609	// DoAlert 0, "re-do this and do a better job of filling the fake calibration data"
610
611	DoAlert 0, "Calibration waves are read in from the data file"
612
613	// V_fMakeFakeCalibrationWaves()
614	End
615
616
617
618	//
619	// TESTING ONLY
620	//
621	// orientation does not matter, there are 48 tubes in each bank
622	// so dimension (3,48) for everything.
623	//
624	// -- but the orientation does indicate TB vs LR, which has implications for
625	// the (fictional) dimension of the pixel along the tube axis, at least as far
626	// as for making the fake coefficients.
627	//
628	Function V_fMakeFakeCalibrationWaves()
629
630	Variable ii,pixSize
631	String detStr,fname="RAW",orientation
632
633	for(ii=0;ii<ItemsInList(ksDetectorListNoB);ii+=1)
634	detStr = StringFromList(ii, ksDetectorListNoB, ";")
635	// Wave w = V_getDetectorDataW(fname,detStr)
636	Make/O/D/N=(3,48) $("root:Packages:NIST:VSANS:RAW:entry:instrument:detector_"+detStr+":spatial_calibration")
637	Wave calib = $("root:Packages:NIST:VSANS:RAW:entry:instrument:detector_"+detStr+":spatial_calibration")
638	// !!!! this overwrites what is there
639
640	orientation = V_getDet_tubeOrientation(fname,detStr)
641	if(cmpstr(orientation,"vertical")==0)
642	// this is vertical tube data dimensioned as (Ntubes,Npix)
643	pixSize = 8.4 //V_getDet_y_pixel_size(fname,detStr)
644
645	elseif(cmpstr(orientation,"horizontal")==0)
646	// this is data (horizontal) dimensioned as (Npix,Ntubes)
647	pixSize = 4 //V_getDet_x_pixel_size(fname,detStr)
648
649	else
650	DoAlert 0,"Orientation not correctly passed in NonLinearCorrection(). No correction done."
651	endif
652
653	calib[0][] = -(128/2)pixSize //approx (n/2)pixSixe
654	calib[1][] = pixSize
655	calib[2][] = 2e-4
656
657	endfor
658
659	// now fake calibration for "B"
660	Wave cal_x = V_getDet_cal_x("RAW","B")
661	Wave cal_y = V_getDet_cal_y("RAW","B")
662
663	cal_x = .34 // mm, ignore the other 2 values
664	cal_y = .34 // mm
665	return(0)
666	End
667
668	//
669	// (DONE)
670	// x- MUST VERIFY the definition of SDD and how (if) setback is written to the data files
671	// x- currently I'm assuming that the SDD is the "nominal" value which is correct for the
672	// L/R panels, but is not correct for the T/B panels (must add in the setback)
673	//
674	//
675	//
676	// data_realDistX, Y must be previously generated from running NonLinearCorrection()
677	//
678	// call with:
679	// fname as the WORK folder, "RAW"
680	// detStr = detector String, "FL"
681	// destPath = path to destination WORK folder ("root:Packages:NIST:VSANS:"+folder)
682	//
683	Function V_Detector_CalcQVals(fname,detStr,destPath)
684	String fname,detStr,destPath
685
686	String orientation
687	Variable xCtr,yCtr,lambda,sdd
688
689	// get all of the geometry information
690	orientation = V_getDet_tubeOrientation(fname,detStr)
691
692
693	sdd = V_getDet_ActualDistance(fname,detStr) //sdd derived, including setback [cm]
694
695	// this is the ctr in pixels --xx-- (now it is in cm!)
696	// xCtr = V_getDet_beam_center_x(fname,detStr)
697	// yCtr = V_getDet_beam_center_y(fname,detStr)
698	// this is ctr in mm
699	xCtr = V_getDet_beam_center_x_mm(fname,detStr)
700	yCtr = V_getDet_beam_center_y_mm(fname,detStr)
701	lambda = V_getWavelength(fname)
702	Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
703	Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
704
705	// make the new waves
706	Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qTot_"+detStr)
707	Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qx_"+detStr)
708	Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qy_"+detStr)
709	Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qz_"+detStr)
710	Wave qTot = $(destPath + ":entry:instrument:detector_"+detStr+":qTot_"+detStr)
711	Wave qx = $(destPath + ":entry:instrument:detector_"+detStr+":qx_"+detStr)
712	Wave qy = $(destPath + ":entry:instrument:detector_"+detStr+":qy_"+detStr)
713	Wave qz = $(destPath + ":entry:instrument:detector_"+detStr+":qz_"+detStr)
714
715	// calculate all of the q-values
716	// sdd is passed in [cm]
717
718
719	// after adding in the 680x1656 back detector, load time was 7.8s, without multithreading
720	// with multithreading, 1.9s
721	// qTot = V_CalcQval(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
722	// qx = V_CalcQX(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
723	// qy = V_CalcQY(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
724	// qz = V_CalcQZ(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
725
726	MultiThread qTot = V_CalcQval(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
727	MultiThread qx = V_CalcQX(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
728	MultiThread qy = V_CalcQY(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
729	MultiThread qz = V_CalcQZ(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
730
731	return(0)
732	End
733
734
735	//function to calculate the overall q-value, given all of the necesary trig inputs
736	//
737	// (DONE)
738	// x- verify the calculation (accuracy - in all input conditions)
739	// x- verify the units of everything here, it's currently all jumbled and wrong... and repeated
740	// x- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
741	// to each pixel
742	//
743	//sdd is in [cm]
744	// distX and distY are in [mm]
745	//wavelength is in Angstroms
746	//
747	//returned magnitude of Q is in 1/Angstroms
748	//
749	ThreadSafe Function V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
750	Variable xaxval,yaxval,xctr,yctr,sdd,lam
751	Wave distX,distY
752
753	Variable dx,dy,qval,two_theta,dist
754
755
756	dx = (distX[xaxval][yaxval] - xctr) //delta x in mm
757	dy = (distY[xaxval][yaxval] - yctr) //delta y in mm
758	dist = sqrt(dx^2 + dy^2)
759
760	dist /= 10 // convert mm to cm
761
762	two_theta = atan(dist/sdd)
763
764	qval = 4Pi/lamsin(two_theta/2)
765
766	return qval
767	End
768
769	//calculates just the q-value in the x-direction on the detector
770	// (DONE)
771	// x- verify the calculation (accuracy - in all input conditions)
772	// x- verify the units of everything here, it's currently all jumbled and wrong... and repeated
773	// x- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
774	// to each pixel
775	//
776	//
777	// this properly accounts for qz
778	//
779	ThreadSafe Function V_CalcQX(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
780	Variable xaxval,yaxval,xctr,yctr,sdd,lam
781	Wave distX,distY
782
783	Variable qx,qval,phi,dx,dy,dist,two_theta
784
785	qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
786
787
788	dx = (distX[xaxval][yaxval] - xctr) //delta x in mm
789	dy = (distY[xaxval][yaxval] - yctr) //delta y in mm
790	phi = V_FindPhi(dx,dy)
791
792	//get scattering angle to project onto flat detector => Qr = qval*cos(theta)
793	dist = sqrt(dx^2 + dy^2)
794	dist /= 10 // convert mm to cm
795
796	two_theta = atan(dist/sdd)
797
798	qx = qvalcos(two_theta/2)cos(phi)
799
800	return qx
801	End
802
803	//calculates just the q-value in the y-direction on the detector
804	// (DONE)
805	// x- verify the calculation (accuracy - in all input conditions)
806	// x- verify the units of everything here, it's currently all jumbled and wrong... and repeated
807	// x- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
808	// to each pixel
809	//
810	//
811	// this properly accounts for qz
812	//
813	ThreadSafe Function V_CalcQY(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
814	Variable xaxval,yaxval,xctr,yctr,sdd,lam
815	Wave distX,distY
816
817	Variable qy,qval,phi,dx,dy,dist,two_theta
818
819	qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
820
821
822	dx = (distX[xaxval][yaxval] - xctr) //delta x in mm
823	dy = (distY[xaxval][yaxval] - yctr) //delta y in mm
824	phi = V_FindPhi(dx,dy)
825
826	//get scattering angle to project onto flat detector => Qr = qval*cos(theta)
827	dist = sqrt(dx^2 + dy^2)
828	dist /= 10 // convert mm to cm
829
830	two_theta = atan(dist/sdd)
831
832	qy = qvalcos(two_theta/2)sin(phi)
833
834	return qy
835	End
836
837	//calculates just the q-value in the z-direction on the detector
838	// (DONE)
839	// x- verify the calculation (accuracy - in all input conditions)
840	// x- verify the units of everything here, it's currently all jumbled and wrong... and repeated
841	// x- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
842	// to each pixel
843	//
844	// not actually used for any calculations, but here for completeness if anyone asks, or for 2D data export
845	//
846	// this properly accounts for qz, because it is qz
847	//
848	ThreadSafe Function V_CalcQZ(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
849	Variable xaxval,yaxval,xctr,yctr,sdd,lam
850	Wave distX,distY
851
852	Variable qz,qval,phi,dx,dy,dist,two_theta
853
854	qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
855
856
857	dx = (distX[xaxval][yaxval] - xctr) //delta x in mm
858	dy = (distY[xaxval][yaxval] - yctr) //delta y in mm
859
860	//get scattering angle to project onto flat detector => Qr = qval*cos(theta)
861	dist = sqrt(dx^2 + dy^2)
862	dist /= 10 // convert mm to cm
863
864	two_theta = atan(dist/sdd)
865
866	qz = qval*sin(two_theta/2)
867
868	return qz
869	End
870
871
872	//
873	// (DONE)
874	// x- VERIFY calculations
875	// x- This is the actual solid angle per pixel, not a ratio vs. some "unit SA"
876	// Do I just correct for the different area vs. the "nominal" central area?
877	// x- decide how to implement - YES - directly change the data values (as was done in the past)
878	// 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
879	// would need this to be applied before exporting)
880	// x- do I keep a wave note indicating that this correction has been applied to the data
881	// so that it can be "un-applied"? NO
882	// x- do I calculate theta from geometry directly, or get it from Q (Assuming it's present?)
883	// (YES just from geometry, since I need SDD and dx and dy values...)
884	//
885	//
886	Function V_SolidAngleCorrection(w,w_err,fname,detStr,destPath)
887	Wave w,w_err
888	String fname,detStr,destPath
889
890	Variable sdd,xCtr,yCtr,lambda
891
892	// get all of the geometry information
893	// orientation = V_getDet_tubeOrientation(fname,detStr)
894	sdd = V_getDet_ActualDistance(fname,detStr)
895
896
897	// this is ctr in mm
898	xCtr = V_getDet_beam_center_x_mm(fname,detStr)
899	yCtr = V_getDet_beam_center_y_mm(fname,detStr)
900	lambda = V_getWavelength(fname)
901
902	SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr)
903
904	Wave data_realDistX = data_realDistX
905	Wave data_realDistY = data_realDistY
906
907	Duplicate/O w solid_angle,tmp_theta,tmp_dist //in the current df
908
909	//// calculate the scattering angle
910	// dx = (distX - xctr) //delta x in mm
911	// dy = (distY - yctr) //delta y in mm
912	tmp_dist = sqrt((data_realDistX - xctr)^2 + (data_realDistY - yctr)^2)
913
914	tmp_dist /= 10 // convert mm to cm
915	// sdd is in [cm]
916
917	tmp_theta = atan(tmp_dist/sdd) //this is two_theta, the scattering angle
918
919	Variable ii,jj,numx,numy,dx,dy
920	numx = DimSize(tmp_theta,0)
921	numy = DimSize(tmp_theta,1)
922
923	for(ii=0 ;ii<numx;ii+=1)
924	for(jj=0;jj<numy;jj+=1)
925
926	if(ii==0) //do a forward difference if ii==0
927	dx = (data_realDistX[ii+1][jj] - data_realDistX[ii][jj]) //delta x for the pixel
928	else
929	dx = (data_realDistX[ii][jj] - data_realDistX[ii-1][jj]) //delta x for the pixel
930	endif
931
932
933	if(jj==0)
934	dy = (data_realDistY[ii][jj+1] - data_realDistY[ii][jj]) //delta y for the pixel
935	else
936	dy = (data_realDistY[ii][jj] - data_realDistY[ii][jj-1]) //delta y for the pixel
937	endif
938
939	dx /= 10
940	dy /= 10 // convert mm to cm (since sdd is in cm)
941	solid_angle[ii][jj] = dx*dy //this is in cm^2
942	endfor
943	endfor
944
945	// to cover up any issues w/negative dx or dy
946	solid_angle = abs(solid_angle)
947
948	// solid_angle correction
949	// == dxdycos^3/sdd^2
950	solid_angle *= (cos(tmp_theta))^3
951	solid_angle /= sdd^2
952
953	// Here it is! Apply the correction to the intensity (I divide -- to get the counts per solid angle!!)
954	w /= solid_angle
955
956	//
957	// correctly apply the correction to the error wave (assume a perfect value?)
958	w_err /= solid_angle //
959
960	// DONE x- clean up after I'm satisfied computations are correct
961	KillWaves/Z tmp_theta,tmp_dist
962
963	return(0)
964	end
965
966
967	////////////
968	// TODO: all of below is untested code
969	// copied from SANS
970	//
971	//
972	// NOV 2017
973	// Currently, this is not called from any VSANS routines. it is only referenced
974	// from V_Add_raw_to_work(), which would add two VSANS raw data files together. This has
975	// not yet been implemented. I am only keeping this function around to be sure that
976	// if/when V_Add_raw_to_work() is implemented, all of the functionality of V_DetCorr() is
977	// properly duplicated.
978	//
979	//
980	//
981	//performs solid angle and non-linear detector corrections to raw data as it is "added" to a work folder
982	//function is called by Raw_to_work() and Add_raw_to_work() functions
983	//works on the actual data array, assumes that is is already on LINEAR scale
984	//
985	Function V_DetCorr(data,data_err,realsread,doEfficiency,doTrans)
986	Wave data,data_err,realsread
987	Variable doEfficiency,doTrans
988
989	DoAlert 0,"This has not yet been updated for VSANS"
990
991	Variable xcenter,ycenter,x0,y0,sx,sx3,sy,sy3,xx0,yy0
992	Variable ii,jj,dtdist,dtdis2
993	Variable xi,xd,yd,rad,ratio,domega,xy
994	Variable lambda,trans,trans_err,lat_err,tmp_err,lat_corr
995
996	// Print "...doing jacobian and non-linear corrections"
997
998	NVAR pixelsX = root:myGlobals:gNPixelsX
999	NVAR pixelsY = root:myGlobals:gNPixelsY
1000
1001	//set up values to send to auxiliary trig functions
1002	xcenter = pixelsX/2 + 0.5 // == 64.5 for 128x128 Ordela
1003	ycenter = pixelsY/2 + 0.5 // == 64.5 for 128x128 Ordela
1004
1005	x0 = realsread[16]
1006	y0 = realsread[17]
1007	sx = realsread[10]
1008	sx3 = realsread[11]
1009	sy = realsread[13]
1010	sy3 = realsread[14]
1011
1012	dtdist = 1000*realsread[18] //sdd in mm
1013	dtdis2 = dtdist^2
1014
1015	lambda = realsRead[26]
1016	trans = RealsRead[4]
1017	trans_err = RealsRead[41] //new, March 2011
1018
1019
1020	//waves to contain repeated function calls
1021	Make/O/N=(pixelsX) fyy,xx,yy //Assumes square detector !!!
1022	ii=0
1023	do
1024	xi = ii
1025	// fyy[ii] = dc_fy(ii+1,sy,sy3,ycenter)
1026	// xx[ii] = dc_fxn(ii+1,sx,sx3,xcenter)
1027	// yy[ii] = dc_fym(ii+1,sy,sy3,ycenter)
1028	ii+=1
1029	while(ii<pixelsX)
1030
1031	Make/O/N=(pixelsX,pixelsY) SolidAngle // testing only
1032
1033	ii=0
1034	do
1035	xi = ii
1036	// xd = dc_fx(ii+1,sx,sx3,xcenter)-xx0
1037	jj=0
1038	do
1039	yd = fyy[jj]-yy0
1040	//rad is the distance of pixel ij from the sample
1041	//domega is the ratio of the solid angle of pixel ij versus center pixel
1042	// product xy = 1 for a detector with a linear spatial response (modern Ordela)
1043	// solid angle calculated, dW^3 >=1, so multiply data to raise measured values to correct values.
1044	rad = sqrt(dtdis2 + xd^2 + yd^2)
1045	domega = rad/dtdist
1046	ratio = domega^3
1047	xy = xx[ii]*yy[jj]
1048
1049	data[ii][jj] = xyratio
1050
1051	solidAngle[ii][jj] = xy*ratio //testing only
1052	data_err[ii][jj] = xyratio //error propagation assumes that SA and Jacobian are exact, so simply scale error
1053
1054
1055	// correction factor for detector efficiency JBG memo det_eff_cor2.doc 3/20/07
1056	// correction inserted 11/2007 SRK
1057	// large angle detector efficiency is >= 1 and will "bump up" the measured value at the highest angles
1058	// so divide here to get the correct answer (5/22/08 SRK)
1059	if(doEfficiency)
1060	// data[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)
1061	// data_err[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)
1062	// solidAngle[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd) //testing only
1063	endif
1064
1065	// large angle transmission calculation is <= 1 and will "bump down" the measured value at the highest angles
1066	// so divide here to get the correct answer
1067	if(doTrans)
1068
1069	if(trans<0.1 && ii==0 && jj==0)
1070	Print "*transmission is less than 0.1* and is a significant correction"
1071	endif
1072
1073	if(trans==0)
1074	if(ii==0 && jj==0)
1075	Print "*transmission is ZERO* and has been reset to 1.0 for the averaging calculation"
1076	endif
1077	trans = 1
1078	endif
1079
1080	// pass in the transmission error, and the error in the correction is returned as the last parameter
1081
1082	// lat_corr = V_LargeAngleTransmissionCorr(trans,dtdist,xd,yd,trans_err,lat_err) //moved from 1D avg SRK 11/2007
1083
1084	data[ii][jj] /= lat_corr //divide by the correction factor
1085	//
1086	//
1087	//
1088	// relative errors add in quadrature
1089	tmp_err = (data_err[ii][jj]/lat_corr)^2 + (lat_err/lat_corr)^2data[ii][jj]data[ii][jj]/lat_corr^2
1090	tmp_err = sqrt(tmp_err)
1091
1092	data_err[ii][jj] = tmp_err
1093
1094	// solidAngle[ii][jj] = lat_err
1095
1096
1097	//solidAngle[ii][jj] = LargeAngleTransmissionCorr(trans,dtdist,xd,yd) //testing only
1098	endif
1099
1100	jj+=1
1101	while(jj<pixelsX)
1102	ii+=1
1103	while(ii<pixelsX)
1104
1105	//clean up waves
1106
1107	Return(0)
1108	End
1109
1110
1111	//
1112	// Large angle transmission correction
1113	//
1114	// DIVIDE the intensity by this correction to get the right answer
1115	//
1116	//
1117	// Apply the large angle transmssion correction as the data is converted to WORK
1118	// so that whether the data is saved as 2D or 1D, the correction has properly been done.
1119	//
1120	// This is, however, a SAMPLE dependent calculation, not purely instrument geometry.
1121	//
1122	Function V_LargeAngleTransmissionCorr(w,w_err,fname,detStr,destPath)
1123	Wave w,w_err
1124	String fname,detStr,destPath
1125
1126	Variable sdd,xCtr,yCtr,trans,trans_err,uval
1127
1128	// get all of the geometry information
1129	// orientation = V_getDet_tubeOrientation(fname,detStr)
1130	sdd = V_getDet_ActualDistance(fname,detStr)
1131
1132	// this is ctr in mm
1133	xCtr = V_getDet_beam_center_x_mm(fname,detStr)
1134	yCtr = V_getDet_beam_center_y_mm(fname,detStr)
1135	trans = V_getSampleTransmission(fname)
1136	trans_err = V_getSampleTransError(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 lat_corr,tmp_theta,tmp_dist,lat_err,tmp_err //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 scattering angle
1154
1155	Variable ii,jj,numx,numy,dx,dy,cos_th,arg,tmp
1156	numx = DimSize(tmp_theta,0)
1157	numy = DimSize(tmp_theta,1)
1158
1159
1160	//optical thickness
1161	uval = -ln(trans) //use natural logarithm
1162
1163	for(ii=0 ;ii<numx;ii+=1)
1164	for(jj=0;jj<numy;jj+=1)
1165
1166	cos_th = cos(tmp_theta[ii][jj])
1167	arg = (1-cos_th)/cos_th
1168
1169	// a Taylor series around uval*arg=0 only needs about 4 terms for very good accuracy
1170	// correction= 1 - 0.5uvalarg + (uvalarg)^2/6 - (uvalarg)^3/24 + (uval*arg)^4/120
1171	// OR
1172	if((uval<0.01) \|\| (cos_th>0.99))
1173	//small arg, approx correction
1174	lat_corr[ii][jj] = 1-0.5uvalarg
1175	else
1176	//large arg, exact correction
1177	lat_corr[ii][jj] = (1-exp(-uvalarg))/(uvalarg)
1178	endif
1179
1180	// (DONE)
1181	// x- properly calculate and apply the 2D error propagation
1182	if(trans == 1)
1183	lat_err[ii][jj] = 0 //no correction, no error
1184	else
1185	//sigT, calculated from the Taylor expansion
1186	tmp = (1/trans)(arg/2-arg^2/3uval+arg^3/8uval^2-arg^4/30uval^3)
1187	tmp *= tmp
1188	tmp *= trans_err^2
1189	tmp = sqrt(tmp) //sigT
1190
1191	lat_err[ii][jj] = tmp
1192	endif
1193
1194
1195	endfor
1196	endfor
1197
1198
1199
1200	// Here it is! Apply the correction to the intensity (divide -- to get the proper correction)
1201	w /= lat_corr
1202
1203	// relative errors add in quadrature to the current 2D error
1204	tmp_err = (w_err/lat_corr)^2 + (lat_err/lat_corr)^2ww/lat_corr^2
1205	tmp_err = sqrt(tmp_err)
1206
1207	w_err = tmp_err
1208
1209
1210	// DONE x- clean up after I'm satisfied computations are correct
1211	KillWaves/Z tmp_theta,tmp_dist,tmp_err,lat_err
1212
1213	return(0)
1214	end
1215
1216
1217
1218	//
1219	//test procedure, not called anymore
1220	Proc V_AbsoluteScaling(type,c0,c1,c2,c3,c4,c5,I_err)
1221	String type
1222	Variable c0=1,c1=0.1,c2=0.95,c3=0.1,c4=1,c5=32.0,I_err=0.32
1223	Prompt type,"WORK data type",popup,"CAL;COR;SAM"
1224	Prompt c0, "Sample Transmission"
1225	Prompt c1, "Sample Thickness (cm)"
1226	Prompt c2, "Standard Transmission"
1227	Prompt c3, "Standard Thickness (cm)"
1228	Prompt c4, "I(0) from standard fit (normalized to 1E8 monitor cts)"
1229	Prompt c5, "Standard Cross-Section (cm-1)"
1230	Prompt I_err, "error in I(q=0) (one std dev)"
1231
1232	Variable err
1233	//call the function to do the math
1234	//data from "type" will be scaled and deposited in ABS
1235	err = V_Absolute_Scale(type,c0,c1,c2,c3,c4,c5,I_err)
1236
1237	if(err)
1238	Abort "Error in V_Absolute_Scale()"
1239	endif
1240
1241	//contents are always dumped to ABS
1242	type = "ABS"
1243
1244	//need to update the display with "data" from the correct dataFolder
1245	//reset the current display type to "type"
1246	SVAR gCurDispType = root:Packages:NIST:VSANS:Globals:gCurDispType
1247	gCurDispType = Type
1248
1249	V_UpdateDisplayInformation(Type)
1250
1251	End
1252
1253	//
1254	//
1255	// kappa comes in as s_izero, so be sure to use 1/kappa_err
1256	//
1257	//convert the "type" data to absolute scale using the given standard information
1258	//s_ is the standard
1259	//w_ is the "work" file
1260	//both are work files and should already be normalized to 10^8 monitor counts
1261	Function V_Absolute_Scale(type,w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err)
1262	String type
1263	Variable w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err
1264
1265
1266	Variable defmon = 1e8,w_moncount,s1,s2,s3,s4
1267	Variable scale,trans_err
1268	Variable err,ii
1269	String detStr
1270
1271	// be sure that the starting data exists
1272	err = V_WorkDataExists(type)
1273	if(err==1)
1274	return(err)
1275	endif
1276
1277	//copy from current dir (type) to ABS
1278	V_CopyHDFToWorkFolder(type,"ABS")
1279
1280	// TODO: -- which monitor to use? Here, I think it should already be normalized to 10^8
1281	//
1282	// w_moncount = V_getMonitorCount(type) //monitor count in "type"
1283
1284	w_moncount = V_getBeamMonNormData(type)
1285
1286
1287	if(w_moncount == 0)
1288	//zero monitor counts will give divide by zero ---
1289	DoAlert 0,"Total monitor count in data file is zero. No rescaling of data"
1290	Return(1) //report error
1291	Endif
1292
1293	//calculate scale factor
1294	s1 = defmon/w_moncount // monitor count (s1 should be 1)
1295	s2 = s_thick/w_thick
1296	s3 = s_trans/w_trans
1297	s4 = s_cross/s_izero
1298	scale = s1s2s3*s4
1299
1300	trans_err = V_getSampleTransError(type)
1301
1302	// kappa comes in as s_izero, so be sure to use 1/kappa_err
1303
1304	// and now loop through all of the detectors
1305	//do the actual absolute scaling here, modifying the data in ABS
1306	for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
1307	detStr = StringFromList(ii, ksDetectorListAll, ";")
1308	Wave data = V_getDetectorDataW("ABS",detStr)
1309	Wave data_err = V_getDetectorDataErrW("ABS",detStr)
1310
1311	data *= scale
1312	data_err = sqrt(scale^2data_err^2 + scale^2data^2*(kappa_err^2/s_izero^2 +trans_err^2/w_trans^2))
1313	endfor
1314
1315	//********* 15APR02
1316	// DO NOT correct for atenuators here - the COR step already does this, putting all of the data on equal
1317	// footing (zero atten) before doing the subtraction.
1318
1319	Return (0) //no error
1320	End
1321
1322
1323	//
1324	// TODO:
1325	// -- DoAlert 0,"This has not yet been updated for VSANS"
1326	//
1327	//
1328	// match the attenuation of the RAW data to the "type" data
1329	// so that they can be properly added
1330	//
1331	// are the attenuator numbers the same? if so exit
1332	//
1333	// if not, find the attenuator number for type
1334	// - find both attenuation factors
1335	//
1336	// rescale the raw data to match the ratio of the two attenuation factors
1337	// -- adjust the detector count (rw)
1338	// -- the linear data
1339	//
1340	//
1341	Function V_Adjust_RAW_Attenuation(type)
1342	String type
1343
1344	DoAlert 0,"This has not yet been updated for VSANS"
1345
1346	WAVE rw=$("root:Packages:NIST:RAW:realsread")
1347	WAVE linear_data=$("root:Packages:NIST:RAW:linear_data")
1348	WAVE data=$("root:Packages:NIST:RAW:data")
1349	WAVE data_err=$("root:Packages:NIST:RAW:linear_data_error")
1350	WAVE/T tw = $("root:Packages:NIST:RAW:textRead")
1351
1352	WAVE dest_reals=$("root:Packages:NIST:"+type+":realsread")
1353
1354	Variable dest_atten,raw_atten,tol
1355	Variable lambda,raw_atten_err,raw_AttenFactor,dest_attenFactor,dest_atten_err
1356	String fileStr
1357
1358	dest_atten = dest_reals[3]
1359	raw_atten = rw[3]
1360
1361	tol = 0.1 // within 0.1 atten units is OK
1362	if(abs(dest_atten - raw_atten) < tol )
1363	return(0)
1364	endif
1365
1366	fileStr = tw[3]
1367	lambda = rw[26]
1368	// TODO access correct values
1369	raw_AttenFactor = 1//AttenuationFactor(fileStr,lambda,raw_atten,raw_atten_err)
1370	dest_AttenFactor = 1//AttenuationFactor(fileStr,lambda,dest_atten,dest_atten_err)
1371
1372	rw[2] *= dest_AttenFactor/raw_AttenFactor
1373	linear_data *= dest_AttenFactor/raw_AttenFactor
1374
1375	// to keep "data" and linear_data in sync
1376	data = linear_data
1377
1378	return(0)
1379	End
1380
1381	//
1382	// testing procedure, called from a menu selection
1383	//
1384	Proc V_DIV_a_Workfile(type)
1385	String type
1386	Prompt type,"WORK data type",popup,"SAM;EMP;BGD;ADJ;"
1387
1388	//macro will take whatever is in SELECTED folder and DIVide it by the current
1389	//contents of the DIV folder - the function will check for existence
1390	//before proceeding
1391
1392	Abort "This has not yet been updated for VSANS"
1393
1394	Variable err
1395	err = V_DIVCorrection(type) //returns err = 1 if data doesn't exist in specified folders
1396
1397	if(err)
1398	Abort "error in V_DIVCorrection()"
1399	endif
1400
1401	//contents are NOT always dumped to CAL, but are in the new type folder
1402
1403	String newTitle = "WORK_"+type
1404	DoWindow/F VSANS_Data
1405	DoWindow/T VSANS_Data, newTitle
1406	KillStrings/Z newTitle
1407
1408	//need to update the display with "data" from the correct dataFolder
1409	//reset the current displaytype to "type"
1410	String/G root:Packages:NIST:VSANS:Globals:gCurDispType=Type
1411
1412	V_UpdateDisplayInformation(type)
1413
1414	End
1415
1416
1417	//
1418	// TODO:
1419	// x- DoAlert 0,"This has not yet been updated for VSANS"
1420	// -- how is the error propagation handled? Be sure it is calculated correctly when DIV is generated
1421	// and is applied correctly here...
1422	//
1423	//function will divide the contents of "workType" folder with the contents of
1424	//the DIV folder + detStr
1425	// all data is linear scale for the calculation
1426	//
1427	Function V_DIVCorrection(data,data_err,detStr,workType)
1428	Wave data,data_err
1429	String detStr,workType
1430
1431	//check for existence of data in type and DIV
1432	// if the desired data doesn't exist, let the user know, and abort
1433	String destPath=""
1434
1435	if(WaveExists(data) == 0)
1436	Print "The data wave does not exist in V_DIVCorrection()"
1437	Return(1) //error condition
1438	Endif
1439
1440	//check for DIV
1441	// if the DIV workfile doesn't exist, let the user know,and abort
1442	// !! be sure to check first, before trying to access the wave
1443
1444	// WAVE/Z div_data = $("root:Packages:NIST:VSANS:DIV:entry:instrument:detector_"+detStr+":data")
1445	if(WaveExists($("root:Packages:NIST:VSANS:DIV:entry:instrument:detector_"+detStr+":data")) == 0)
1446	Print "The DIV wave does not exist in V_DIVCorrection()"
1447	Return(1) //error condition
1448	Endif
1449	if(WaveExists($("root:Packages:NIST:VSANS:DIV:entry:instrument:detector_"+detStr+":linear_data_error")) == 0)
1450	Print "The DIV error wave does not exist in V_DIVCorrection()"
1451	Return(1) //error condition
1452	Endif
1453	//files exist, proceed
1454
1455	WAVE/Z div_data_err = V_getDetectorDataErrW("DIV",detStr)
1456	WAVE/Z div_data = V_getDetectorDataW("DIV",detStr)
1457
1458
1459
1460	// do the error propagation first, since data is changed by the correction
1461	data_err = sqrt(data_err^2/div_data^2 + div_data_err^2 * data^2/div_data^4 )
1462
1463	// then the correction
1464	data /= div_data
1465
1466
1467	Return(0)
1468	End
1469
1470
1471	//////////////////////////

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

Download in other formats:

Original Format