source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_WorkFolderUtils.ipf @ 1145

Last change on this file since 1145 was 1145, checked in by srkline, 3 years ago

added a line in to handle data folders that might be written to HDF with a space in the name. This is OK in HDF5, but not handled well in Igor. This was found in the temperature_env field.

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