Context Navigation

VC_DetectorBinning_Utils.ipf @ 1249

Last change on this file since 1249 was 1249, checked in by srkline, 2 years ago
minor changes - can't find the diff command in TortoiseSVN...
File size: 47.9 KB

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

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

Download in other formats:

Original Format