Context Navigation

V_DetectorCorrections.ipf @ 1073

Last change on this file since 1073 was 1073, checked in by srkline, 5 years ago

lots of changes here:
many little fixes to clean up TODO items and marke them DONE

changed the handling of the panel "gap" to split the gap evenly. Q-calculations have been re-verified with this change.

re-named the list of "bin Type" values, and added a few more choices. Streamlined how the averaging and plotting works with this list so that it can be more easily modified as different combinations of binning are envisioned. This resulted in a lot of excess code being cut out and replaced with cleaner logic. This change has also been verified to work as intended.

Attenuation is now always calculated from the table. The table also by (NEW) definition has values for the white beam (one waelength) and graphite (multiple possible wavelengths) where the wavelengths are artificially scaled (*1000) or *1e6) so that the interpolations can be done internally without the need for multiple attenuator tables.

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

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

Download in other formats:

Original Format