Context Navigation

VC_DetectorBinning_Utils.ipf @ 1117

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

extensive changes to accomodate 1x1 binning of the HighRes? detector. It is implemented as a global flag. Currently only 4x4 and 1x1 are allowed. 1x1 has never been tested in reality, only simulated data - so my assumed dimensions may not be correct. look for TODOHIGHRES in the file for places that may need to be updated for different file dimensions. Testing of the simulated data is proving to be excruciatingly slow, but passable for a test. Speed optimization will be needed if this is the final solution. Memory management will also be an issue since every "copy" of the highRes matrix is enormous. Carry as few of these around as possible in the future to keep the experiment size to something reasonable.

File size: 51.7 KB

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

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

Download in other formats:

Original Format