Context Navigation

V_WorkFolderUtils.ipf @ 1055

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

LOTS of changes to accommodate the beam center being reported in cm rather than pixels. Required a lot of changes to VCALC (to fill in simulated data), and in the reading and interpreting of data for display, and most importantly, the calculation of q.

There may still be a few residual bugs with this. I am still re-testing with new sample data sets.

File size: 35.0 KB

Line
1	#pragma rtGlobals=3 // Use modern global access method and strict wave access.
2	#pragma version=1.0
3	#pragma IgorVersion=6.1
4
5	//*******************
6	// Vers 1.0 JAN2016
7	//
8	//*******************
9	// VSANS Utility procedures for handling of workfiles (each is housed in a separate datafolder)
10	//
11	// - adding RAW data to a workfile
12	// -- this conversion applies the detector corrections
13	// -- Raw_to_work(newType) IS THE MAJOR ROUTINE TO APPLY DETECTOR CORRECTIONS
14	//
15	//
16	// - copying workfiles to another folder
17	//
18	// - absolute scaling
19	//
20	// - (no) the WorkFile Math panel for simple image math (not done - maybe in the future?)
21	// -
22	// - (no) adding work.drk data without normalizing to monitor counts (the case not currently handled)
23	//***************************
24
25	//
26	// Functions used for manipulation of the local Igor "WORK" folder
27	// structure as raw data is displayed and processed.
28	//
29	//
30
31
32	//
33	//Entry procedure from main panel
34	//
35	Proc V_CopyWorkFolder(oldType,newType)
36	String oldType,newType
37	Prompt oldType,"Source WORK data type",popup,"SAM;EMP;BGD;DIV;COR;CAL;RAW;ABS;STO;SUB;DRK;"
38	Prompt newType,"Destination WORK data type",popup,"SAM;EMP;BGD;DIV;COR;CAL;RAW;ABS;STO;SUB;DRK;"
39
40	// data folder "old" will be copied to "new" (either kills/copies or will overwrite)
41	V_CopyHDFToWorkFolder(oldtype,newtype)
42	End
43
44	//
45	// copy what is needed for data processing (not the DAS_logs)
46	// from the RawVSANS storage folder to the local WORK folder as needed
47	//
48	// TODO -- at what stage do I make copies of data in linear/log forms for data display?
49	// -- do I need to make copies, if I'm displaying with the lookup wave (no copy needed if this works)
50	// -- when do I make the 2D error waves?
51	//
52	// TODO - decide what exactly I need to copy over. May be best to copy all, and delete
53	// what I know that I don't need
54	//
55	// TODO !!! DuplicateDataFolder will FAIL - in the base case of RAW data files, the
56	// data is actually in use - so it will fail every time. need an alternate solution. in SANS,
57	// there are a limited number of waves to carry over, so Dupliate/O is used for rw, tw, data, etc.
58	//
59	// TODO : I also need a list of what is generated during processing that may be hanging around - that I need to
60	// be sure to get rid of - like the calibration waves, solidAngle, etc.
61	//
62	// hdfDF is the name only of the data in storage. May be full file name with extension (clean as needed)
63	// type is the destination WORK folder for the copy
64	//
65	Function V_CopyHDFToWorkFolder(fromStr,toStr)
66	String fromStr,toStr
67
68	String fromDF, toDF
69
70	// make the DF paths - source and destination
71	fromDF = "root:Packages:NIST:VSANS:"+fromStr
72	toDF = "root:Packages:NIST:VSANS:"+toStr
73
74	// // make a copy of the file name for my own use, since it's not in the file
75	// String/G $(toDF+":file_name") = root:
76
77	// copy the folders
78	KillDataFolder/Z $toDF //DuplicateDataFolder will not overwrite, so Kill
79
80	if(V_flag == 0) // kill DF was OK
81	DuplicateDataFolder $("root:Packages:NIST:VSANS:"+fromStr),$("root:Packages:NIST:VSANS:"+toStr)
82
83	// I can delete these if they came along with RAW
84	// DAS_logs
85	// top-level copies of data (duplicate links, these should not be present in a proper NICE file)
86	KillDataFolder/Z $(toDF+":entry:DAS_logs")
87	KillDataFolder/Z $(toDF+":entry:data")
88	KillDataFolder/Z $(toDF+":entry:data_B")
89	KillDataFolder/Z $(toDF+":entry:data_ML")
90	KillDataFolder/Z $(toDF+":entry:data_MR")
91	KillDataFolder/Z $(toDF+":entry:data_MT")
92	KillDataFolder/Z $(toDF+":entry:data_MB")
93	KillDataFolder/Z $(toDF+":entry:data_FL")
94	KillDataFolder/Z $(toDF+":entry:data_FR")
95	KillDataFolder/Z $(toDF+":entry:data_FT")
96	KillDataFolder/Z $(toDF+":entry:data_FB")
97
98	return(0)
99	else
100
101	V_KillWavesFullTree($toDF,toStr,0,"",1) // this will traverse the whole tree, trying to kill what it can
102
103	if(DataFolderExists("root:Packages:NIST:VSANS:"+toStr) == 0) //if the data folder (RAW, SAM, etc.) was just killed?
104	NewDataFolder/O $("root:Packages:NIST:VSANS:"+toStr)
105	endif
106
107	// need to do this the hard way, duplicate/O recursively
108	// see V_CopyToWorkFolder()
109
110	// everything on the top level
111	V_DuplicateDataFolder($(fromDF+":entry"),fromStr,toStr,0,"",0) //no recursion here
112	// control
113	V_DuplicateDataFolder($(fromDF+":entry:control"),fromStr,toStr,0,"",1) //yes recursion here
114	// instrument
115	V_DuplicateDataFolder($(fromDF+":entry:instrument"),fromStr,toStr,0,"",1) //yes recursion here
116	// reduction
117	V_DuplicateDataFolder($(fromDF+":entry:reduction"),fromStr,toStr,0,"",1) //yes recursion here
118	// sample
119	V_DuplicateDataFolder($(fromDF+":entry:sample"),fromStr,toStr,0,"",1) //yes recursion here
120	// user
121	V_DuplicateDataFolder($(fromDF+":entry:user"),fromStr,toStr,0,"",1) //yes recursion here
122
123	endif
124
125	return(0)
126	end
127
128
129
130	Function V_KillWavesInFolder(folderStr)
131	String folderStr
132
133	if(DataFolderExists(folderStr) && strlen(folderStr) != 0)
134	SetDataFolder $folderStr
135	KillWaves/A/Z
136	endif
137
138	SetDataFolder root:
139	return(0)
140	end
141
142	////////
143	// see the help entry for IndexedDir for help on (possibly) how to do this faster
144	// -- see Function ScanDirectories(pathName, printDirNames)
145	//
146
147
148	// from IgorExchange On July 17th, 2011 jtigor
149	// started from "Recursively List Data Folder Contents"
150	// Posted July 15th, 2011 by hrodstein
151	//
152	//
153	//
154	Proc V_CopyWorkFolderProc(dataFolderStr, fromStr, toStr, level, sNBName, recurse)
155	String dataFolderStr="root:Packages:NIST:VSANS:RAW"
156	String fromStr = "RAW"
157	String toStr="SAM"
158	Variable level=0
159	String sNBName="DataFolderTree"
160	Variable recurse = 1
161
162	V_DuplicateDataFolder($dataFolderStr, fromStr, toStr, level, sNBName, recurse)
163
164
165	end
166
167	// ListDataFolder(dfr, level)
168	// Recursively lists objects in data folder.
169	// Pass data folder path for dfr and 0 for level.
170	// pass level == 0 for the first call
171	// sNBName = "" prints nothing. any name will generate a notebook
172	//
173	// recurse == 0 will do only the specified folder, anything else will recurse all levels
174	// toStr is the string name of the top-level folder only, not the full path
175	//
176	//
177	Function V_DuplicateDataFolder(dfr, fromStr, toStr, level, sNBName,recurse)
178	DFREF dfr
179	String fromStr
180	String toStr
181	Variable level // Pass 0 to start
182	String sNBName
183	Variable recurse
184
185	String name
186	String dfName
187	String sString
188
189	String toDF = ""
190
191	if (level == 0) // this is the data folder, generate if needed in the destination
192	name = GetDataFolder(1, dfr)
193	// sPrintf sString, "%s (data folder)\r", name
194	toDF = ReplaceString(fromStr,name,toStr,1) // case-sensitive replace
195	sprintf sString, "NewDataFolder/O %s\r",toDF
196	NewDataFolder/O $(RemoveEnding(toDF,":")) // remove trailing semicolon if it's there
197
198	V_WriteBrowserInfo_test(sString, 1, sNBName)
199	endif
200
201	dfName = GetDataFolder(1, dfr)
202	toDF = ReplaceString(fromStr,dfName,toStr,1) // case-sensitive replace
203	Variable i
204
205	String indentStr = "\t"
206	for(i=0; i<level; i+=1)
207	indentStr += "\t"
208	endfor
209
210	Variable numWaves = CountObjectsDFR(dfr, 1)
211	for(i=0; i<numWaves; i+=1)
212	name = GetIndexedObjNameDFR(dfr, 1, i)
213	//
214	// wave type does not matter now. Duplicate does not care
215	//
216	sPrintf sString, "Duplicate/O %s, %s\r",dfName+name,toDF+name
217	Duplicate/O $(dfName+name),$(toDF+name)
218
219	V_WriteBrowserInfo_test(sString, 2, sNBName)
220	endfor
221
222	Variable numNumericVariables = CountObjectsDFR(dfr, 2)
223	for(i=0; i<numNumericVariables; i+=1)
224	name = GetIndexedObjNameDFR(dfr, 2, i)
225	sPrintf sString, "%s%s (numeric variable)\r", indentStr, name
226	V_WriteBrowserInfo_test(sString, 3, sNBName)
227	endfor
228
229	Variable numStringVariables = CountObjectsDFR(dfr, 3)
230	for(i=0; i<numStringVariables; i+=1)
231	name = GetIndexedObjNameDFR(dfr, 3, i)
232	sPrintf sString, "%s%s (string variable)\r", indentStr, name
233	V_WriteBrowserInfo_test(sString, 4, sNBName)
234	endfor
235
236	if(recurse)
237	Variable numDataFolders = CountObjectsDFR(dfr, 4)
238	for(i=0; i<numDataFolders; i+=1)
239	name = GetIndexedObjNameDFR(dfr, 4, i)
240	// sPrintf sString, "%s%s (data folder)\r", indentStr, name
241	dfName = GetDataFolder(1, dfr)
242
243	toDF = ReplaceString(fromStr,dfName,toStr,1) // case-sensitive replace
244	sprintf sString, "NewDataFolder/O %s\r",toDF+name
245	NewDataFolder/O $(toDF+name)
246
247
248	V_WriteBrowserInfo_test(sString, 1, sNBName)
249	DFREF childDFR = dfr:$(name)
250	V_DuplicateDataFolder(childDFR, fromStr, toStr, level+1, sNBName, recurse)
251	endfor
252	endif
253
254
255	End
256
257
258	// ListDataFolder(dfr, level)
259	// Recursively lists objects in data folder.
260	// Pass data folder path for dfr and 0 for level.
261	// pass level == 0 for the first call
262	// sNBName = "" prints nothing. any name will generate a notebook
263	//
264	// recurse == 0 will do only the specified folder, anything else will recurse all levels
265	// toStr is the string name of the top-level folder only, not the full path
266	//
267	//
268	Function V_KillWavesFullTree(dfr, fromStr, level, sNBName,recurse)
269	DFREF dfr
270	String fromStr
271	// String toStr
272	Variable level // Pass 0 to start
273	String sNBName
274	Variable recurse
275
276	String name
277	String dfName
278	String sString
279
280	String toDF = ""
281
282	if (level == 0) // this is the data folder, generate if needed in the destination
283	name = GetDataFolder(1, dfr)
284	sPrintf sString, "%s (data folder)\r", name
285	// toDF = ReplaceString(fromStr,name,toStr,1) // case-sensitive replace
286	// sprintf sString, "NewDataFolder/O %s\r",toDF
287	// NewDataFolder/O $(RemoveEnding(toDF,":")) // remove trailing semicolon if it's there
288
289	V_WriteBrowserInfo_test(sString, 1, sNBName)
290	endif
291
292	dfName = GetDataFolder(1, dfr)
293	// toDF = ReplaceString(fromStr,dfName,toStr,1) // case-sensitive replace
294	Variable i
295
296	String indentStr = "\t"
297	for(i=0; i<level; i+=1)
298	indentStr += "\t"
299	endfor
300
301	Variable numWaves = CountObjectsDFR(dfr, 1)
302	for(i=0; i<numWaves; i+=1)
303	name = GetIndexedObjNameDFR(dfr, 1, i)
304	//
305	// wave type does not matter now. Kill does not care
306	//
307	sPrintf sString, "Killing %s\r",dfName+name
308	KillWaves/Z $(dfName+name)
309
310	V_WriteBrowserInfo_test(sString, 2, sNBName)
311	endfor
312
313	// now kill the data folder if possible
314	KillDataFolder/Z $dfName
315
316
317	Variable numNumericVariables = CountObjectsDFR(dfr, 2)
318	for(i=0; i<numNumericVariables; i+=1)
319	name = GetIndexedObjNameDFR(dfr, 2, i)
320	sPrintf sString, "%s%s (numeric variable)\r", indentStr, name
321	V_WriteBrowserInfo_test(sString, 3, sNBName)
322	endfor
323
324	Variable numStringVariables = CountObjectsDFR(dfr, 3)
325	for(i=0; i<numStringVariables; i+=1)
326	name = GetIndexedObjNameDFR(dfr, 3, i)
327	sPrintf sString, "%s%s (string variable)\r", indentStr, name
328	V_WriteBrowserInfo_test(sString, 4, sNBName)
329	endfor
330
331	if(recurse)
332	Variable numDataFolders = CountObjectsDFR(dfr, 4)
333	for(i=0; i<numDataFolders; i+=1)
334	name = GetIndexedObjNameDFR(dfr, 4, i)
335	sPrintf sString, "%s%s (data folder)\r", indentStr, name
336	dfName = GetDataFolder(1, dfr)
337
338	// toDF = ReplaceString(fromStr,dfName,toStr,1) // case-sensitive replace
339	// sprintf sString, "NewDataFolder/O %s\r",toDF+name
340	// NewDataFolder/O $(toDF+name)
341
342
343	V_WriteBrowserInfo_test(sString, 1, sNBName)
344	DFREF childDFR = dfr:$(name)
345	V_KillWavesFullTree(childDFR, fromStr, level+1, sNBName, recurse)
346	endfor
347	endif
348
349
350	End
351
352
353
354	Function V_WriteBrowserInfo_test(sString, vType, sNBName)
355	String sString
356	Variable vType
357	String sNBName
358
359	if(strlen(sNBName) == 0)
360	// print sString
361	return 0
362	endif
363	DoWindow $sNBName
364	if(V_flag != 1)
365	NewNoteBook/F=0 /N=$sNBName /V=1 as sNBName
366	else
367	DoWindow/F $sNBName
368	endif
369	Notebook $sNBName selection={endOfFile, endOfFile}
370	if(vType == 1) // a data folder
371	// Notebook $sNBName fstyle=1
372	Notebook $sNBName text=sString
373	// Notebook $sNBName fstyle=-1
374	else
375	Notebook $sNBName text=sString
376	endif
377
378	End
379
380	///////////////////////////////
381
382
383	//
384	// given the folder, duplicate the data -> linear_data and generate the error
385	// TODO
386	// -- do I want to use different names here? If it turns out that I don't need to drag a copy of
387	// the data around as "linear_data", then I can eliminate that, and rename the error wave
388	// -- be sure the data is either properly written as 2D in the file, or converted to 2D before
389	// duplicating here
390	// -- ? do I recast to DP here. Probably necessary since I'm doing a DP calculation, but Redimension
391	// is done in the Raw_to_Work step too. very confusing.
392	Function V_MakeDataError(folderStr)
393	String folderStr
394
395	SetDataFolder $folderStr
396	Wave data=data
397	Duplicate/O data linear_data // at this point, the data is still the raw data, and is linear_data
398
399	// proper error for counting statistics, good for low count values too
400	// rather than just sqrt(n)
401	// see N. Gehrels, Astrophys. J., 303 (1986) 336-346, equation (7)
402	// for S = 1 in eq (7), this corresponds to one sigma error bars
403	Duplicate/O linear_data linear_data_error
404	linear_data_error = 1 + sqrt(linear_data + 0.75)
405	//
406
407	SetDataFolder root:
408	return(0)
409	End
410
411
412	/////////////////////
413
414
415
416	//testing procedure
417	// TODO -- can't duplicate this with another proceudre, but if I change the name of the variable
418	// "newType" to "type", then when Raw_to_work() gets to CopyHDFToWorkFolder(), the KillDataFolder/Z
419	// line fails (but reports no error), then DuplicateDataFolder fails, and reports an error. Trying
420	// to simplify this condition, I can't reproduce the bug for WM...
421	Proc V_Convert_to_Workfile(newtype, doadd)
422	String newtype,doadd
423	Prompt newtype,"WORK data type",popup,"SAM;EMP;BGD;ADJ;"
424	Prompt doadd,"Add to current WORK contents?",popup,"No;Yes;"
425
426	//macro will take whatever is in RAW folder and "ADD" it to the folder specified
427	//in the popup menu
428
429	//"add" = yes/no, don't add to previous runs
430	//switch here - two separate functions to avoid (my) confusion
431	Variable err// = Raw_to_work(newtype)
432	if(cmpstr(doadd,"No")==0)
433	//don't add to prev work contents, copy RAW contents to work and convert
434	err = V_Raw_to_work(newtype)
435	else
436	//yes, add RAW to the current work folder contents
437	Abort "Adding RAW data files is currently unsupported"
438	err = V_Add_raw_to_work(newtype)
439	endif
440
441	String newTitle = "WORK_"+newtype
442	DoWindow/F VSANS_Data
443	DoWindow/T VSANS_Data, newTitle
444	KillStrings/Z newTitle
445
446	//need to update the display with "data" from the correct dataFolder
447	V_UpdateDisplayInformation(newtype)
448
449	End
450
451
452	//
453	// THIS IS THE MAJOR ROUTINE TO APPLY DATA CORRECTIONS
454	//
455	//will copy the current contents of the RAW folder to the newType work folder
456	//and do the geometric corrections and normalization to monitor counts
457	//(the function Add_Raw_to_work(type) adds multiple runs together - and is LOW priority)
458	//
459	//the current display type is updated to newType (global)
460	//
461	Function V_Raw_to_work(newType)
462	String newType
463
464	Variable deadTime,defmon,total_mon,total_det,total_trn,total_numruns,total_rtime
465	Variable ii,jj,itim,cntrate,dscale,scale,uscale
466	String destPath
467
468	String fname = newType
469	String detStr
470	Variable ctTime
471
472	//initialize values before normalization
473	total_mon=0
474	total_det=0
475	total_trn=0
476	total_numruns=0
477	total_rtime=0
478
479	//Not adding multiple runs, so wipe out the old contents of the work folder and
480	// replace with the contents of raw
481
482	destPath = "root:Packages:NIST:VSANS:" + newType
483
484	//copy from current dir (RAW) to work, defined by newType
485	V_CopyHDFToWorkFolder("RAW",newType)
486
487	// now work with the waves from the destination folder.
488
489	// apply corrections ---
490	// switches to control what is done, don't do the transmission correction for the BGD measurement
491	// start with the DIV correction, before conversion to mm
492	// then do all of the other corrections, order doesn't matter.
493	// rescaling to default monitor counts however, must be LAST.
494
495	// each correction must loop over each detector. tedious.
496
497	// (0) Redimension the data waves in the destination folder
498	// so that they are DP, not integer
499	// TODO
500	// -- currently only redimensioning the data and linear_data_error - What else???
501	// -- ?? some of this is done at load time for RAW data. shouldn't be an issue to re-do the redimension
502	for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
503	detStr = StringFromList(ii, ksDetectorListAll, ";")
504	Wave w = V_getDetectorDataW(fname,detStr)
505	Wave w_err = V_getDetectorDataErrW(fname,detStr)
506	Redimension/D w,w_err
507	endfor
508
509
510	// (1) DIV correction
511	// do this in terms of pixels.
512	// TODO : This must also exist at the time the first work folder is generated.
513	// So it must be in the user folder at the start of the experiment, and defined.
514	NVAR gDoDIVCor = root:Packages:NIST:VSANS:Globals:gDoDIVCor
515	if (gDoDIVCor == 1)
516	// need extra check here for file existence
517	// if not in DIV folder, load.
518	// if unable to load, skip correction and report error (Alert?) (Ask to Load?)
519	Print "Doing DIV correction"// for "+ detStr
520	for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
521	detStr = StringFromList(ii, ksDetectorListAll, ";")
522	Wave w = V_getDetectorDataW(fname,detStr)
523	Wave w_err = V_getDetectorDataErrW(fname,detStr)
524
525	V_DIVCorrection(w,w_err,detStr,newType)
526	endfor
527	else
528	Print "DIV correction NOT DONE" // not an error since correction was unchecked
529	endif
530
531	// (2) non-linear correction
532	// TODO:
533	// x- the "B" detector is calculated in its own routines
534	// -- document what is generated here:
535	// **in each detector folder: data_realDistX and data_realDistY (2D waves of the mm? position of the pixel)
536	// x- these spatial calculations ARE DONE as the RAW data is loaded. It allows the RAW
537	// data to be properly displayed, but without all of the (complete) set of detector corrections
538	// * the corrected distances are calculated into arrays, but nothing is done with them yet
539	// * there is enough information now to calculate the q-arrays, so it is done now
540	// - other corrections may modify the data, this calculation does NOT modify the data
541	NVAR gDoNonLinearCor = root:Packages:NIST:VSANS:Globals:gDoNonLinearCor
542	// generate a distance matrix for each of the detectors
543	if (gDoNonLinearCor == 1)
544	Print "Doing Non-linear correction"// for "+ detStr
545	for(ii=0;ii<ItemsInList(ksDetectorListNoB);ii+=1)
546	detStr = StringFromList(ii, ksDetectorListNoB, ";")
547	Wave w = V_getDetectorDataW(fname,detStr)
548	// Wave w_err = V_getDetectorDataErrW(fname,detStr)
549	Wave w_calib = V_getDetTube_spatialCalib(fname,detStr)
550	Variable tube_width = V_getDet_tubeWidth(fname,detStr)
551	V_NonLinearCorrection(fname,w,w_calib,tube_width,detStr,destPath)
552
553	//(2.4) Convert the beam center values from pixels to mm
554	// TODO -- there needs to be a permanent location for these values??
555	//
556	// TODO
557	// -- the beam center value in mm needs to be present - it is used in calculation of Qvalues
558	// -- but having both the same is wrong...
559	// -- the pixel value is needed for display of the panels
560	if(kBCTR_CM)
561	//V_ConvertBeamCtr_to_mm(folder,detStr,destPath)
562	//
563
564	Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
565	Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
566	WAVE x_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
567	WAVE y_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
568	x_mm[0] = V_getDet_beam_center_x(fname,detStr) * 10 // convert cm to mm
569	y_mm[0] = V_getDet_beam_center_y(fname,detStr) * 10 // convert cm to mm
570
571	// TODO:::
572	// now I need to convert the beam center in mm to pixels
573	// and have some rational place to look for it...
574	V_ConvertBeamCtr_to_pix(fname,detStr,destPath)
575	else
576	// beam center is in pixels, so use the old routine
577	V_ConvertBeamCtr_to_mm(fname,detStr,destPath)
578	endif
579
580	// (2.5) Calculate the q-values
581	// calculating q-values can't be done unless the non-linear corrections are calculated
582	// so go ahead and put it in this loop.
583	// TODO :
584	// -- make sure that everything is present before the calculation
585	// -- beam center must be properly defined in terms of real distance
586	// -- distances/zero location/ etc. must be clearly documented for each detector
587	// ** this assumes that NonLinearCorrection() has been run to generate data_RealDistX and Y
588	// ** this routine Makes the waves QTot, qx, qy, qz in each detector folder.
589	//
590	V_Detector_CalcQVals(fname,detStr,destPath)
591
592	endfor
593
594	//"B" is separate
595	Wave w = V_getDetectorDataW(fname,"B")
596	Wave cal_x = V_getDet_cal_x(fname,"B")
597	Wave cal_y = V_getDet_cal_y(fname,"B")
598	V_NonLinearCorrection_B(fname,w,cal_x,cal_y,"B",destPath)
599	V_ConvertBeamCtr_to_mmB(fname,"B",destPath)
600	V_Detector_CalcQVals(fname,"B",destPath)
601
602	else
603	Print "Non-linear correction NOT DONE"
604	endif
605
606
607	// (3) dead time correction
608	// TODO:
609	// x- remove the hard-wired test - done
610	// -- test for correct operation
611	// x- loop over all of the detectors
612	// x- B detector is a special case (do separately, then loop over NoB)
613	// -- this DOES alter the data
614	// -- verify the error propagation (not done yet)
615	//
616	Variable countRate
617	NVAR gDoDeadTimeCor = root:Packages:NIST:VSANS:Globals:gDoDeadTimeCor
618	if (gDoDeadTimeCor == 1)
619	Print "Doing DeadTime correction"// for "+ detStr
620	for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
621	detStr = StringFromList(ii, ksDetectorListAll, ";")
622	Wave w = V_getDetectorDataW(fname,detStr)
623	Wave w_err = V_getDetectorDataErrW(fname,detStr)
624	ctTime = V_getCount_time(fname)
625
626	if(cmpstr(detStr,"B") == 0)
627	Variable b_dt = V_getDetector_deadtime_B(fname,detStr)
628	// do the correction for the back panel
629	countRate = sum(w,-inf,inf)/ctTime //use sum of detector counts
630
631	w = w/(1-countRate*b_dt)
632	w_err = w_err/(1-countRate*b_dt)
633
634	else
635	// do the corrections for 8 tube panels
636	Wave w_dt = V_getDetector_deadtime(fname,detStr)
637	V_DeadTimeCorrectionTubes(w,w_err,w_dt,ctTime)
638
639	endif
640	endfor
641
642	else
643	Print "Dead Time correction NOT DONE"
644	endif
645
646
647	// (4) solid angle correction
648	// TODO -- this currently calculates the correction factor AND applys it to the data
649	// -- as a result, the data values are very large since they are divided by a very small
650	// solid angle per pixel. But all of the count values are now on the basis of
651	// counts/(solid angle) --- meaning that they can all be binned together for I(q)
652	// -and- TODO - this will need to be taken into account for absolute scaling (this part is already done)
653	//
654	// the solid angle correction is calculated for ALL detector panels.
655	NVAR gDoSolidAngleCor = root:Packages:NIST:VSANS:Globals:gDoSolidAngleCor
656	if (gDoSolidAngleCor == 1)
657	Print "Doing Solid Angle correction"// for "+ detStr
658	for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
659	detStr = StringFromList(ii, ksDetectorListAll, ";")
660	Wave w = V_getDetectorDataW(fname,detStr)
661	Wave w_err = V_getDetectorDataErrW(fname,detStr)
662	// any other dimensions to pass in?
663	V_SolidAngleCorrection(w,w_err,fname,detStr,destPath)
664
665	endfor
666	else
667	Print "Solid Angle correction NOT DONE"
668	endif
669
670
671	// (5) angle-dependent tube shadowing
672	// TODO:
673	// -- not sure about this correction yet...
674	//
675	NVAR gDoTubeShadowCor = root:Packages:NIST:VSANS:Globals:gDoTubeShadowCor
676	if (gDoTubeShadowCor == 1)
677	Print "(stub)Tube Shadow correction"
678	else
679	Print "Tube shadowing correction NOT DONE"
680	endif
681
682	// (6) angle dependent transmission correction
683	// TODO:
684	// x- this still needs to be filled in
685	// -- still some debate of when/where in the corrections that this is best applied
686	// - do it here, and it's done whether the output is 1D or 2D
687	// - do it later (where SAMPLE information is used) since this section is ONLY instrument-specific
688	// -- verify that the calculation is correct
689	// -- verify that the error propagation (in 2D) is correct
690	//
691	NVAR gDoTrans = root:Packages:NIST:VSANS:Globals:gDoTransmissionCor
692	if (gDoTrans == 1)
693	Print "Doing Large-angle transmission correction"// for "+ detStr
694	for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
695	detStr = StringFromList(ii, ksDetectorListAll, ";")
696	Wave w = V_getDetectorDataW(fname,detStr)
697	Wave w_err = V_getDetectorDataErrW(fname,detStr)
698
699	V_LargeAngleTransmissionCorr(w,w_err,fname,detStr,destPath)
700
701	endfor
702	else
703	Print "Sample Transmission correction NOT DONE"
704	endif
705
706	// (7) normalize to default monitor counts
707	// TODO(DONE) x- each detector is rescaled separately, but the rescaling factor is global (only one monitor!)
708	// TODO -- but there are TWO monitors - so how to switch?
709	// --- AND, there is also /entry/control/monitor_counts !!! Which one is the correct value? Which will NICE write
710	// TODO -- what do I really need to save?
711
712	NVAR gDoMonitorNormalization = root:Packages:NIST:VSANS:Globals:gDoMonitorNormalization
713	if (gDoMonitorNormalization == 1)
714
715	Variable monCount,savedMonCount
716	defmon=1e8 //default monitor counts
717	monCount = V_getMonitorCount(fname) // TODO -- this is read in since VCALC fakes this on output
718	// monCount = V_getBeamMonNormData(fname) // TODO -- I think this is the real one to read
719	savedMonCount = monCount
720	scale = defMon/monCount // scale factor to MULTIPLY data by to rescale to defmon
721
722	// write to newType=fname will put new values in the destination WORK folder
723	V_writeBeamMonNormSaved_count(fname,savedMonCount) // save the true count
724	V_writeBeamMonNormData(fname,defMon) // mon ct is now 10^8
725
726	// // TODO
727	// // the low efficiency monitor, expect to use this for white beam mode
728	// -- need code switch here to determine which monitor to use.
729	//
730	// V_getBeamMonLowData(fname)
731	// V_getBeamMonLowSaved_count(fname)
732
733	for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
734	detStr = StringFromList(ii, ksDetectorListAll, ";")
735	Wave w = V_getDetectorDataW(fname,detStr)
736	Wave w_err = V_getDetectorDataErrW(fname,detStr)
737
738	// do the calculation right here. It's a simple scaling and not worth sending to another function.
739	//scale the data and error to the default montor counts
740
741	//
742	w *= scale
743	w_err *= scale //assumes total monitor count is so large there is essentially no error
744
745	// TODO
746	// -- do I want to update and save the integrated detector count?
747	Variable integratedCount = V_getDet_IntegratedCount(fname,detStr)
748	V_writeDet_IntegratedCount(fname,detStr,integratedCount*scale)
749
750	endfor
751	else
752	Print "Monitor Normalization correction NOT DONE"
753	endif
754
755
756	// (not done) angle dependent efficiency correction
757	NVAR doEfficiency = root:Packages:NIST:VSANS:Globals:gDoDetectorEffCor
758
759	//
760	///// these lines are if files are added together, not done (yet) for VSANS
761	//
762	//update totals to put in the work header (at the end of the function)
763	// total_mon += realsread[0]
764	//
765	// total_det += dscale*realsread[2]
766	//
767	// total_trn += realsread[39]
768	// total_rtime += integersread[2]
769	// total_numruns +=1
770	//
771
772	//
773	//reset the current displaytype to "newtype"
774	String/G root:Packages:NIST:VSANS:Globals:gCurDispType=newType
775
776	//return to root folder (redundant)
777	SetDataFolder root:
778
779	Return(0)
780	End
781
782
783	//will "ADD" the current contents of the RAW folder to the newType work folder
784	//and will ADD the RAW contents to the existing content of the newType folder
785	// - used when adding multiple runs together
786	//(the function Raw_to_work(type) makes a fresh workfile)
787	//
788	//the current display type is updated to newType (global)
789	Function V_Add_raw_to_work(newType)
790	String newType
791
792	// NEW OCT 2014
793	// this corrects for adding raw data files with different attenuation
794	// does nothing if the attenuation of RAW and destination are the same
795	NVAR doAdjustRAW_Atten = root:Packages:NIST:gDoAdjustRAW_Atten
796	if(doAdjustRAW_Atten)
797	V_Adjust_RAW_Attenuation(newType)
798	endif
799
800	String destPath=""
801
802	// if the desired workfile doesn't exist, let the user know, and just make a new one
803	if(WaveExists($("root:Packages:NIST:"+newType + ":data")) == 0)
804	Print "There is no old work file to add to - a new one will be created"
805	//call Raw_to_work(), then return from this function
806	V_Raw_to_Work(newType)
807	Return(0) //does not generate an error - a single file was converted to work.newtype
808	Endif
809
810	NVAR pixelsX = root:myGlobals:gNPixelsX
811	NVAR pixelsY = root:myGlobals:gNPixelsY
812
813	//now make references to data in newType folder
814	DestPath="root:Packages:NIST:"+newType
815	WAVE data=$(destPath +":linear_data") // these wave references point to the EXISTING work data
816	WAVE data_copy=$(destPath +":data") // these wave references point to the EXISTING work data
817	WAVE dest_data_err=$(destPath +":linear_data_error") // these wave references point to the EXISTING work data
818	WAVE/T textread=$(destPath + ":textread")
819	WAVE integersread=$(destPath + ":integersread")
820	WAVE realsread=$(destPath + ":realsread")
821
822	Variable deadTime,defmon,total_mon,total_det,total_trn,total_numruns,total_rtime
823	Variable ii,jj,itim,cntrate,dscale,scale,uscale,wrk_beamx,wrk_beamy,xshift,yshift
824
825
826	defmon=1e8 //default monitor counts
827
828	//Yes, add to previous run(s) in work, that does exist
829	//use the actual monitor count run.savmon rather than the normalized monitor count
830	//in run.moncnt and unscale the work data
831
832	total_mon = realsread[1] //saved monitor count
833	uscale = total_mon/defmon //unscaling factor
834	total_det = uscale*realsread[2] //unscaled detector count
835	total_trn = uscale*realsread[39] //unscaled trans det count
836	total_numruns = integersread[3] //number of runs in workfile
837	total_rtime = integersread[2] //total counting time in workfile
838	//retrieve workfile beamcenter
839	wrk_beamx = realsread[16]
840	wrk_beamy = realsread[17]
841	//unscale the workfile data in "newType"
842	//
843	//check for log-scaling and adjust if necessary
844	// should not be needed now - using display flag instead
845	// ConvertFolderToLinearScale(newType)
846	//
847	//then unscale the data array
848	data *= uscale
849	dest_data_err *= uscale
850
851	//DetCorr() has not been applied to the data in RAW , do it now in a local reference to the raw data
852	WAVE raw_data = $"root:Packages:NIST:RAW:linear_data"
853	WAVE raw_data_err = $"root:Packages:NIST:RAW:linear_data_error"
854	WAVE raw_reals = $"root:Packages:NIST:RAW:realsread"
855	WAVE/T raw_text = $"root:Packages:NIST:RAW:textread"
856	WAVE raw_ints = $"root:Packages:NIST:RAW:integersread"
857
858	//check for log-scaling of the raw data - make sure it's linear
859	// should not be needed now - using display flag instead
860	// ConvertFolderToLinearScale("RAW")
861
862	// switches to control what is done, don't do the transmission correction for the BGD measurement
863	NVAR doEfficiency = root:Packages:NIST:gDoDetectorEffCor
864	NVAR gDoTrans = root:Packages:NIST:gDoTransmissionCor
865	Variable doTrans = gDoTrans
866	if(cmpstr("BGD",newtype) == 0)
867	doTrans = 0 //skip the trans correction for the BGD file but don't change the value of the global
868	endif
869
870	V_DetCorr(raw_data,raw_data_err,raw_reals,doEfficiency,doTrans) //applies correction to raw_data, and overwrites it
871
872	//deadtime corrections to raw data
873	// TODO - do the tube correction for dead time now
874	deadTime = 1//DetectorDeadtime(raw_text[3],raw_text[9],dateAndTimeStr=raw_text[1],dtime=raw_reals[48]) //pick the correct detector deadtime, switch on date too
875	itim = raw_ints[2]
876	cntrate = sum(raw_data,-inf,inf)/itim //080802 use data sum, rather than scaler value
877	dscale = 1/(1-deadTime*cntrate)
878
879	#if (exists("ILL_D22")==6)
880	Variable tubeSum
881	// for D22 detector might need to use cntrate/128 as it is the tube response
882	for(ii=0;ii<pixelsX;ii+=1)
883	//sum the counts in each tube
884	tubeSum = 0
885	for(jj=0;jj<pixelsY;jj+=1)
886	tubeSum += data[jj][ii]
887	endfor
888	// countrate in tube ii
889	cntrate = tubeSum/itim
890	// deadtime scaling in tube ii
891	dscale = 1/(1-deadTime*cntrate)
892	// multiply data[ii][] by the dead time
893	raw_data[][ii] *= dscale
894	raw_data_err[][ii] *= dscale
895	endfor
896	#else
897	// dead time correction on all other RAW data, including NCNR
898	raw_data *= dscale
899	raw_data_err *= dscale
900	#endif
901
902	//update totals by adding RAW values to the local ones (write to work header at end of function)
903	total_mon += raw_reals[0]
904
905	total_det += dscale*raw_reals[2]
906
907	total_trn += raw_reals[39]
908	total_rtime += raw_ints[2]
909	total_numruns +=1
910
911	//do the beamcenter shifting if there is a mismatch
912	//and then add the two data sets together, changing "data" since it is the workfile data
913	xshift = raw_reals[16] - wrk_beamx
914	yshift = raw_reals[17] - wrk_beamy
915
916	If((xshift != 0) \|\| (yshift != 0))
917	DoAlert 1,"Do you want to ignore the beam center mismatch?"
918	if(V_flag==1)
919	xshift=0
920	yshift=0
921	endif
922	endif
923
924	If((xshift == 0) && (yshift == 0)) //no shift, just add them
925	data += raw_data //deadtime correction has already been done to the raw data
926	dest_data_err = sqrt(dest_data_err^2 + raw_data_err^2) // error of the sum
927	Endif
928
929	//scale the data to the default montor counts
930	scale = defmon/total_mon
931	data *= scale
932	dest_data_err *= scale
933
934	// keep "data" and linear_data in sync in the destination folder
935	data_copy = data
936
937	//all is done, except for the bookkeeping of updating the header info in the work folder
938	textread[1] = date() + " " + time() //date + time stamp
939	integersread[3] = total_numruns //numruns = more than one
940	realsread[1] = total_mon //save the true monitor count
941	realsread[0] = defmon //monitor ct = defmon
942	integersread[2] = total_rtime // total counting time
943	realsread[2] = scale*total_det //scaled detector counts
944	realsread[39] = scale*total_trn //scaled transmission counts
945
946	//Add the added raw filename to the list of files in the workfile
947	String newfile = ";" + raw_text[0]
948	SVAR oldList = $(destPath + ":fileList")
949	String/G $(destPath + ":fileList") = oldList + newfile
950
951	//reset the current display type to "newtype"
952	SVAR gCurDispType = root:Packages:NIST:VSANS:Globals:gCurDispType
953	gCurDispType = newType
954
955	//return to root folder (redundant)
956	SetDataFolder root:
957
958	Return(0)
959	End
960
961
962	//used for adding DRK (beam shutter CLOSED) data to a workfile
963	//force the monitor count to 1, since it's irrelevant
964	// run data through normal "add" step, then unscale default monitor counts
965	//to get the data back on a simple time basis
966	//
967	Function V_Raw_to_Work_NoNorm(type)
968	String type
969
970	WAVE reals=$("root:Packages:NIST:RAW:realsread")
971	reals[1]=1 //true monitor counts, still in raw
972	V_Raw_to_work(type)
973	//data is now in "type" folder
974	WAVE data=$("root:Packages:NIST:"+type+":linear_data")
975	WAVE data_copy=$("root:Packages:NIST:"+type+":data")
976	WAVE data_err=$("root:Packages:NIST:"+type+":linear_data_error")
977	WAVE new_reals=$("root:Packages:NIST:"+type+":realsread")
978
979	Variable norm_mon,tot_mon,scale
980
981	norm_mon = new_reals[0] //should be 1e8
982	tot_mon = new_reals[1] //should be 1
983	scale= norm_mon/tot_mon
984
985	data /= scale //unscale the data
986	data_err /= scale
987
988	// to keep "data" and linear_data in sync
989	data_copy = data
990
991	return(0)
992	End
993
994	//used for adding DRK (beam shutter CLOSED) data to a workfile
995	//force the monitor count to 1, since it's irrelevant
996	// run data through normal "add" step, then unscale default monitor counts
997	//to get the data back on a simple time basis
998	//
999	Function V_Add_Raw_to_Work_NoNorm(type)
1000	String type
1001
1002	WAVE reals=$("root:Packages:NIST:RAW:realsread")
1003	reals[1]=1 //true monitor counts, still in raw
1004	V_Add_Raw_to_work(type)
1005	//data is now in "type" folder
1006	WAVE data=$("root:Packages:NIST:"+type+":linear_data")
1007	WAVE data_copy=$("root:Packages:NIST:"+type+":data")
1008	WAVE data_err=$("root:Packages:NIST:"+type+":linear_data_error")
1009	WAVE new_reals=$("root:Packages:NIST:"+type+":realsread")
1010
1011	Variable norm_mon,tot_mon,scale
1012
1013	norm_mon = new_reals[0] //should be 1e8
1014	tot_mon = new_reals[1] //should be equal to the number of runs (1 count per run)
1015	scale= norm_mon/tot_mon
1016
1017	data /= scale //unscale the data
1018	data_err /= scale
1019
1020	// to keep "data" and linear_data in sync
1021	data_copy = data
1022
1023	return(0)
1024	End
1025

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

Download in other formats:

Original Format