Context Navigation

V_WorkFolderUtils.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: 36.2 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	// gFileList is above "entry" which is my additions
111	SVAR fileList_dest = $("root:Packages:NIST:VSANS:"+toStr+":gFileList")
112	SVAR fileList_tmp = $("root:Packages:NIST:VSANS:"+fromStr+":gFileList")
113	fileList_dest = fileList_tmp
114
115	// everything on the top level
116	V_DuplicateDataFolder($(fromDF+":entry"),fromStr,toStr,0,"",0) //no recursion here
117	// control
118	V_DuplicateDataFolder($(fromDF+":entry:control"),fromStr,toStr,0,"",1) //yes recursion here
119	// instrument
120	V_DuplicateDataFolder($(fromDF+":entry:instrument"),fromStr,toStr,0,"",1) //yes recursion here
121	// reduction
122	V_DuplicateDataFolder($(fromDF+":entry:reduction"),fromStr,toStr,0,"",1) //yes recursion here
123	// sample
124	V_DuplicateDataFolder($(fromDF+":entry:sample"),fromStr,toStr,0,"",1) //yes recursion here
125	// user
126	V_DuplicateDataFolder($(fromDF+":entry:user"),fromStr,toStr,0,"",1) //yes recursion here
127
128	endif
129
130	return(0)
131	end
132
133
134
135	Function V_KillWavesInFolder(folderStr)
136	String folderStr
137
138	if(DataFolderExists(folderStr) && strlen(folderStr) != 0)
139	SetDataFolder $folderStr
140	KillWaves/A/Z
141	endif
142
143	SetDataFolder root:
144	return(0)
145	end
146
147	////////
148	// see the help entry for IndexedDir for help on (possibly) how to do this faster
149	// -- see Function ScanDirectories(pathName, printDirNames)
150	//
151
152
153	// from IgorExchange On July 17th, 2011 jtigor
154	// started from "Recursively List Data Folder Contents"
155	// Posted July 15th, 2011 by hrodstein
156	//
157	//
158	//
159	Proc V_CopyWorkFolderProc(dataFolderStr, fromStr, toStr, level, sNBName, recurse)
160	String dataFolderStr="root:Packages:NIST:VSANS:RAW"
161	String fromStr = "RAW"
162	String toStr="SAM"
163	Variable level=0
164	String sNBName="DataFolderTree"
165	Variable recurse = 1
166
167	V_DuplicateDataFolder($dataFolderStr, fromStr, toStr, level, sNBName, recurse)
168
169
170	end
171
172	// ListDataFolder(dfr, level)
173	// Recursively lists objects in data folder.
174	// Pass data folder path for dfr and 0 for level.
175	// pass level == 0 for the first call
176	// sNBName = "" prints nothing. any name will generate a notebook
177	//
178	// recurse == 0 will do only the specified folder, anything else will recurse all levels
179	// toStr is the string name of the top-level folder only, not the full path
180	//
181	//
182	Function V_DuplicateDataFolder(dfr, fromStr, toStr, level, sNBName,recurse)
183	DFREF dfr
184	String fromStr
185	String toStr
186	Variable level // Pass 0 to start
187	String sNBName
188	Variable recurse
189
190	String name
191	String dfName
192	String sString
193
194	String toDF = ""
195
196	if (level == 0) // this is the data folder, generate if needed in the destination
197	name = GetDataFolder(1, dfr)
198	// sPrintf sString, "%s (data folder)\r", name
199	toDF = ReplaceString(fromStr,name,toStr,1) // case-sensitive replace
200	sprintf sString, "NewDataFolder/O %s\r",toDF
201	NewDataFolder/O $(RemoveEnding(toDF,":")) // remove trailing semicolon if it's there
202
203	V_WriteBrowserInfo_test(sString, 1, sNBName)
204	endif
205
206	dfName = GetDataFolder(1, dfr)
207	toDF = ReplaceString(fromStr,dfName,toStr,1) // case-sensitive replace
208	Variable i
209
210	String indentStr = "\t"
211	for(i=0; i<level; i+=1)
212	indentStr += "\t"
213	endfor
214
215	Variable numWaves = CountObjectsDFR(dfr, 1)
216	for(i=0; i<numWaves; i+=1)
217	name = GetIndexedObjNameDFR(dfr, 1, i)
218	//
219	// wave type does not matter now. Duplicate does not care
220	//
221	sPrintf sString, "Duplicate/O %s, %s\r",dfName+name,toDF+name
222	Duplicate/O $(dfName+name),$(toDF+name)
223
224	V_WriteBrowserInfo_test(sString, 2, sNBName)
225	endfor
226
227	Variable numNumericVariables = CountObjectsDFR(dfr, 2)
228	for(i=0; i<numNumericVariables; i+=1)
229	name = GetIndexedObjNameDFR(dfr, 2, i)
230	sPrintf sString, "%s%s (numeric variable)\r", indentStr, name
231	V_WriteBrowserInfo_test(sString, 3, sNBName)
232	endfor
233
234	Variable numStringVariables = CountObjectsDFR(dfr, 3)
235	for(i=0; i<numStringVariables; i+=1)
236	name = GetIndexedObjNameDFR(dfr, 3, i)
237	sPrintf sString, "%s%s (string variable)\r", indentStr, name
238	V_WriteBrowserInfo_test(sString, 4, sNBName)
239	endfor
240
241	if(recurse)
242	Variable numDataFolders = CountObjectsDFR(dfr, 4)
243	for(i=0; i<numDataFolders; i+=1)
244	name = GetIndexedObjNameDFR(dfr, 4, i)
245	// sPrintf sString, "%s%s (data folder)\r", indentStr, name
246	dfName = GetDataFolder(1, dfr)
247
248	toDF = ReplaceString(fromStr,dfName,toStr,1) // case-sensitive replace
249	sprintf sString, "NewDataFolder/O %s\r",toDF+name
250	NewDataFolder/O $(toDF+name)
251
252
253	V_WriteBrowserInfo_test(sString, 1, sNBName)
254	DFREF childDFR = dfr:$(name)
255	V_DuplicateDataFolder(childDFR, fromStr, toStr, level+1, sNBName, recurse)
256	endfor
257	endif
258
259
260	End
261
262
263	// ListDataFolder(dfr, level)
264	// Recursively lists objects in data folder.
265	// Pass data folder path for dfr and 0 for level.
266	// pass level == 0 for the first call
267	// sNBName = "" prints nothing. any name will generate a notebook
268	//
269	// recurse == 0 will do only the specified folder, anything else will recurse all levels
270	// toStr is the string name of the top-level folder only, not the full path
271	//
272	//
273	Function V_KillWavesFullTree(dfr, fromStr, level, sNBName,recurse)
274	DFREF dfr
275	String fromStr
276	// String toStr
277	Variable level // Pass 0 to start
278	String sNBName
279	Variable recurse
280
281	String name
282	String dfName
283	String sString
284
285	String toDF = ""
286
287	if (level == 0) // this is the data folder, generate if needed in the destination
288	name = GetDataFolder(1, dfr)
289	sPrintf sString, "%s (data folder)\r", name
290	// toDF = ReplaceString(fromStr,name,toStr,1) // case-sensitive replace
291	// sprintf sString, "NewDataFolder/O %s\r",toDF
292	// NewDataFolder/O $(RemoveEnding(toDF,":")) // remove trailing semicolon if it's there
293
294	V_WriteBrowserInfo_test(sString, 1, sNBName)
295	endif
296
297	dfName = GetDataFolder(1, dfr)
298	// toDF = ReplaceString(fromStr,dfName,toStr,1) // case-sensitive replace
299	Variable i
300
301	String indentStr = "\t"
302	for(i=0; i<level; i+=1)
303	indentStr += "\t"
304	endfor
305
306	Variable numWaves = CountObjectsDFR(dfr, 1)
307	for(i=0; i<numWaves; i+=1)
308	name = GetIndexedObjNameDFR(dfr, 1, i)
309	//
310	// wave type does not matter now. Kill does not care
311	//
312	sPrintf sString, "Killing %s\r",dfName+name
313	KillWaves/Z $(dfName+name)
314
315	V_WriteBrowserInfo_test(sString, 2, sNBName)
316	endfor
317
318	// now kill the data folder if possible
319	KillDataFolder/Z $dfName
320
321
322	Variable numNumericVariables = CountObjectsDFR(dfr, 2)
323	for(i=0; i<numNumericVariables; i+=1)
324	name = GetIndexedObjNameDFR(dfr, 2, i)
325	sPrintf sString, "%s%s (numeric variable)\r", indentStr, name
326	V_WriteBrowserInfo_test(sString, 3, sNBName)
327	endfor
328
329	Variable numStringVariables = CountObjectsDFR(dfr, 3)
330	for(i=0; i<numStringVariables; i+=1)
331	name = GetIndexedObjNameDFR(dfr, 3, i)
332	sPrintf sString, "%s%s (string variable)\r", indentStr, name
333	V_WriteBrowserInfo_test(sString, 4, sNBName)
334	endfor
335
336	if(recurse)
337	Variable numDataFolders = CountObjectsDFR(dfr, 4)
338	for(i=0; i<numDataFolders; i+=1)
339	name = GetIndexedObjNameDFR(dfr, 4, i)
340	sPrintf sString, "%s%s (data folder)\r", indentStr, name
341	dfName = GetDataFolder(1, dfr)
342
343	// toDF = ReplaceString(fromStr,dfName,toStr,1) // case-sensitive replace
344	// sprintf sString, "NewDataFolder/O %s\r",toDF+name
345	// NewDataFolder/O $(toDF+name)
346
347
348	V_WriteBrowserInfo_test(sString, 1, sNBName)
349	DFREF childDFR = dfr:$(name)
350	V_KillWavesFullTree(childDFR, fromStr, level+1, sNBName, recurse)
351	endfor
352	endif
353
354
355	End
356
357
358
359	Function V_WriteBrowserInfo_test(sString, vType, sNBName)
360	String sString
361	Variable vType
362	String sNBName
363
364	if(strlen(sNBName) == 0)
365	// print sString
366	return 0
367	endif
368	DoWindow $sNBName
369	if(V_flag != 1)
370	NewNoteBook/F=0 /N=$sNBName /V=1 as sNBName
371	else
372	DoWindow/F $sNBName
373	endif
374	Notebook $sNBName selection={endOfFile, endOfFile}
375	if(vType == 1) // a data folder
376	// Notebook $sNBName fstyle=1
377	Notebook $sNBName text=sString
378	// Notebook $sNBName fstyle=-1
379	else
380	Notebook $sNBName text=sString
381	endif
382
383	End
384
385	///////////////////////////////
386
387
388	//
389	// given the folder, duplicate the data -> linear_data and generate the error
390	//
391	// x- do I want to use different names here? If it turns out that I don't need to drag a copy of
392	// the data around as "linear_data", then I can eliminate that, and rename the error wave
393	// x- be sure the data is either properly written as 2D in the file, or converted to 2D before
394	// duplicating here
395	// x- ? do I recast to DP here? No- V_MakeDataWaves_DP() is called directly prior to this call, so data
396	// coming in is already DP
397	Function V_MakeDataError(folderStr)
398	String folderStr
399
400	SetDataFolder $folderStr
401	Wave data=data
402	Duplicate/O data linear_data // at this point, the data is still the raw data, and is linear_data
403
404	// proper error for counting statistics, good for low count values too
405	// rather than just sqrt(n)
406	// see N. Gehrels, Astrophys. J., 303 (1986) 336-346, equation (7)
407	// for S = 1 in eq (7), this corresponds to one sigma error bars
408	Duplicate/O linear_data linear_data_error
409	linear_data_error = 1 + sqrt(linear_data + 0.75)
410
411	Duplicate/O linear_data_error data_error
412	//
413
414	SetDataFolder root:
415	return(0)
416	End
417
418
419	/////////////////////
420
421
422
423	//testing procedure
424	// TODO -- can't duplicate this with another proceudre, but if I change the name of the variable
425	// "newType" to "type", then when Raw_to_work() gets to CopyHDFToWorkFolder(), the KillDataFolder/Z
426	// line fails (but reports no error), then DuplicateDataFolder fails, and reports an error. Trying
427	// to simplify this condition, I can't reproduce the bug for WM...
428	Proc V_Convert_to_Workfile(newtype, doadd)
429	String newtype,doadd
430	Prompt newtype,"WORK data type",popup,"SAM;EMP;BGD;ADJ;"
431	Prompt doadd,"Add to current WORK contents?",popup,"No;Yes;"
432
433	//macro will take whatever is in RAW folder and "ADD" it to the folder specified
434	//in the popup menu
435
436	//"add" = yes/no, don't add to previous runs
437	//switch here - two separate functions to avoid (my) confusion
438	Variable err// = Raw_to_work(newtype)
439	if(cmpstr(doadd,"No")==0)
440	//don't add to prev work contents, copy RAW contents to work and convert
441	err = V_Raw_to_work(newtype)
442	else
443	//yes, add RAW to the current work folder contents
444	//Abort "Adding RAW data files is currently unsupported"
445	err = V_Add_raw_to_work(newtype)
446	endif
447
448	String newTitle = "WORK_"+newtype
449	DoWindow/F VSANS_Data
450	DoWindow/T VSANS_Data, newTitle
451	KillStrings/Z newTitle
452
453	//need to update the display with "data" from the correct dataFolder
454	V_UpdateDisplayInformation(newtype)
455
456	End
457
458
459	//
460	// THIS IS THE MAJOR ROUTINE TO APPLY DATA CORRECTIONS
461	//
462	//will copy the current contents of the RAW folder to the newType work folder
463	//and do the geometric corrections and normalization to monitor counts
464	//(the function Add_Raw_to_work(type) adds multiple runs together - and is LOW priority)
465	//
466	// JUL 2018
467	// now removes a constant kReadNoiseLevel from the back detector (this is a constant value, not time
468	// or count dependent). It does not appear to be dependent on gain, and is hoepfully stable over time.
469	//
470	//the current display type is updated to newType (global)
471	//
472	Function V_Raw_to_work(newType)
473	String newType
474
475	Variable deadTime,defmon,total_mon,total_det,total_trn,total_numruns,total_rtime
476	Variable ii,jj,itim,cntrate,dscale,scale,uscale
477	String destPath
478
479	String fname = newType
480	String detStr
481	Variable ctTime
482
483	//initialize values before normalization
484	total_mon=0
485	total_det=0
486	total_trn=0
487	total_numruns=0
488	total_rtime=0
489
490	//Not adding multiple runs, so wipe out the old contents of the work folder and
491	// replace with the contents of raw
492
493	destPath = "root:Packages:NIST:VSANS:" + newType
494
495	//copy from current dir (RAW) to work, defined by newType
496	V_CopyHDFToWorkFolder("RAW",newType)
497
498	// now work with the waves from the destination folder.
499
500	// apply corrections ---
501	// switches to control what is done, don't do the transmission correction for the BGD measurement
502	// start with the DIV correction, before conversion to mm
503	// then do all of the other corrections, order doesn't matter.
504	// rescaling to default monitor counts however, must be LAST.
505
506	// each correction must loop over each detector. tedious.
507
508	//except for removing the read noise of the back detector
509	NVAR gIgnoreDetB = root:Packages:NIST:VSANS:Globals:gIgnoreDetB
510
511	if(gIgnoreDetB == 0)
512	Wave w = V_getDetectorDataW(fname,"B")
513
514	NVAR gHighResBinning = root:Packages:NIST:VSANS:Globals:gHighResBinning
515
516	switch(gHighResBinning)
517	case 1:
518	w -= kReadNoiseLevel_bin1 // a constant value
519
520	MatrixFilter /N=11 /P=1 median w // /P=n flag sets the number of passes (default is 1 pass)
521
522	Print "* median noise filter 11x11 applied to the back detector (1 pass) *"
523	break
524	case 4:
525	w -= kReadNoiseLevel_bin4 // a constant value
526
527	MatrixFilter /N=3 /P=3 median w // /P=n flag sets the number of passes (default is 1 pass)
528
529	Print "* median noise filter 3x3 applied to the back detector (3 passes) *"
530	break
531	default:
532	Abort "No binning case matches in V_Raw_to_Work"
533	endswitch
534
535	endif
536
537
538	// (0) Redimension the data waves in the destination folder
539	// so that they are DP, not integer
540	// TODO
541	// -- currently only redimensioning the data and linear_data_error - anything else???
542	// -- ?? some of this is done at load time for RAW data. Not a problem to re-do the redimension
543	for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
544	detStr = StringFromList(ii, ksDetectorListAll, ";")
545	Wave w = V_getDetectorDataW(fname,detStr)
546	Wave w_err = V_getDetectorDataErrW(fname,detStr)
547	Redimension/D w,w_err
548	endfor
549
550
551	// (1) DIV correction
552	// do this in terms of pixels.
553	// TODO : This must also exist at the time the first work folder is generated.
554	// So it must be in the user folder at the start of the experiment, and defined.
555	NVAR gDoDIVCor = root:Packages:NIST:VSANS:Globals:gDoDIVCor
556	if (gDoDIVCor == 1)
557	// need extra check here for file existence
558	// if not in DIV folder, load.
559	// if unable to load, skip correction and report error (Alert?) (Ask to Load?)
560	Print "Doing DIV correction"// for "+ detStr
561	for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
562	detStr = StringFromList(ii, ksDetectorListAll, ";")
563	Wave w = V_getDetectorDataW(fname,detStr)
564	Wave w_err = V_getDetectorDataErrW(fname,detStr)
565
566	V_DIVCorrection(w,w_err,detStr,newType)
567	endfor
568	else
569	Print "DIV correction NOT DONE" // not an error since correction was unchecked
570	endif
571
572	// (2) non-linear correction
573	// TODO:
574	// x- the "B" detector is calculated in its own routines
575	// -- document what is generated here:
576	// **in each detector folder: data_realDistX and data_realDistY (2D waves of the [mm] position of each pixel)
577	// x- these spatial calculations ARE DONE as the RAW data is loaded. It allows the RAW
578	// data to be properly displayed, but without all of the (complete) set of detector corrections
579	// * the corrected distances are calculated into arrays, but nothing is done with them yet
580	// * there is enough information now to calculate the q-arrays, so it is done now
581	// - other corrections may modify the data, this calculation does NOT modify the data
582	NVAR gDoNonLinearCor = root:Packages:NIST:VSANS:Globals:gDoNonLinearCor
583	// generate a distance matrix for each of the detectors
584	if (gDoNonLinearCor == 1)
585	Print "Doing Non-linear correction"// for "+ detStr
586	for(ii=0;ii<ItemsInList(ksDetectorListNoB);ii+=1)
587	detStr = StringFromList(ii, ksDetectorListNoB, ";")
588	Wave w = V_getDetectorDataW(fname,detStr)
589	// Wave w_err = V_getDetectorDataErrW(fname,detStr)
590	Wave w_calib = V_getDetTube_spatialCalib(fname,detStr)
591	Variable tube_width = V_getDet_tubeWidth(fname,detStr)
592	V_NonLinearCorrection(fname,w,w_calib,tube_width,detStr,destPath)
593
594	//(2.4) Convert the beam center values from pixels to mm
595	// TODO -- there needs to be a permanent location for these values??
596	//
597	// TODO
598	// -- the beam center value in mm needs to be present - it is used in calculation of Qvalues
599	// -- but having both the same is wrong...
600	// -- the pixel value is needed for display of the panels
601	if(kBCTR_CM)
602	//V_ConvertBeamCtrPix_to_mm(folder,detStr,destPath)
603	//
604
605	Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
606	Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
607	WAVE x_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
608	WAVE y_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
609	x_mm[0] = V_getDet_beam_center_x(fname,detStr) * 10 // convert cm to mm
610	y_mm[0] = V_getDet_beam_center_y(fname,detStr) * 10 // convert cm to mm
611
612	// TODO:::
613	// now I need to convert the beam center in mm to pixels
614	// and have some rational place to look for it...
615	V_ConvertBeamCtr_to_pix(fname,detStr,destPath)
616	else
617	// beam center is in pixels, so use the old routine
618	V_ConvertBeamCtrPix_to_mm(fname,detStr,destPath)
619	endif
620
621	// (2.5) Calculate the q-values
622	// calculating q-values can't be done unless the non-linear corrections are calculated
623	// so go ahead and put it in this loop.
624	// TODO :
625	// -- make sure that everything is present before the calculation
626	// -- beam center must be properly defined in terms of real distance
627	// -- distances/zero location/ etc. must be clearly documented for each detector
628	// ** this assumes that NonLinearCorrection() has been run to generate data_RealDistX and Y
629	// ** this routine Makes the waves QTot, qx, qy, qz in each detector folder.
630	//
631	V_Detector_CalcQVals(fname,detStr,destPath)
632
633	endfor
634
635	//"B" is separate
636	detStr = "B"
637	Wave w = V_getDetectorDataW(fname,detStr)
638	Wave cal_x = V_getDet_cal_x(fname,detStr)
639	Wave cal_y = V_getDet_cal_y(fname,detStr)
640
641	V_NonLinearCorrection_B(fname,w,cal_x,cal_y,detStr,destPath)
642
643	// "B" is always naturally defined in terms of a pixel center. This can be converted to mm,
644	// but the experiment will measure pixel x,y - just like ordela detectors.
645
646	// if(kBCTR_CM)
647	//
648	// Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
649	// Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
650	// WAVE x_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
651	// WAVE y_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
652	// x_mm[0] = V_getDet_beam_center_x(fname,detStr) * 10 // convert cm to mm
653	// y_mm[0] = V_getDet_beam_center_y(fname,detStr) * 10 // convert cm to mm
654	//
655	// // now I need to convert the beam center in mm to pixels
656	// // and have some rational place to look for it...
657	// V_ConvertBeamCtr_to_pixB(fname,detStr,destPath)
658	// else
659	// beam center is in pixels, so use the old routine
660	V_ConvertBeamCtrPix_to_mmB(fname,detStr,destPath)
661
662	// endif
663	V_Detector_CalcQVals(fname,detStr,destPath)
664
665	else
666	Print "Non-linear correction NOT DONE"
667	endif
668
669
670	// (3) dead time correction
671	// DONE:
672	// x- test for correct operation
673	// x- loop over all of the detectors
674	// x- B detector is a special case (do separately, then loop over NoB)
675	// x- this DOES alter the data
676	// x- verify the error propagation (not done yet)
677	//
678	Variable countRate
679	NVAR gDoDeadTimeCor = root:Packages:NIST:VSANS:Globals:gDoDeadTimeCor
680	if (gDoDeadTimeCor == 1)
681	Print "Doing DeadTime correction"// for "+ detStr
682	for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
683	detStr = StringFromList(ii, ksDetectorListAll, ";")
684	Wave w = V_getDetectorDataW(fname,detStr)
685	Wave w_err = V_getDetectorDataErrW(fname,detStr)
686	ctTime = V_getCount_time(fname)
687
688	if(cmpstr(detStr,"B") == 0)
689	Variable b_dt = V_getDetector_deadtime_B(fname,detStr)
690	// do the correction for the back panel
691	countRate = sum(w,-inf,inf)/ctTime //use sum of detector counts
692
693	w = w/(1-countRate*b_dt)
694	w_err = w_err/(1-countRate*b_dt)
695
696	else
697	// do the corrections for 8 tube panels
698	Wave w_dt = V_getDetector_deadtime(fname,detStr)
699	V_DeadTimeCorrectionTubes(w,w_err,w_dt,ctTime)
700
701	endif
702	endfor
703
704	else
705	Print "Dead Time correction NOT DONE"
706	endif
707
708
709	// (4) solid angle correction
710	// -- this currently calculates the correction factor AND applys it to the data
711	// -- as a result, the data values are very large since they are divided by a very small
712	// solid angle per pixel. But all of the count values are now on the basis of
713	// counts/(solid angle) --- meaning that they can all be binned together for I(q)
714	// -and- - this is taken into account for absolute scaling (this part is already done)
715	//
716	// the solid angle correction is calculated for ALL detector panels.
717	NVAR gDoSolidAngleCor = root:Packages:NIST:VSANS:Globals:gDoSolidAngleCor
718	if (gDoSolidAngleCor == 1)
719	Print "Doing Solid Angle correction"// for "+ detStr
720	for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
721	detStr = StringFromList(ii, ksDetectorListAll, ";")
722	Wave w = V_getDetectorDataW(fname,detStr)
723	Wave w_err = V_getDetectorDataErrW(fname,detStr)
724	// any other dimensions to pass in?
725	V_SolidAngleCorrection(w,w_err,fname,detStr,destPath)
726
727	endfor
728	else
729	Print "Solid Angle correction NOT DONE"
730	endif
731
732
733	// (5) angle-dependent tube shadowing
734	// TODO:
735	// -- not sure about this correction yet...
736	//
737	NVAR gDoTubeShadowCor = root:Packages:NIST:VSANS:Globals:gDoTubeShadowCor
738	if (gDoTubeShadowCor == 1)
739	Print "(stub)Tube Shadow correction"
740	else
741	Print "Tube shadowing correction NOT DONE"
742	endif
743
744	// (6) angle dependent transmission correction
745	// TODO:
746	// x- this still needs to be filled in
747	// -- still some debate of when/where in the corrections that this is best applied
748	// - do it here, and it's done whether the output is 1D or 2D
749	// - do it later (where SAMPLE information is used) since this section is ONLY instrument-specific
750	// -- verify that the calculation is correct
751	// -- verify that the error propagation (in 2D) is correct
752	//
753	NVAR gDoTrans = root:Packages:NIST:VSANS:Globals:gDoTransmissionCor
754	if (gDoTrans == 1)
755	Print "Doing Large-angle transmission correction"// for "+ detStr
756	for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
757	detStr = StringFromList(ii, ksDetectorListAll, ";")
758	Wave w = V_getDetectorDataW(fname,detStr)
759	Wave w_err = V_getDetectorDataErrW(fname,detStr)
760
761	V_LargeAngleTransmissionCorr(w,w_err,fname,detStr,destPath)
762
763	endfor
764	else
765	Print "Sample Transmission correction NOT DONE"
766	endif
767
768	// (7) normalize to default monitor counts
769	// TODO(DONE) x- each detector is rescaled separately, but the rescaling factor is global (only one monitor!)
770	// TODO -- but there are TWO monitors - so how to switch?
771	// --- AND, there is also /entry/control/monitor_counts !!! Which one is the correct value? Which will NICE write
772	// TODO -- what do I really need to save?
773
774	NVAR gDoMonitorNormalization = root:Packages:NIST:VSANS:Globals:gDoMonitorNormalization
775	if (gDoMonitorNormalization == 1)
776
777	Variable monCount,savedMonCount
778	defmon=1e8 //default monitor counts
779	// monCount = V_getControlMonitorCount(fname) // TODO -- this is read in since VCALC fakes this on output
780	monCount = V_getBeamMonNormData(fname) // TODO -- I think this is the real one to read
781	savedMonCount = monCount
782	scale = defMon/monCount // scale factor to MULTIPLY data by to rescale to defmon
783
784	// write to newType=fname will put new values in the destination WORK folder
785	V_writeBeamMonNormSaved_count(fname,savedMonCount) // save the true count
786	V_writeBeamMonNormData(fname,defMon) // mon ct is now 10^8
787
788	// // TODO
789	// // the low efficiency monitor, expect to use this for white beam mode
790	// -- need code switch here to determine which monitor to use.
791	//
792	// V_getBeamMonLowData(fname)
793	// V_getBeamMonLowSaved_count(fname)
794
795	for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
796	detStr = StringFromList(ii, ksDetectorListAll, ";")
797	Wave w = V_getDetectorDataW(fname,detStr)
798	Wave w_err = V_getDetectorDataErrW(fname,detStr)
799
800	// do the calculation right here. It's a simple scaling and not worth sending to another function.
801	//scale the data and error to the default monitor counts
802
803	//
804	w *= scale
805	w_err *= scale //assumes total monitor count is so large there is essentially no error
806
807	// TODO
808	// -- do I want to update and save the integrated detector count?
809	Variable integratedCount = V_getDet_IntegratedCount(fname,detStr)
810	V_writeDet_IntegratedCount(fname,detStr,integratedCount*scale)
811
812	endfor
813	else
814	Print "Monitor Normalization correction NOT DONE"
815	endif
816
817
818	// (not done) angle dependent efficiency correction
819	NVAR doEfficiency = root:Packages:NIST:VSANS:Globals:gDoDetectorEffCor
820
821	//
822	//reset the current displaytype to "newtype"
823	String/G root:Packages:NIST:VSANS:Globals:gCurDispType=newType
824
825	//return to root folder (redundant)
826	SetDataFolder root:
827
828	Return(0)
829	End
830
831	//
832	//will "ADD" the current contents of the RAW folder to the newType work folder
833	//
834	// - used when adding multiple runs together
835	//(the function V_Raw_to_work(type) makes a fresh workfile)
836	//
837	//the current display type is updated to newType (global)
838	Function V_Add_raw_to_work(newType)
839	String newType
840
841	String destPath=""
842	// if the desired workfile doesn't exist, let the user know, and just make a new one
843	if(WaveExists($("root:Packages:NIST:VSANS:"+newType + ":entry:instrument:detector_FL:data")) == 0)
844	Print "There is no old work file to add to - a new one will be created"
845	//call V_Raw_to_work(), then return from this function
846	V_Raw_to_Work(newType)
847	Return(0) //does not generate an error - a single file was converted to work.newtype
848	Endif
849
850
851	// convert the RAW data to a WORK file.
852	// this will do all of the necessary corrections to the data
853	// put this in some separate work folder that can be cleaned out at the end (ADJ)
854	String tmpType="ADJ"
855
856	//this step removes the read noise from the back so that neither added file will have this constant
857	V_Raw_to_Work(tmpType)
858
859	//now make references to data in newType folder
860	DestPath="root:Packages:NIST:VSANS:"+newType
861
862
863	/////////////////
864	//fields that need to be added together
865	// entry block
866	// collection_time V_getCollectionTime(fname) V_putCollectionTime(fname,val)
867
868	// instrument block
869	// beam_monitor_norm
870	// data (this will be 1e8) V_getBeamMonNormData(fname) V_putBeamMonNormData(fname,val)
871	// saved_count (this is the original monitor count) V_getBeamMonNormSaved_count(fname) V_putBeamMonNormSaved_count(fname,val)
872
873	// for each detector
874	// data V_getDetectorDataW(fname,detStr)
875	// integrated_count V_getDet_IntegratedCount(fname,detStr) V_putDet_IntegratedCount(fname,detStr,val)
876	// linear_data V_getDetectorLinearDataW(fname,detStr)
877	// RECALCULATE (or add properly) linear_data_error V_getDetectorDataErrW(fname,detStr)
878
879
880	// control block (these may not actually be used?)
881	// count_time V_getCount_time(fname) V_putCount_time(fname,val)
882	// detector_counts V_getDetector_counts(fname) V_putDetector_counts(fname,val)
883	// monitor_counts V_getControlMonitorCount(fname) V_putControlMonitorCount(fname,val)
884
885	// sample block - nothing
886	// reduction block - nothing
887	// user block - nothing
888
889	// ?? need to add the file name to a list of what was actually added - so it will be saved with I(q)
890	//
891	////////////////////
892
893	// total_mon = realsread[1] //saved monitor count
894	// uscale = total_mon/defmon //unscaling factor
895	// total_det = uscale*realsread[2] //unscaled detector count
896	// total_trn = uscale*realsread[39] //unscaled trans det count
897	// total_numruns = integersread[3] //number of runs in workfile
898	// total_rtime = integersread[2] //total counting time in workfile
899
900
901	String detStr
902
903	Variable saved_mon_dest,scale_dest,saved_mon_tmp,scale_tmp
904	Variable collection_time_dest,collection_time_tmp,count_time_dest,count_time_tmp
905	Variable detCount_dest,detCount_tmp,det_integrated_ct_dest,det_integrated_ct_tmp
906	Variable ii,new_scale,defMon
907
908	defMon=1e8 //default monitor counts
909
910	// find the scaling factors, one for each folder
911	saved_mon_dest = V_getBeamMonNormSaved_count(newType)
912	scale_dest = saved_mon_dest/defMon //un-scaling factor
913
914	saved_mon_tmp = V_getBeamMonNormSaved_count(tmpType)
915	scale_tmp = saved_mon_tmp/defMon //un-scaling factor
916
917	new_scale = defMon / (saved_mon_dest+saved_mon_tmp)
918
919
920	// get the count time for each (two locations)
921	collection_time_dest = V_getCollectionTime(newType)
922	collection_time_tmp = V_getCollectionTime(tmpType)
923
924	count_time_dest = V_getCount_time(newType)
925	count_time_tmp = V_getCount_time(tmpType)
926
927	detCount_dest = V_getDetector_counts(newType)
928	detCount_tmp = V_getDetector_counts(tmpType)
929
930	// update the fields that are not in the detector blocks
931	// in entry
932	V_putCollectionTime(newType,collection_time_dest+collection_time_tmp)
933
934	// in control block
935	V_putCount_time(newType,count_time_dest+count_time_tmp)
936	V_putDetector_counts(newType,detCount_dest+detCount_tmp)
937	V_putControlMonitorCount(newType,saved_mon_dest+saved_mon_tmp)
938
939
940	// now loop over all of the detector panels
941	// data
942	// data_err
943	// integrated count
944	// linear_data
945	for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
946	detStr = StringFromList(ii, ksDetectorListAll, ";")
947
948	Wave data_dest = V_getDetectorDataW(newType,detStr)
949	Wave data_err_dest = V_getDetectorDataErrW(newType,detStr)
950	Wave linear_data_dest = V_getDetectorLinearDataW(newType,detStr)
951	det_integrated_ct_dest = V_getDet_IntegratedCount(newType,detStr)
952
953	Wave data_tmp = V_getDetectorDataW(tmpType,detStr)
954	Wave data_err_tmp = V_getDetectorDataErrW(tmpType,detStr)
955	Wave linear_data_tmp = V_getDetectorLinearDataW(tmpType,detStr)
956	det_integrated_ct_tmp = V_getDet_IntegratedCount(tmpType,detStr)
957
958	// unscale the data arrays
959	data_dest *= scale_dest
960	data_err_dest *= scale_dest
961	linear_data_dest *= scale_dest
962
963	data_tmp *= scale_tmp
964	data_err_tmp *= scale_tmp
965	linear_data_tmp *= scale_tmp
966
967	// add them together, the dest is a wave so it is automatically changed in the "dest" folder
968	V_putDet_IntegratedCount(tmpType,detStr,det_integrated_ct_dest+det_integrated_ct_tmp)
969	data_dest += data_tmp
970	data_err_dest = sqrt(data_err_dest^2 + data_err_tmp^2) // add in quadrature
971	linear_data_dest += linear_data_tmp
972
973	// now rescale the data_dest to the monitor counts
974	data_dest *= new_scale
975	data_err_dest *= new_scale
976	linear_data_dest *= new_scale
977
978	endfor
979
980
981	//Add the added raw filename to the list of files in the workfile
982	SVAR fileList_dest = $("root:Packages:NIST:VSANS:"+newType+":gFileList")
983	SVAR fileList_tmp = $("root:Packages:NIST:VSANS:"+tmpType+":gFileList")
984
985	fileList_dest += ";" + fileList_tmp
986
987	//reset the current display type to "newtype"
988	SVAR gCurDispType = root:Packages:NIST:VSANS:Globals:gCurDispType
989	gCurDispType = newType
990
991	//return to root folder (redundant)
992	SetDataFolder root:
993
994	Return(0)
995	End
996
997
998	//used for adding DRK (beam shutter CLOSED) data to a workfile
999	//force the monitor count to 1, since it's irrelevant
1000	// run data through normal "add" step, then unscale default monitor counts
1001	//to get the data back on a simple time basis
1002	//
1003	Function V_Raw_to_Work_NoNorm(type)
1004	String type
1005
1006	WAVE reals=$("root:Packages:NIST:RAW:realsread")
1007	reals[1]=1 //true monitor counts, still in raw
1008	V_Raw_to_work(type)
1009	//data is now in "type" folder
1010	WAVE data=$("root:Packages:NIST:"+type+":linear_data")
1011	WAVE data_copy=$("root:Packages:NIST:"+type+":data")
1012	WAVE data_err=$("root:Packages:NIST:"+type+":linear_data_error")
1013	WAVE new_reals=$("root:Packages:NIST:"+type+":realsread")
1014
1015	Variable norm_mon,tot_mon,scale
1016
1017	norm_mon = new_reals[0] //should be 1e8
1018	tot_mon = new_reals[1] //should be 1
1019	scale= norm_mon/tot_mon
1020
1021	data /= scale //unscale the data
1022	data_err /= scale
1023
1024	// to keep "data" and linear_data in sync
1025	data_copy = data
1026
1027	return(0)
1028	End
1029
1030	//used for adding DRK (beam shutter CLOSED) data to a workfile
1031	//force the monitor count to 1, since it's irrelevant
1032	// run data through normal "add" step, then unscale default monitor counts
1033	//to get the data back on a simple time basis
1034	//
1035	Function V_Add_Raw_to_Work_NoNorm(type)
1036	String type
1037
1038	WAVE reals=$("root:Packages:NIST:RAW:realsread")
1039	reals[1]=1 //true monitor counts, still in raw
1040	V_Add_Raw_to_work(type)
1041	//data is now in "type" folder
1042	WAVE data=$("root:Packages:NIST:"+type+":linear_data")
1043	WAVE data_copy=$("root:Packages:NIST:"+type+":data")
1044	WAVE data_err=$("root:Packages:NIST:"+type+":linear_data_error")
1045	WAVE new_reals=$("root:Packages:NIST:"+type+":realsread")
1046
1047	Variable norm_mon,tot_mon,scale
1048
1049	norm_mon = new_reals[0] //should be 1e8
1050	tot_mon = new_reals[1] //should be equal to the number of runs (1 count per run)
1051	scale= norm_mon/tot_mon
1052
1053	data /= scale //unscale the data
1054	data_err /= scale
1055
1056	// to keep "data" and linear_data in sync
1057	data_copy = data
1058
1059	return(0)
1060	End
1061

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

Download in other formats:

Original Format