Context Navigation

VC_DetectorBinning_Utils.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: 47.7 KB

Line
1	#pragma rtGlobals=3 // Use modern global access method and strict wave access.
2
3	/////////////////////////
4	//
5	// Utility functions to:
6	// calculate Q, Qx, Qy, Qz
7	// fill the detector panels with simulated data (the model functions are here)
8	// bin the 2D detector to 1D I(Q) based on Q and deltaQ (bin width)
9	//
10	/////////////////////////
11
12
13
14
15	// x- hard wired for a sphere - change this to allow minimal selections and altering of coefficients
16	// x- add the "fake" 2D simulation to fill the panels which are then later averaged as I(Q)
17	//
18	// NOTE - this is a VCALC only routine, so it's not been made completely generic
19	//
20	Function FillPanel_wModelData(det,qTot,type)
21	Wave det,qTot
22	String type
23
24	// SetDataFolder root:Packages:NIST:VSANS:VCALC:Front
25
26	// q-values and detector arrays already allocated and calculated
27	Duplicate/O det tmpInten,tmpSig,prob_i
28
29	Variable imon,trans,thick,sdd,pixSizeX,pixSizeY,sdd_offset
30
31	//imon = VC_BeamIntensity()*CountTime
32	imon = VCALC_getImon() //TODO: currently from the panel, not calculated
33	trans = 0.8
34	thick = 0.1
35
36	// need SDD
37	// need pixel dimensions
38	// nominal sdd in cm, offset in cm, want result in cm !
39	sdd = VCALC_getSDD(type) + VCALC_getTopBottomSDDSetback(type) // result is sdd in [cm]
40
41	pixSizeX = VCALC_getPixSizeX(type) // cm
42	pixSizeY = VCALC_getPixSizeY(type)
43
44
45	//?? pick the function from a popup on the panel? (bypass the analysis panel, or maybe it's better to
46	// keep the panel to keep people used to using it.)
47	// peak @ 0.1 ~ AgBeh
48	// Make/O/D coef_BroadPeak = {1e-9, 3, 20, 100.0, 0.1,3,0.1}
49	//
50	// peak @ 0.015 in middle of middle detector, maybe not "real" vycor, but that is to be resolved
51	// Make/O/D coef_BroadPeak = {1e-9, 3, 20, 500.0, 0.015,3,0.1}
52	//
53	//
54	Variable addEmpBgd=0
55
56
57	String funcStr = VCALC_getModelFunctionStr()
58	strswitch(funcStr)
59	case "Big Debye":
60	tmpInten = VC_Debye(100,3000,0.0001,qTot[p][q])
61	break
62	case "Big Sphere":
63	tmpInten = VC_SphereForm(1,900,1e-6,0.01,qTot[p][q])
64	break
65	case "Debye":
66	tmpInten = VC_Debye(10,300,0.0001,qTot[p][q])
67	break
68	case "Sphere":
69	tmpInten = VC_SphereForm(1,60,1e-6,0.001,qTot[p][q])
70	break
71	case "AgBeh":
72	tmpInten = VC_BroadPeak(1e-11,3,20,100.0,0.1,3,0.1,qTot[p][q])
73	break
74	case "Vycor":
75	tmpInten = VC_BroadPeak(1e-9,3,20,500.0,0.015,3,0.1,qTot[p][q])
76	break
77	case "Empty Cell":
78	tmpInten = VC_EC_Empirical(2.2e-12,3.346,0.0065,9.0,0.016,qTot[p][q])
79	break
80	case "Blocked Beam":
81	tmpInten = VC_BlockedBeam(0.01,qTot[p][q])
82	break
83	case "Debye +":
84	tmpInten = VC_Debye(10,300,0.0001,qTot[p][q])
85	addEmpBgd = 1
86	break
87	case "AgBeh +":
88	tmpInten = VC_BroadPeak(1e-11,3,20,100.0,0.1,3,0.1,qTot[p][q])
89	addEmpBgd = 1
90	break
91	case "Empty Cell +":
92	tmpInten = VC_EC_Empirical(2.2e-12,3.346,0.0065,9.0,0.016,qTot[p][q])
93	tmpInten += VC_BlockedBeam(0.01,qTot[p][q])
94	break
95	default:
96	tmpInten = VC_Debye(10,300,0.1,qTot[p][q])
97	endswitch
98
99
100	if(addEmpBgd == 1)
101	tmpInten += VC_EC_Empirical(2.2e-12,3.346,0.0065,9.0,0.016,qTot[p][q])
102	tmpInten += VC_BlockedBeam(0.01,qTot[p][q])
103	endif
104
105
106	// x- this is faked to get around the singularity at the center of the back detector
107	//
108	//
109	if(cmpstr(type,"B") == 0)
110	Variable nx,ny,px,py
111	nx = VCALC_get_nPix_X(type)
112	ny = VCALC_get_nPix_Y(type)
113	px = trunc(nx/2)
114	py = trunc(ny/2)
115
116	tmpInten[px][py] = (tmpInten[px][py+1] + tmpInten[px][py-1])/2
117	endif
118
119
120
121	///////////////
122	// // calculate the scattering cross section simply to be able to estimate the transmission
123	// Variable sig_sas=0
124	//
125	// // remember that the random deviate is the coherent portion ONLY - the incoherent background is
126	// // subtracted before the calculation.
127	// CalculateRandomDeviate(funcUnsmeared,$coefStr,wavelength,"root:Packages:NIST:SAS:ran_dev",sig_sas)
128	//
129	// if(sig_sas > 100)
130	// DoAlert 0,"SAS cross section > 100. Estimates of multiple scattering are unreliable. Choosing a model with a well-defined Rg may help"
131	// endif
132	//
133	// // calculate the multiple scattering fraction for display (10/2009)
134	// Variable ii,nMax=10,tau
135	// mScat=0
136	// tau = thick*sig_sas
137	// // this sums the normalized scattering P', so the result is the fraction of multiply coherently scattered
138	// // neutrons out of those that were scattered
139	// for(ii=2;ii<nMax;ii+=1)
140	// mScat += tau^(ii)/factorial(ii)
141	//// print tau^(ii)/factorial(ii)
142	// endfor
143	// estTrans = exp(-1thicksig_sas) //thickness and sigma both in units of cm
144	// mscat *= (estTrans)/(1-estTrans)
145	//
146	//// if(mScat > 0.1) // Display warning
147	//
148	// Print "Sig_sas = ",sig_sas
149	////////////////////
150
151	prob_i = transthickpixSizeXpixSizeY/(sdd)^2tmpInten //probability of a neutron in q-bin(i)
152
153	tmpInten = (imon)*prob_i //tmpInten is not the model calculation anymore!!
154
155
156	/// **** can I safely assume a Gaussian error in the count rate??
157	tmpSig = sqrt(tmpInten) // corrected based on John's memo, from 8/9/99
158
159	tmpInten += gnoise(tmpSig)
160	tmpInten = (tmpInten[p][q] < 0) ? 0 : tmpInten[p][q] // MAR 2013 -- is this the right thing to do
161	tmpInten = trunc(tmpInten)
162
163
164	det = tmpInten
165
166	// if I want "absolute" scale -- then I lose the integer nature of the detector (but keep the random)
167	// det /= transthickpixSizeXpixSizeY/(sdd)^2imon
168
169
170	KillWaves/Z tmpInten,tmpSig,prob_i
171	SetDataFolder root:
172
173	return(0)
174	End
175
176
177	// For a given detector panel, calculate the q-values
178	// -work with everything as arrays
179	// Input needed:
180	// detector data
181	// detector type (LRTB?)
182	// beam center (may be off the detector)
183	// SDD
184	// lambda
185	//
186	// pixel dimensions for detector type (global constants)
187	// - data dimensions read directly from array
188	//
189	// --What is calculated:
190	// array of Q
191	// array of qx,qy,qz
192	// array of error already exists
193	//
194	//
195	// -- sdd in cm
196	// -- lambda in Angstroms
197	Function VC_Detector_2Q(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY)
198	Wave data,qTot,qx,qy,qz
199	Variable xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY
200
201	// loop over the array and calculate the values - this is done as a wave assignment
202	// TODO -- be sure that it's p,q -- or maybe p+1,q+1 as used in WriteQIS.ipf
203	qTot = VC_CalcQval(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
204	qx = VC_CalcQX(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
205	qy = VC_CalcQY(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
206	qz = VC_CalcQZ(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
207
208	return(0)
209	End
210
211
212	// for testing, a version that will calculate the q-arrays for VCALC based on whatever nonlinear coefficients
213	// exist in the RAW data folder
214	//
215	// reverts to the "regular" linear detector if waves not found or a flag is set
216	//
217	// need to call the VSANS V_CalcQval routines (these use the real-space distance, not pixel dims)
218	//
219	// ***** everything passed in is [cm], except for wavelength [A]
220	//
221	// **** TODO :: calibration constants are still in [mm]
222	//
223	//
224	// TODO:
225	// -- tube width is hard-wired in
226	//
227	//
228	Function VC_Detector_2Q_NonLin(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY,detStr)
229	Wave data,qTot,qx,qy,qz
230	Variable xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY
231	String detStr
232
233	String destPath = "root:Packages:NIST:VSANS:VCALC"
234
235	// be sure that the real distance waves exist
236	// TODO -- this may not be the best location?
237
238	// calibration waves do not exist yet, so make some fake ones '
239	// do I count on the orientation as an input, or do I just figure it out on my own?
240	String orientation
241	Variable dimX,dimY
242	dimX = DimSize(data,0)
243	dimY = DimSize(data,1)
244	if(dimX > dimY)
245	orientation = "horizontal"
246	else
247	orientation = "vertical"
248	endif
249
250	if(cmpstr(orientation,"vertical")==0)
251	Make/O/D/N=(3,48) tmpCalib
252	// for the "tall" L/R banks
253	tmpCalib[0][] = -512
254	tmpCalib[1][] = 8
255	tmpCalib[2][] = 0
256	else
257	Make/O/D/N=(3,48) tmpCalib
258	// for the "short" T/B banks
259	tmpCalib[0][] = -256
260	tmpCalib[1][] = 4
261	tmpCalib[2][] = 0
262	endif
263	// override if back panel
264	if(cmpstr(detStr,"B") == 0)
265	// and for the back detector "B"
266	Make/O/D/N=3 tmpCalibX,tmpCalibY
267	tmpCalibX[0] = VCALC_getPixSizeX(detStr)*10 // pixel size in mm VCALC_getPixSizeX(detStr) is [cm]
268	tmpCalibX[1] = 1
269	tmpcalibX[2] = 10000
270	tmpCalibY[0] = VCALC_getPixSizeY(detStr)*10 // pixel size in mm VCALC_getPixSizeX(detStr) is [cm]
271	tmpCalibY[1] = 1
272	tmpcalibY[2] = 10000
273	endif
274
275	// Wave w_calib = V_getDetTube_spatialCalib("VCALC",detStr)
276	Variable tube_width = 8.4 // TODO: UNITS!!! Hard-wired value in [mm]
277	if(cmpstr(detStr,"B") == 0)
278	V_NonLinearCorrection_B("VCALC",data,tmpCalibX,tmpCalibY,detStr,destPath)
279	else
280	V_NonLinearCorrection("VCALC",data,tmpCalib,tube_width,detStr,destPath)
281	endif
282
283	Wave/Z data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
284	Wave/Z data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
285	NVAR gUseNonLinearDet = root:Packages:NIST:VSANS:VCALC:gUseNonLinearDet
286
287	if(kBCTR_CM)
288	if(gUseNonLinearDet && WaveExists(data_realDistX) && WaveExists(data_realDistY))
289	// beam ctr is in cm already
290
291	// calculate all of the q-values
292	qTot = V_CalcQval(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
293	qx = V_CalcQX(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
294	qy = V_CalcQY(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
295	qz = V_CalcQZ(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
296
297	// Print "det, x_mm, y_mm ",detStr,num2str(newX),num2str(newY)
298	// Print "det, x_pix, y_pix ",detStr,num2str(xCtr),num2str(yCtr)
299	else
300	// do the q-calculation using linear detector
301	//VC_Detector_2Q(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY)
302	qTot = V_CalcQval(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
303	qx = V_CalcQX(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
304	qy = V_CalcQY(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
305	qz = V_CalcQZ(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
306	endif
307
308
309	else
310	// using the old calculation with beam center in pixels
311	if(gUseNonLinearDet && WaveExists(data_realDistX) && WaveExists(data_realDistY))
312	// convert the beam centers to mm
313	// String orientation
314	Variable newX,newY
315	dimX = DimSize(data_realDistX,0)
316	dimY = DimSize(data_realDistX,1)
317	if(dimX > dimY)
318	orientation = "horizontal"
319	else
320	orientation = "vertical"
321	endif
322
323
324	//
325	if(cmpstr(orientation,"vertical")==0)
326	// this is data dimensioned as (Ntubes,Npix)
327	newX = tube_width*xCtr
328	newY = data_realDistY[0][yCtr]
329	else
330	// this is data (horizontal) dimensioned as (Npix,Ntubes)
331	newX = data_realDistX[xCtr][0]
332	newY = tube_width*yCtr
333	endif
334
335	//if detector "B", different calculation for the centers (not tubes)
336	if(cmpstr(detStr,"B")==0)
337	newX = data_realDistX[xCtr][0]
338	newY = data_realDistY[0][yCtr]
339	//newX = xCtr
340	//newY = yCtr
341	endif
342
343	// calculate all of the q-values
344	qTot = V_CalcQval(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
345	qx = V_CalcQX(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
346	qy = V_CalcQY(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
347	qz = V_CalcQZ(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
348
349	// Print "det, x_mm, y_mm ",detStr,num2str(newX),num2str(newY)
350	// Print "det, x_pix, y_pix ",detStr,num2str(xCtr),num2str(yCtr)
351	else
352	// do the q-calculation using linear detector
353	VC_Detector_2Q(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY)
354	endif
355
356	endif
357
358	KillWaves/Z tmpCalib,tmpCalibX,tmpCalibY
359
360	return(0)
361	End
362
363
364	//////////////////////
365	// NOTE: The Q calculations are different than what is in GaussUtils in that they take into
366	// accout the different x/y pixel sizes and the beam center not being on the detector -
367	// off a different edge for each LRTB type
368	/////////////////////
369
370	//function to calculate the overall q-value, given all of the necesary trig inputs
371	//and are in detector coordinates (1,128) rather than axis values
372	//the pixel locations need not be integers, reals are ok inputs
373	//sdd is in [cm]
374	//wavelength is in Angstroms
375	//
376	//returned magnitude of Q is in 1/Angstroms
377	//
378	//
379	Function VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
380	Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
381
382	Variable dx,dy,qval,two_theta,dist
383
384
385	dx = (xaxval - xctr)*pixSizeX //delta x in cm
386	dy = (yaxval - yctr)*pixSizeY //delta y in cm
387	dist = sqrt(dx^2 + dy^2)
388
389	two_theta = atan(dist/sdd)
390
391	qval = 4Pi/lamsin(two_theta/2)
392
393	return qval
394	End
395
396	//calculates just the q-value in the x-direction on the detector
397	//input/output is the same as CalcQval()
398	//ALL inputs are in detector coordinates
399	//
400	//sdd is in [cm]
401	//wavelength is in Angstroms
402	//
403	// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst)
404	// now properly accounts for qz
405	//
406	Function VC_CalcQX(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
407	Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
408
409	Variable qx,qval,phi,dx,dy,dist,two_theta
410
411	qval = VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
412
413	// sdd *=100 //convert to cm
414	dx = (xaxval - xctr)*pixSizeX //delta x in cm
415	dy = (yaxval - yctr)*pixSizeY //delta y in cm
416	phi = V_FindPhi(dx,dy)
417
418	//get scattering angle to project onto flat detector => Qr = qval*cos(theta)
419	dist = sqrt(dx^2 + dy^2)
420	two_theta = atan(dist/sdd)
421
422	qx = qvalcos(two_theta/2)cos(phi)
423
424	return qx
425	End
426
427	//calculates just the q-value in the y-direction on the detector
428	//input/output is the same as CalcQval()
429	//ALL inputs are in detector coordinates
430	//sdd is in [cm]
431	//wavelength is in Angstroms
432	//
433	// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst)
434	// now properly accounts for qz
435	//
436	Function VC_CalcQY(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
437	Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
438
439	Variable dy,qval,dx,phi,qy,dist,two_theta
440
441	qval = VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
442
443	// sdd *=100 //convert to cm
444	dx = (xaxval - xctr)*pixSizeX //delta x in cm
445	dy = (yaxval - yctr)*pixSizeY //delta y in cm
446	phi = V_FindPhi(dx,dy)
447
448	//get scattering angle to project onto flat detector => Qr = qval*cos(theta)
449	dist = sqrt(dx^2 + dy^2)
450	two_theta = atan(dist/sdd)
451
452	qy = qvalcos(two_theta/2)sin(phi)
453
454	return qy
455	End
456
457	//calculates just the z-component of the q-vector, not measured on the detector
458	//input/output is the same as CalcQval()
459	//ALL inputs are in detector coordinates
460	//sdd is in [cm]
461	//wavelength is in Angstroms
462	//
463	// not actually used, but here for completeness if anyone asks
464	//
465	Function VC_CalcQZ(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
466	Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
467
468	Variable dy,qval,dx,phi,qz,dist,two_theta
469
470	qval = VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
471
472	// sdd *=100 //convert to cm
473
474	//get scattering angle to project onto flat detector => Qr = qval*cos(theta)
475	dx = (xaxval - xctr)*pixSizeX //delta x in cm
476	dy = (yaxval - yctr)*pixSizeY //delta y in cm
477	dist = sqrt(dx^2 + dy^2)
478	two_theta = atan(dist/sdd)
479
480	qz = qval*sin(two_theta/2)
481
482	return qz
483	End
484
485	//phi is defined from +x axis, proceeding CCW around [0,2Pi]
486	Threadsafe Function V_FindPhi(vx,vy)
487	variable vx,vy
488
489	variable phi
490
491	phi = atan(vy/vx) //returns a value from -pi/2 to pi/2
492
493	// special cases
494	if(vx==0 && vy > 0)
495	return(pi/2)
496	endif
497	if(vx==0 && vy < 0)
498	return(3*pi/2)
499	endif
500	if(vx >= 0 && vy == 0)
501	return(0)
502	endif
503	if(vx < 0 && vy == 0)
504	return(pi)
505	endif
506
507
508	if(vx > 0 && vy > 0)
509	return(phi)
510	endif
511	if(vx < 0 && vy > 0)
512	return(phi + pi)
513	endif
514	if(vx < 0 && vy < 0)
515	return(phi + pi)
516	endif
517	if( vx > 0 && vy < 0)
518	return(phi + 2*pi)
519	endif
520
521	return(phi)
522	end
523
524	Function VC_SphereForm(scale,radius,delrho,bkg,x)
525	Variable scale,radius,delrho,bkg
526	Variable x
527
528	// variables are:
529	//[0] scale
530	//[1] radius (A)
531	//[2] delrho (A-2)
532	//[3] background (cm-1)
533
534	// Variable scale,radius,delrho,bkg
535	// scale = w[0]
536	// radius = w[1]
537	// delrho = w[2]
538	// bkg = w[3]
539
540
541	// calculates scale * f^2/Vol where f=Vol3delrho*((sin(qr)-qrcos(qr))/qr^3
542	// and is rescaled to give [=] cm^-1
543
544	Variable bes,f,vol,f2
545	//
546	//handle q==0 separately
547	If(x==0)
548	f = 4/3piradius^3delrhodelrhoscale1e8 + bkg
549	return(f)
550	Endif
551
552	// bes = 3(sin(xradius)-xradiuscos(x*radius))/x^3/radius^3
553
554	bes = 3sqrt(pi/(2xradius))BesselJ(1.5,xradius)/(xradius)
555
556	vol = 4pi/3radius^3
557	f = volbesdelrho // [=] A
558	// normalize to single particle volume, convert to 1/cm
559	f2 = f * f / vol * 1.0e8 // [=] 1/cm
560
561	return (scale*f2+bkg) // Scale, then add in the background
562
563	End
564
565	Function VC_Debye(scale,rg,bkg,x)
566	Variable scale,rg,bkg
567	Variable x
568
569	// variables are:
570	//[0] scale factor
571	//[1] radius of gyration [A]
572	//[2] background [cm-1]
573
574	// calculates (scale*debye)+bkg
575	Variable Pq,qr2
576
577	qr2=(x*rg)^2
578	Pq = 2*(exp(-(qr2))-1+qr2)/qr2^2
579
580	//scale
581	Pq *= scale
582	// then add in the background
583	return (Pq+bkg)
584	End
585
586	// a sum of a power law and debye to approximate the scattering from a real empty cell
587	//
588	// make/O/D coef_ECEmp = {2.2e-8,3.346,0.0065,9.0,0.016}
589	//
590	Function VC_EC_Empirical(aa,mm,scale,rg,bkg,x)
591	Variable aa,mm,scale,rg,bkg
592	Variable x
593
594	// variables are:
595	//[0] = A
596	//[1] = power m
597	//[2] scale factor
598	//[3] radius of gyration [A]
599	//[4] background [cm-1]
600
601	Variable Iq
602
603	// calculates (scale*debye)+bkg
604	Variable Pq,qr2
605
606	// if(x*Rg < 1e-3) //added Oct 2008 to avoid numerical errors at low arg values
607	// return(scale+bkg)
608	// endif
609
610	Iq = aa*x^-mm
611
612	qr2=(x*rg)^2
613	Pq = 2*(exp(-(qr2))-1+qr2)/qr2^2
614
615	//scale
616	Pq *= scale
617	// then add the terms up
618	return (Iq + Pq + bkg)
619	End
620
621	// blocked beam
622	//
623	Function VC_BlockedBeam(bkg,x)
624	Variable bkg
625	Variable x
626
627	return (bkg)
628	End
629
630
631	//
632	// a broad peak to simulate silver behenate or vycor
633	//
634	// peak @ 0.1 ~ AgBeh
635	// Make/O/D coef_BroadPeak = {1e-9, 3, 20, 100.0, 0.1,3,0.1}
636	//
637	//
638	// peak @ 0.015 in middle of middle detector, maybe not "real" vycor, but that is to be resolved
639	// Make/O/D coef_BroadPeak = {1e-9, 3, 20, 500.0, 0.015,3,0.1}
640	//
641	//
642	Function VC_BroadPeak(aa,nn,cc,LL,Qzero,mm,bgd,x)
643	Variable aa,nn,cc,LL,Qzero,mm,bgd
644	Variable x
645
646	// variables are:
647	//[0] Porod term scaling
648	//[1] Porod exponent
649	//[2] Lorentzian term scaling
650	//[3] Lorentzian screening length [A]
651	//[4] peak location [1/A]
652	//[5] Lorentzian exponent
653	//[6] background
654
655	// local variables
656	Variable inten, qval
657	// x is the q-value for the calculation
658	qval = x
659	// do the calculation and return the function value
660
661	inten = aa/(qval)^nn + cc/(1 + (abs(qval-Qzero)*LL)^mm) + bgd
662
663	Return (inten)
664
665	End
666
667	//
668	// updated to new folder structure Feb 2016
669	// folderStr = RAW,SAM, VCALC or other
670	// detStr is the panel identifer "ML", etc.
671	//
672	Function V_SetDeltaQ(folderStr,detStr)
673	String folderStr,detStr
674
675	Variable isVCALC
676	if(cmpstr(folderStr,"VCALC") == 0)
677	isVCALC = 1
678	endif
679
680	String folderPath = "root:Packages:NIST:VSANS:"+folderStr
681	String instPath = ":entry:instrument:detector_"
682
683	if(isVCALC)
684	WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr) // 2D detector data
685	else
686	Wave inten = V_getDetectorDataW(folderStr,detStr)
687	endif
688
689	Wave qx = $(folderPath+instPath+detStr+":qx_"+detStr)
690	Wave qy = $(folderPath+instPath+detStr+":qy_"+detStr)
691
692	Variable xDim,yDim,delQ
693
694	xDim=DimSize(inten,0)
695	yDim=DimSize(inten,1)
696
697	if(xDim<yDim)
698	delQ = abs(qx[0][0] - qx[1][0])/2
699	else
700	delQ = abs(qy[0][1] - qy[0][0])/2
701	endif
702
703	// set the global
704	Variable/G $(folderPath+instPath+detStr+":gDelQ_"+detStr) = delQ
705	// Print "SET delQ = ",delQ," for ",type
706
707	return(delQ)
708	end
709
710
711	//TODO -- need a switch here to dispatch to the averaging type
712	Proc VC_BinQxQy_to_1D(folderStr,type)
713	String folderStr
714	String type
715	// Prompt folderStr,"Pick the data folder containing 2D data",popup,getAList(4)
716	// Prompt type,"detector identifier"
717
718
719	VC_fDoBinning_QxQy2D(folderStr, type)
720
721
722	/// this is for a tall, narrow slit mode
723	// VC_fBinDetector_byRows(folderStr,type)
724
725	End
726
727
728	// folderStr is RAW, VCALC, SAM, etc.
729	// type is "B", "FL" for single binning, "FLR", or "MLRTB" or similar if multiple panels are combined
730	//
731	Proc VC_Graph_1D_detType(folderStr,type)
732	String folderStr,type
733
734	SetDataFolder $("root:Packages:NIST:VSANS:"+folderStr)
735
736	Display $("iBin_qxqy"+"_"+type) vs $("qBin_qxqy"+"_"+type)
737	ModifyGraph mirror=2,grid=1,log=1
738	ModifyGraph mode=4,marker=19,msize=2
739	// ErrorBars/T=0 iBin_qxqy Y,wave=(eBin2D_qxqy,eBin2D_qxqy) // for simulations, I don't have 2D uncertainty
740	ErrorBars/T=0 $("iBin_qxqy"+"_"+type) Y,wave=($("eBin_qxqy"+"_"+type),$("eBin_qxqy"+"_"+type))
741	legend
742
743	SetDataFolder root:
744
745	End
746
747
748
749	//////////
750	//
751	// Function that bins a 2D detctor panel into I(q) based on the q-value of the pixel
752	// - each pixel QxQyQz has been calculated beforehand
753	// - if multiple panels are selected to be combined, it is done here during the binning
754	// - the setting of deltaQ step is still a little suspect (TODO)
755	//
756	//
757	// see the equivalent function in PlotUtils2D_v40.ipf
758	//
759	//Function fDoBinning_QxQy2D(inten,qx,qy,qz)
760	//
761	// this has been modified to accept different detector panels and to take arrays
762	// -- type = FL or FR or...other panel identifiers
763	//
764	// TODO "iErr" is not always defined correctly since it doesn't really apply here for data that is not 2D simulation
765	//
766	//
767	// updated Feb2016 to take new folder structure
768	// TODO
769	// -- VERIFY
770	// -- figure out what the best location is to put the averaged data? currently @ top level of WORK folder
771	// but this is a lousy choice.
772	// x- binning is now Mask-aware. If mask is not present, all data is used. If data is from VCALC, all data is used
773	// x- Where do I put the solid angle correction? In here as a weight for each point, or later on as
774	// a blanket correction (matrix multiply) for an entire panel? (Solid Angle correction is done in the
775	// step where data is added to a WORK file (see Raw_to_Work())
776	//
777	//
778	// TODO:
779	// -- some of the input parameters for the resolution calcuation are either assumed (apOff) or are currently
780	// hard-wired. these need to be corrected before even the pinhole resolution is correct
781	// x- resolution calculation is in the correct place. The calculation is done per-panel (specified by TYPE),
782	// and then the unwanted points can be discarded (all 6 columns) as the data is trimmed and concatenated
783	// is separate functions that are resolution-aware.
784	//
785	//
786	// folderStr = WORK folder, type = the binning type (may include multiple detectors)
787	Function VC_fDoBinning_QxQy2D(folderStr,type)
788	String folderStr,type
789
790	Variable nSets = 0
791	Variable xDim,yDim
792	Variable ii,jj
793	Variable qVal,nq,var,avesq,aveisq
794	Variable binIndex,val,isVCALC=0,maskMissing
795
796	String folderPath = "root:Packages:NIST:VSANS:"+folderStr
797	String instPath = ":entry:instrument:detector_"
798	String detStr
799
800	if(cmpstr(folderStr,"VCALC") == 0)
801	isVCALC = 1
802	endif
803
804	detStr = type
805
806	// now switch on the type to determine which waves to declare and create
807	// since there may be more than one panel to step through. There may be two, there may be four
808	//
809
810	// TODO:
811	// -- Solid_Angle -- waves will be present for WORK data other than RAW, but not for RAW
812	//
813	// assume that the mask files are missing unless we can find them. If VCALC data,
814	// then the Mask is missing by definition
815	maskMissing = 1
816
817	strswitch(type) // string switch
818	// case "FL": // execute if case matches expression
819	// case "FR":
820	// detStr = type
821	// if(isVCALC)
822	// WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
823	// WAVE/Z iErr = $("iErr_"+detStr) // 2D errors -- may not exist, especially for simulation
824	// else
825	// Wave inten = V_getDetectorDataW(folderStr,detStr)
826	// Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
827	// Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
828	// if(WaveExists(mask) == 1)
829	// maskMissing = 0
830	// endif
831	//
832	// endif
833	// NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
834	// Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr) // 2D q-values
835	// nSets = 1
836	// break
837
838	// case "FT":
839	// case "FB":
840	// detStr = type
841	// if(isVCALC)
842	// WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
843	// WAVE/Z iErr = $("iErr_"+detStr) // 2D errors -- may not exist, especially for simulation
844	// else
845	// Wave inten = V_getDetectorDataW(folderStr,detStr)
846	// Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
847	// Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
848	// if(WaveExists(mask) == 1)
849	// maskMissing = 0
850	// endif
851	// endif
852	// NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
853	// Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr) // 2D q-values
854	// nSets = 1
855	// break
856
857	// case "ML":
858	// case "MR":
859	// detStr = type
860	// if(isVCALC)
861	// WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
862	// WAVE/Z iErr = $("iErr_"+detStr) // 2D errors -- may not exist, especially for simulation
863	// else
864	// Wave inten = V_getDetectorDataW(folderStr,detStr)
865	// Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
866	// Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
867	// if(WaveExists(mask) == 1)
868	// maskMissing = 0
869	// endif
870	// endif
871	// //TODO:
872	// // -- decide on the proper deltaQ for binning. either nominal value for LR, or one
873	// // determined specifically for that panel (currently using one tube width as deltaQ)
874	// // -- this is repeated multiple times in this switch
875	// NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
876	// Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr) // 2D q-values
877	// nSets = 1
878	// break
879
880	// case "MT":
881	// case "MB":
882	// detStr = type
883	// if(isVCALC)
884	// WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
885	// WAVE/Z iErr = $("iErr_"+detStr) // 2D errors -- may not exist, especially for simulation
886	// else
887	// Wave inten = V_getDetectorDataW(folderStr,detStr)
888	// Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
889	// Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
890	// if(WaveExists(mask) == 1)
891	// maskMissing = 0
892	// endif
893	// endif
894	// NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
895	// Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr) // 2D q-values
896	// nSets = 1
897	// break
898
899	// only one panel, simply pick that panel and move on out of the switch
900	case "FL":
901	case "FR":
902	case "FT":
903	case "FB":
904	case "ML":
905	case "MR":
906	case "MT":
907	case "MB":
908	case "B":
909	if(isVCALC)
910	WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
911	WAVE/Z iErr = $("iErr_"+detStr) // 2D errors -- may not exist, especially for simulation
912	else
913	Wave inten = V_getDetectorDataW(folderStr,detStr)
914	Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
915	Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
916	if(WaveExists(mask) == 1)
917	maskMissing = 0
918	endif
919	endif
920	NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
921	Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr) // 2D q-values
922	nSets = 1
923	break
924
925	case "FLR":
926	// detStr has multiple values now, so unfortuntely, I'm hard-wiring things...
927	// TODO
928	// -- see if I can un-hard-wire some of this below when more than one panel is combined
929	if(isVCALC)
930	WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
931	WAVE/Z iErr = $("iErr_"+"FL") // 2D errors -- may not exist, especially for simulation
932	WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
933	WAVE/Z iErr2 = $("iErr_"+"FR") // 2D errors -- may not exist, especially for simulation
934	else
935	Wave inten = V_getDetectorDataW(folderStr,"FL")
936	Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
937	Wave inten2 = V_getDetectorDataW(folderStr,"FR")
938	Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
939	Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data")
940	Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data")
941	if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
942	maskMissing = 0
943	endif
944	endif
945	NVAR delQ = $(folderPath+instPath+"FL"+":gDelQ_FL")
946
947	Wave qTotal = $(folderPath+instPath+"FL"+":qTot_"+"FL") // 2D q-values
948	Wave qTotal2 = $(folderPath+instPath+"FR"+":qTot_"+"FR") // 2D q-values
949
950	nSets = 2
951	break
952
953	case "FTB":
954	if(isVCALC)
955	WAVE inten = $(folderPath+instPath+"FT"+":det_"+"FT")
956	WAVE/Z iErr = $("iErr_"+"FT") // 2D errors -- may not exist, especially for simulation
957	WAVE inten2 = $(folderPath+instPath+"FB"+":det_"+"FB")
958	WAVE/Z iErr2 = $("iErr_"+"FB") // 2D errors -- may not exist, especially for simulation
959	else
960	Wave inten = V_getDetectorDataW(folderStr,"FT")
961	Wave iErr = V_getDetectorDataErrW(folderStr,"FT")
962	Wave inten2 = V_getDetectorDataW(folderStr,"FB")
963	Wave iErr2 = V_getDetectorDataErrW(folderStr,"FB")
964	Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FT"+":data")
965	Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FB"+":data")
966	if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
967	maskMissing = 0
968	endif
969	endif
970	NVAR delQ = $(folderPath+instPath+"FT"+":gDelQ_FT")
971
972	Wave qTotal = $(folderPath+instPath+"FT"+":qTot_"+"FT") // 2D q-values
973	Wave qTotal2 = $(folderPath+instPath+"FB"+":qTot_"+"FB") // 2D q-values
974
975	nSets = 2
976	break
977
978	case "FLRTB":
979	if(isVCALC)
980	WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
981	WAVE/Z iErr = $("iErr_"+"FL") // 2D errors -- may not exist, especially for simulation
982	WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
983	WAVE/Z iErr2 = $("iErr_"+"FR") // 2D errors -- may not exist, especially for simulation
984	WAVE inten3 = $(folderPath+instPath+"FT"+":det_"+"FT")
985	WAVE/Z iErr3 = $("iErr_"+"FT") // 2D errors -- may not exist, especially for simulation
986	WAVE inten4 = $(folderPath+instPath+"FB"+":det_"+"FB")
987	WAVE/Z iErr4 = $("iErr_"+"FB") // 2D errors -- may not exist, especially for simulation
988	else
989	Wave inten = V_getDetectorDataW(folderStr,"FL")
990	Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
991	Wave inten2 = V_getDetectorDataW(folderStr,"FR")
992	Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
993	Wave inten3 = V_getDetectorDataW(folderStr,"FT")
994	Wave iErr3 = V_getDetectorDataErrW(folderStr,"FT")
995	Wave inten4 = V_getDetectorDataW(folderStr,"FB")
996	Wave iErr4 = V_getDetectorDataErrW(folderStr,"FB")
997	Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data")
998	Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data")
999	Wave/Z mask3 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FT"+":data")
1000	Wave/Z mask4 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FB"+":data")
1001	if(WaveExists(mask) == 1 && WaveExists(mask2) == 1 && WaveExists(mask3) == 1 && WaveExists(mask4) == 1)
1002	maskMissing = 0
1003	endif
1004	endif
1005	NVAR delQ = $(folderPath+instPath+"FL"+":gDelQ_FL")
1006
1007	Wave qTotal = $(folderPath+instPath+"FL"+":qTot_"+"FL") // 2D q-values
1008	Wave qTotal2 = $(folderPath+instPath+"FR"+":qTot_"+"FR") // 2D q-values
1009	Wave qTotal3 = $(folderPath+instPath+"FT"+":qTot_"+"FT") // 2D q-values
1010	Wave qTotal4 = $(folderPath+instPath+"FB"+":qTot_"+"FB") // 2D q-values
1011
1012	nSets = 4
1013	break
1014
1015	case "MLR":
1016	if(isVCALC)
1017	WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
1018	WAVE/Z iErr = $("iErr_"+"ML") // 2D errors -- may not exist, especially for simulation
1019	WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
1020	WAVE/Z iErr2 = $("iErr_"+"MR") // 2D errors -- may not exist, especially for simulation
1021	else
1022	Wave inten = V_getDetectorDataW(folderStr,"ML")
1023	Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
1024	Wave inten2 = V_getDetectorDataW(folderStr,"MR")
1025	Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
1026	Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data")
1027	Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data")
1028	if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
1029	maskMissing = 0
1030	endif
1031	endif
1032	NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
1033
1034	Wave qTotal = $(folderPath+instPath+"ML"+":qTot_"+"ML") // 2D q-values
1035	Wave qTotal2 = $(folderPath+instPath+"MR"+":qTot_"+"MR") // 2D q-values
1036
1037	nSets = 2
1038	break
1039
1040	case "MTB":
1041	if(isVCALC)
1042	WAVE inten = $(folderPath+instPath+"MT"+":det_"+"MT")
1043	WAVE/Z iErr = $("iErr_"+"MT") // 2D errors -- may not exist, especially for simulation
1044	WAVE inten2 = $(folderPath+instPath+"MB"+":det_"+"MB")
1045	WAVE/Z iErr2 = $("iErr_"+"MB") // 2D errors -- may not exist, especially for simulation
1046	else
1047	Wave inten = V_getDetectorDataW(folderStr,"MT")
1048	Wave iErr = V_getDetectorDataErrW(folderStr,"MT")
1049	Wave inten2 = V_getDetectorDataW(folderStr,"MB")
1050	Wave iErr2 = V_getDetectorDataErrW(folderStr,"MB")
1051	Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MT"+":data")
1052	Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MB"+":data")
1053	if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
1054	maskMissing = 0
1055	endif
1056	endif
1057	NVAR delQ = $(folderPath+instPath+"MT"+":gDelQ_MT")
1058
1059	Wave qTotal = $(folderPath+instPath+"MT"+":qTot_"+"MT") // 2D q-values
1060	Wave qTotal2 = $(folderPath+instPath+"MB"+":qTot_"+"MB") // 2D q-values
1061
1062	nSets = 2
1063	break
1064
1065	case "MLRTB":
1066	if(isVCALC)
1067	WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
1068	WAVE/Z iErr = $("iErr_"+"ML") // 2D errors -- may not exist, especially for simulation
1069	WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
1070	WAVE/Z iErr2 = $("iErr_"+"MR") // 2D errors -- may not exist, especially for simulation
1071	WAVE inten3 = $(folderPath+instPath+"MT"+":det_"+"MT")
1072	WAVE/Z iErr3 = $("iErr_"+"MT") // 2D errors -- may not exist, especially for simulation
1073	WAVE inten4 = $(folderPath+instPath+"MB"+":det_"+"MB")
1074	WAVE/Z iErr4 = $("iErr_"+"MB") // 2D errors -- may not exist, especially for simulation
1075	else
1076	Wave inten = V_getDetectorDataW(folderStr,"ML")
1077	Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
1078	Wave inten2 = V_getDetectorDataW(folderStr,"MR")
1079	Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
1080	Wave inten3 = V_getDetectorDataW(folderStr,"MT")
1081	Wave iErr3 = V_getDetectorDataErrW(folderStr,"MT")
1082	Wave inten4 = V_getDetectorDataW(folderStr,"MB")
1083	Wave iErr4 = V_getDetectorDataErrW(folderStr,"MB")
1084	Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data")
1085	Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data")
1086	Wave/Z mask3 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MT"+":data")
1087	Wave/Z mask4 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MB"+":data")
1088	if(WaveExists(mask) == 1 && WaveExists(mask2) == 1 && WaveExists(mask3) == 1 && WaveExists(mask4) == 1)
1089	maskMissing = 0
1090	endif
1091	endif
1092	NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
1093
1094	Wave qTotal = $(folderPath+instPath+"ML"+":qTot_"+"ML") // 2D q-values
1095	Wave qTotal2 = $(folderPath+instPath+"MR"+":qTot_"+"MR") // 2D q-values
1096	Wave qTotal3 = $(folderPath+instPath+"MT"+":qTot_"+"MT") // 2D q-values
1097	Wave qTotal4 = $(folderPath+instPath+"MB"+":qTot_"+"MB") // 2D q-values
1098
1099	nSets = 4
1100	break
1101
1102	default:
1103	nSets = 0
1104	Print "ERROR ---- type is not recognized "
1105	endswitch
1106
1107	// Print "delQ = ",delQ," for ",type
1108
1109	if(nSets == 0)
1110	SetDataFolder root:
1111	return(0)
1112	endif
1113
1114
1115	// RAW data is currently read in and the 2D error wave is correctly generated
1116	// 2D error is propagated through all reduction steps, but I have not
1117	// verified that it is an exact duplication of the 1D error
1118	//
1119	//
1120	//
1121	// IF ther is no 2D error wave present for some reason, make a fake one
1122	if(WaveExists(iErr)==0 && WaveExists(inten) != 0)
1123	Duplicate/O inten,iErr
1124	Wave iErr=iErr
1125	// iErr = 1+sqrt(inten+0.75) // can't use this -- it applies to counts, not intensity (already a count rate...)
1126	iErr = sqrt(inten+0.75) // TODO -- here I'm just using some fictional value
1127	endif
1128	if(WaveExists(iErr2)==0 && WaveExists(inten2) != 0)
1129	Duplicate/O inten2,iErr2
1130	Wave iErr2=iErr2
1131	// iErr2 = 1+sqrt(inten2+0.75) // can't use this -- it applies to counts, not intensity (already a count rate...)
1132	iErr2 = sqrt(inten2+0.75) // TODO -- here I'm just using some fictional value
1133	endif
1134	if(WaveExists(iErr3)==0 && WaveExists(inten3) != 0)
1135	Duplicate/O inten3,iErr3
1136	Wave iErr3=iErr3
1137	// iErr3 = 1+sqrt(inten3+0.75) // can't use this -- it applies to counts, not intensity (already a count rate...)
1138	iErr3 = sqrt(inten3+0.75) // TODO -- here I'm just using some fictional value
1139	endif
1140	if(WaveExists(iErr4)==0 && WaveExists(inten4) != 0)
1141	Duplicate/O inten4,iErr4
1142	Wave iErr4=iErr4
1143	// iErr4 = 1+sqrt(inten4+0.75) // can't use this -- it applies to counts, not intensity (already a count rate...)
1144	iErr4 = sqrt(inten4+0.75) // TODO -- here I'm just using some fictional value
1145	endif
1146
1147	// TODO -- nq will need to be larger, once the back detector is installed
1148	//
1149	// note that the back panel of 320x320 (1mm res) results in 447 data points!
1150	// - so I upped nq to 600
1151
1152	if(cmpstr(type,"B") == 0)
1153	nq = 8000
1154	else
1155	nq=600
1156	endif
1157
1158	//****TODO**** -- where to put the averaged data -- right now, folderStr is forced to ""
1159	// SetDataFolder $("root:"+folderStr) //should already be here, but make sure...
1160	Make/O/D/N=(nq) $(folderPath+":"+"iBin_qxqy"+"_"+type)
1161	Make/O/D/N=(nq) $(folderPath+":"+"qBin_qxqy"+"_"+type)
1162	Make/O/D/N=(nq) $(folderPath+":"+"nBin_qxqy"+"_"+type)
1163	Make/O/D/N=(nq) $(folderPath+":"+"iBin2_qxqy"+"_"+type)
1164	Make/O/D/N=(nq) $(folderPath+":"+"eBin_qxqy"+"_"+type)
1165	Make/O/D/N=(nq) $(folderPath+":"+"eBin2D_qxqy"+"_"+type)
1166
1167	Wave iBin_qxqy = $(folderPath+":"+"iBin_qxqy_"+type)
1168	Wave qBin_qxqy = $(folderPath+":"+"qBin_qxqy"+"_"+type)
1169	Wave nBin_qxqy = $(folderPath+":"+"nBin_qxqy"+"_"+type)
1170	Wave iBin2_qxqy = $(folderPath+":"+"iBin2_qxqy"+"_"+type)
1171	Wave eBin_qxqy = $(folderPath+":"+"eBin_qxqy"+"_"+type)
1172	Wave eBin2D_qxqy = $(folderPath+":"+"eBin2D_qxqy"+"_"+type)
1173
1174
1175	// delQ = abs(sqrt(qx[2]^2+qy[2]^2+qz[2]^2) - sqrt(qx[1]^2+qy[1]^2+qz[1]^2)) //use bins of 1 pixel width
1176	// TODO: not sure if I want to set dQ in x or y direction...
1177	// the short dimension is the 8mm tubes, use this direction as dQ?
1178	// but don't use the corner of the detector, since dQ will be very different on T/B or L/R due to the location of [0,0]
1179	// WRT the beam center. use qx or qy directly. Still not happy with this way...
1180
1181
1182	qBin_qxqy[] = p*delQ
1183	SetScale/P x,0,delQ,"",qBin_qxqy //allows easy binning
1184
1185	iBin_qxqy = 0
1186	iBin2_qxqy = 0
1187	eBin_qxqy = 0
1188	eBin2D_qxqy = 0
1189	nBin_qxqy = 0 //number of intensities added to each bin
1190
1191	// now there are situations of:
1192	// 1 panel
1193	// 2 panels
1194	// 4 panels
1195	//
1196	// this needs to be a double loop now...
1197	// TODO:
1198	// -- the iErr (=2D) wave and accumulation of error is NOT CALCULATED CORRECTLY YET
1199	// -- verify the 2D error propagation by reducing it to 1D error
1200	//
1201	//
1202	// The 1D error does not use iErr, and IS CALCULATED CORRECTLY
1203	//
1204	// x- the solid angle per pixel will be present for WORK data other than RAW, but not for RAW
1205
1206	//
1207	// if any of the masks don't exist, display the error, and proceed with the averaging, using all data
1208	if(maskMissing == 1)
1209	Print "Mask file not found for at least one detector - so all data is used"
1210	endif
1211
1212	NVAR gIgnoreDetB = root:Packages:NIST:VSANS:Globals:gIgnoreDetB
1213	if(gIgnoreDetB && cmpstr(type,"B") == 0)
1214	maskMissing = 1
1215	Print "Mask skipped for B due to possible mismatch (Panel B ignored in preferences)"
1216	endif
1217
1218	Variable mask_val
1219	// use set 1 (no number) only
1220	if(nSets >= 1)
1221	xDim=DimSize(inten,0)
1222	yDim=DimSize(inten,1)
1223
1224	for(ii=0;ii<xDim;ii+=1)
1225	for(jj=0;jj<yDim;jj+=1)
1226	//qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1227	qVal = qTotal[ii][jj]
1228	binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1229	val = inten[ii][jj]
1230
1231	if(isVCALC \|\| maskMissing) // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
1232	mask_val = 0
1233	else
1234	mask_val = mask[ii][jj]
1235	endif
1236	if (numType(val)==0 && mask_val == 0) //count only the good points, ignore Nan or Inf
1237	iBin_qxqy[binIndex] += val
1238	iBin2_qxqy[binIndex] += val*val
1239	eBin2D_qxqy[binIndex] += iErr[ii][jj]*iErr[ii][jj]
1240	nBin_qxqy[binIndex] += 1
1241	endif
1242	endfor
1243	endfor
1244
1245	endif
1246
1247	// add in set 2 (set 1 already done)
1248	if(nSets >= 2)
1249	xDim=DimSize(inten2,0)
1250	yDim=DimSize(inten2,1)
1251
1252	for(ii=0;ii<xDim;ii+=1)
1253	for(jj=0;jj<yDim;jj+=1)
1254	//qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1255	qVal = qTotal2[ii][jj]
1256	binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1257	val = inten2[ii][jj]
1258
1259	if(isVCALC \|\| maskMissing)
1260	mask_val = 0
1261	else
1262	mask_val = mask2[ii][jj]
1263	endif
1264	if (numType(val)==0 && mask_val == 0) //count only the good points, ignore Nan or Inf
1265	iBin_qxqy[binIndex] += val
1266	iBin2_qxqy[binIndex] += val*val
1267	eBin2D_qxqy[binIndex] += iErr2[ii][jj]*iErr2[ii][jj]
1268	nBin_qxqy[binIndex] += 1
1269	endif
1270	endfor
1271	endfor
1272
1273	endif
1274
1275	// add in set 3 and 4 (set 1 and 2 already done)
1276	if(nSets == 4)
1277	xDim=DimSize(inten3,0)
1278	yDim=DimSize(inten3,1)
1279
1280	for(ii=0;ii<xDim;ii+=1)
1281	for(jj=0;jj<yDim;jj+=1)
1282	//qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1283	qVal = qTotal3[ii][jj]
1284	binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1285	val = inten3[ii][jj]
1286
1287	if(isVCALC \|\| maskMissing)
1288	mask_val = 0
1289	else
1290	mask_val = mask3[ii][jj]
1291	endif
1292	if (numType(val)==0 && mask_val == 0) //count only the good points, ignore Nan or Inf
1293	iBin_qxqy[binIndex] += val
1294	iBin2_qxqy[binIndex] += val*val
1295	eBin2D_qxqy[binIndex] += iErr3[ii][jj]*iErr3[ii][jj]
1296	nBin_qxqy[binIndex] += 1
1297	endif
1298	endfor
1299	endfor
1300
1301
1302	xDim=DimSize(inten4,0)
1303	yDim=DimSize(inten4,1)
1304
1305	for(ii=0;ii<xDim;ii+=1)
1306	for(jj=0;jj<yDim;jj+=1)
1307	//qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1308	qVal = qTotal4[ii][jj]
1309	binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1310	val = inten4[ii][jj]
1311
1312	if(isVCALC \|\| maskMissing)
1313	mask_val = 0
1314	else
1315	mask_val = mask4[ii][jj]
1316	endif
1317	if (numType(val)==0 && mask_val == 0) //count only the good points, ignore Nan or Inf
1318	iBin_qxqy[binIndex] += val
1319	iBin2_qxqy[binIndex] += val*val
1320	eBin2D_qxqy[binIndex] += iErr4[ii][jj]*iErr4[ii][jj]
1321	nBin_qxqy[binIndex] += 1
1322	endif
1323	endfor
1324	endfor
1325
1326	endif
1327
1328
1329	// after looping through all of the data on the panels, calculate errors on I(q),
1330	// just like in CircSectAve.ipf
1331	// TODO:
1332	// -- 2D Errors were (maybe) properly acculumated through reduction, so this loop of calculations is NOT VERIFIED (yet)
1333	// x- the error on the 1D intensity, is correctly calculated as the standard error of the mean.
1334	for(ii=0;ii<nq;ii+=1)
1335	if(nBin_qxqy[ii] == 0)
1336	//no pixels in annuli, data unknown
1337	iBin_qxqy[ii] = 0
1338	eBin_qxqy[ii] = 1
1339	eBin2D_qxqy[ii] = NaN
1340	else
1341	if(nBin_qxqy[ii] <= 1)
1342	//need more than one pixel to determine error
1343	iBin_qxqy[ii] /= nBin_qxqy[ii]
1344	eBin_qxqy[ii] = 1
1345	eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1346	else
1347	//assume that the intensity in each pixel in annuli is normally distributed about mean...
1348	// -- this is correctly calculating the error as the standard error of the mean, as
1349	// was always done for SANS as well.
1350	iBin_qxqy[ii] /= nBin_qxqy[ii]
1351	avesq = iBin_qxqy[ii]^2
1352	aveisq = iBin2_qxqy[ii]/nBin_qxqy[ii]
1353	var = aveisq-avesq
1354	if(var<=0)
1355	eBin_qxqy[ii] = 1e-6
1356	else
1357	eBin_qxqy[ii] = sqrt(var/(nBin_qxqy[ii] - 1))
1358	endif
1359	// and calculate as it is propagated pixel-by-pixel
1360	eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1361	endif
1362	endif
1363	endfor
1364
1365	eBin2D_qxqy = sqrt(eBin2D_qxqy) // as equation (3) of John's memo
1366
1367	// find the last non-zero point, working backwards
1368	val=nq
1369	do
1370	val -= 1
1371	while((nBin_qxqy[val] == 0) && val > 0)
1372
1373	// print val, nBin_qxqy[val]
1374	DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1375
1376	if(val == 0)
1377	// all the points were deleted
1378	return(0)
1379	endif
1380
1381
1382	// since the beam center is not always on the detector, many of the low Q bins will have zero pixels
1383	// find the first non-zero point, working forwards
1384	val = -1
1385	do
1386	val += 1
1387	while(nBin_qxqy[val] == 0)
1388	DeletePoints 0, val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1389
1390	// ?? there still may be a point in the q-range that gets zero pixel contribution - so search this out and get rid of it
1391	val = numpnts(nBin_qxqy)-1
1392	do
1393	if(nBin_qxqy[val] == 0)
1394	DeletePoints val, 1, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1395	endif
1396	val -= 1
1397	while(val>0)
1398
1399	// TODO:
1400	// -- This is where I calculate the resolution in SANS (see CircSectAve)
1401	// -- use the isVCALC flag to exclude VCALC from the resolution calculation if necessary
1402	// -- from the top of the function, folderStr = work folder, type = "FLRTB" or other type of averaging
1403	//
1404	nq = numpnts(qBin_qxqy)
1405	Make/O/D/N=(nq) $(folderPath+":"+"sigmaQ_qxqy"+"_"+type)
1406	Make/O/D/N=(nq) $(folderPath+":"+"qBar_qxqy"+"_"+type)
1407	Make/O/D/N=(nq) $(folderPath+":"+"fSubS_qxqy"+"_"+type)
1408	Wave sigmaq = $(folderPath+":"+"sigmaQ_qxqy_"+type)
1409	Wave qbar = $(folderPath+":"+"qBar_qxqy_"+type)
1410	Wave fsubs = $(folderPath+":"+"fSubS_qxqy_"+type)
1411
1412
1413	Variable ret1,ret2,ret3
1414	Variable lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses
1415
1416	// TODO: check the units of all of the inputs
1417
1418	// lambda = wavelength [A]
1419	lambda = V_getWavelength(folderStr)
1420
1421	// lambdaWidth = [dimensionless]
1422	lambdaWidth = V_getWavelength_spread(folderStr)
1423
1424	// DDet = detector pixel resolution [cm] **assumes square pixel
1425	DDet = 0.8 // TODO -- this is hard-wired
1426
1427	// apOff = sample aperture to sample distance [cm]
1428	apOff = 10 // TODO -- this is hard-wired
1429
1430	// S1 = source aperture diameter [mm]
1431	S1 = str2num(V_getSourceAp_size(folderStr))
1432
1433	// S2 = sample aperture diameter [cm]
1434	// as of 3/2018, the "internal" sample aperture is not in use, only the external
1435	// TODO : verify the units on the Ap2 (external)
1436	// sample aperture 1(internal) is set to report "12.7 mm" as a STRING
1437	// sample aperture 2(external) reports the number typed in...
1438	//
1439	// so I'm trusting [cm] is in the file
1440	S2 = V_getSampleAp2_size(folderStr)*10 // sample ap 1 or 2? 2 = the "external", convert to [mm]
1441
1442	// L1 = source to sample distance [m]
1443	L1 = V_getSourceAp_distance(folderStr)/100
1444
1445	// L2 = sample to detector distance [m]
1446	// take the first two characters of the "type" to get the correct distance.
1447	// if the type is say, MLRTB, then the implicit assumption in combining all four panels is that the resolution
1448	// is not an issue for the slightly different distances.
1449	L2 = V_getDet_ActualDistance(folderStr,type[0,1])/100 //convert cm to m
1450
1451	// BS = beam stop diameter [mm]
1452	//TODO:? which BS is in? carr2, carr3, none?
1453	// -- need to check the num_beamstops field, then description, then shape/size or shape/height and shape/width
1454	//
1455	// TODO: the values in the file are incorrect!!! BS = 1000 mm diameter!!!
1456	// BS = V_getBeamStopC2_size(folderStr) // Units are [mm]
1457	BS = 25.4 //TODO hard-wired value
1458
1459	// del_r = step size [mm] = binWidth*(mm/pixel)
1460	del_r = 1*8 // TODO: 8mm/pixel hard-wired
1461
1462	// usingLenses = flag for lenses = 0 if no lenses, non-zero if lenses are in-beam
1463	usingLenses = 0
1464
1465	if(cmpstr(detStr,"FL")==0)
1466	Print "(FL) Resolution lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses"
1467	Print lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses
1468	endif
1469
1470
1471	// TODO:
1472	// this is the point where I need to switch on the different collimation types (white beam, slit, Xtal, etc)
1473	// to calculate the correct resolution, or fill the waves with the correct "flags"
1474	//
1475
1476	// For white beam data, the wavelength distribution can't be represented as a gaussian, but all of the other
1477	// geometric corrections still apply. Passing zero for the lambdaWidth will return the geometry contribution,
1478	// as long as the wavelength can be handled separately. It appears to be correct to do as a double integral,
1479	// with the inner(lambda) calculated first, then the outer(geometry).
1480	//
1481
1482
1483	ii=0
1484	do
1485	V_getResolution(qBin_qxqy[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,ret1,ret2,ret3)
1486	sigmaq[ii] = ret1
1487	qbar[ii] = ret2
1488	fsubs[ii] = ret3
1489	ii+=1
1490	while(ii<nq)
1491
1492
1493
1494	SetDataFolder root:
1495
1496	return(0)
1497	End
1498
1499

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

Download in other formats:

Original Format