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

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

some addtions:

DIV file generation is now aware of the High resolution detector, but the procedures are still awaiting data for testing.

Read Noise file can now be read in and stored in the (RAW) folder ReadNoise?. This is not a work folder and the data isnot changed from the RAW state. This image is then subtracted from other raw data as it is converted to a work file (SAM, EMP, etc.) Previously, only a constant value was subtracted. If the ReadNoise? data is not present, the constant will be subtracted. There is a menu option to load the ReadNoise? data.

File size: 37.4 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// May 2019
473// If present, a detector image of the back detector containing the read noise (non-uniform values)
474//  is subtracted rather than a single constant value.
475//
476//the current display type is updated to newType (global)
477//
478Function V_Raw_to_work(newType)
479        String newType
480       
481        Variable deadTime,defmon,total_mon,total_det,total_trn,total_numruns,total_rtime
482        Variable ii,jj,itim,cntrate,dscale,scale,uscale
483        String destPath
484       
485        String fname = newType
486        String detStr
487        Variable ctTime
488
489        //initialize values before normalization
490        total_mon=0
491        total_det=0
492        total_trn=0
493        total_numruns=0
494        total_rtime=0
495       
496        //Not adding multiple runs, so wipe out the old contents of the work folder and
497        // replace with the contents of raw
498
499        destPath = "root:Packages:NIST:VSANS:" + newType
500       
501        //copy from current dir (RAW) to work, defined by newType
502        V_CopyHDFToWorkFolder("RAW",newType)
503       
504        // now work with the waves from the destination folder.
505       
506        // apply corrections ---
507        // switches to control what is done, don't do the transmission correction for the BGD measurement
508        // start with the DIV correction, before conversion to mm
509        // then do all of the other corrections, order doesn't matter.
510        // rescaling to default monitor counts however, must be LAST.
511
512// each correction must loop over each detector. tedious.
513
514        //except for removing the read noise of the back detector
515        NVAR gIgnoreDetB = root:Packages:NIST:VSANS:Globals:gIgnoreDetB
516
517        if(gIgnoreDetB == 0)
518                Wave w = V_getDetectorDataW(fname,"B")
519                // I hate to hard-wire this, but the data must be in this specific location...
520                Wave/Z w_ReadNoise = $("root:Packages:NIST:VSANS:ReadNoise:entry:instrument:detector_B:data")
521//              Wave/Z w_ReadNoise = V_getDetectorDataW("ReadNoise","B")
522               
523                if(WaveExists(w_ReadNoise))
524                        w -= w_ReadNoise
525                        Print "Subtracting ReadNoise Array"
526                else
527                        NVAR gHighResBinning = root:Packages:NIST:VSANS:Globals:gHighResBinning
528                       
529                        switch(gHighResBinning)
530                                case 1:
531                                        w -= kReadNoiseLevel_bin1               // a constant value
532                                       
533        //                              MatrixFilter /N=11 /P=1 median w                        //              /P=n flag sets the number of passes (default is 1 pass)                 
534        //                              Print "*** median noise filter 11x11 applied to the back detector (1 pass) ***"
535                                        Print "*** 1x1 binning - subtracted ReadNoise Constant - No Filter ***"
536                                        break
537                                case 4:
538                                        w -= kReadNoiseLevel_bin4               // a constant value
539                                       
540        //                              MatrixFilter /N=3 /P=3 median w                 //              /P=n flag sets the number of passes (default is 1 pass)                         
541        //                              Print "*** median noise filter 3x3 applied to the back detector (3 passes) ***"
542                                        Print "*** 4x4 binning - subtracted ReadNoise Constant - No Filter ***"
543                                        break
544                                default:
545                                        Abort "No binning case matches in V_Raw_to_Work"
546                        endswitch       
547                endif           //waveExists
548        endif           // using det B
549       
550       
551        // (0) Redimension the data waves in the destination folder
552        //     so that they are DP, not integer
553        // TODO
554        // -- currently only redimensioning the data and linear_data_error - anything else???
555        // -- ?? some of this is done at load time for RAW data. Not a problem to re-do the redimension
556        for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
557                detStr = StringFromList(ii, ksDetectorListAll, ";")
558                Wave w = V_getDetectorDataW(fname,detStr)
559                Wave w_err = V_getDetectorDataErrW(fname,detStr)
560                Redimension/D w,w_err
561        endfor
562       
563       
564        // (1) DIV correction
565        // do this in terms of pixels.
566        // TODO : This must also exist at the time the first work folder is generated.
567        //   So it must be in the user folder at the start of the experiment, and defined.
568        NVAR gDoDIVCor = root:Packages:NIST:VSANS:Globals:gDoDIVCor
569        if (gDoDIVCor == 1)
570                // need extra check here for file existence
571                // if not in DIV folder, load.
572                // if unable to load, skip correction and report error (Alert?) (Ask to Load?)
573                Print "Doing DIV correction"// for "+ detStr
574                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
575                        detStr = StringFromList(ii, ksDetectorListAll, ";")
576                       
577                        if(cmpstr(detStr,"B") == 0  && gIgnoreDetB == 1)
578                                // do nothing
579                        else
580                                Wave w = V_getDetectorDataW(fname,detStr)
581                                Wave w_err = V_getDetectorDataErrW(fname,detStr)
582                               
583                                V_DIVCorrection(w,w_err,detStr,newType)
584                        endif
585                endfor
586        else
587                Print "DIV correction NOT DONE"         // not an error since correction was unchecked
588        endif
589       
590        // (2) non-linear correction   
591        // TODO:
592        // x-  the "B" detector is calculated in its own routines
593        // -- document what is generated here:
594        //    **in each detector folder: data_realDistX and data_realDistY (2D waves of the [mm] position of each pixel)
595        // x- these spatial calculations ARE DONE as the RAW data is loaded. It allows the RAW
596        //    data to be properly displayed, but without all of the (complete) set of detector corrections
597        // * the corrected distances are calculated into arrays, but nothing is done with them yet
598        // * there is enough information now to calculate the q-arrays, so it is done now
599        // - other corrections may modify the data, this calculation does NOT modify the data
600        NVAR gDoNonLinearCor = root:Packages:NIST:VSANS:Globals:gDoNonLinearCor
601        // generate a distance matrix for each of the detectors
602        if (gDoNonLinearCor == 1)
603                Print "Doing Non-linear correction"// for "+ detStr
604                for(ii=0;ii<ItemsInList(ksDetectorListNoB);ii+=1)
605                        detStr = StringFromList(ii, ksDetectorListNoB, ";")
606                        Wave w = V_getDetectorDataW(fname,detStr)
607//                      Wave w_err = V_getDetectorDataErrW(fname,detStr)
608                        Wave w_calib = V_getDetTube_spatialCalib(fname,detStr)
609                        Variable tube_width = V_getDet_tubeWidth(fname,detStr)
610                        V_NonLinearCorrection(fname,w,w_calib,tube_width,detStr,destPath)
611
612                        //(2.4) Convert the beam center values from pixels to mm
613                        // TODO -- there needs to be a permanent location for these values??
614                        //
615                                // TODO
616                                // -- the beam center value in mm needs to be present - it is used in calculation of Qvalues
617                                // -- but having both the same is wrong...
618                                // -- the pixel value is needed for display of the panels
619                                if(kBCTR_CM)
620                                        //V_ConvertBeamCtrPix_to_mm(folder,detStr,destPath)
621                                        //
622       
623                                        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
624                                        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
625                                        WAVE x_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
626                                        WAVE y_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
627                                        x_mm[0] = V_getDet_beam_center_x(fname,detStr) * 10             // convert cm to mm
628                                        y_mm[0] = V_getDet_beam_center_y(fname,detStr) * 10             // convert cm to mm
629                                       
630                                        // TODO:::
631                                // now I need to convert the beam center in mm to pixels
632                                // and have some rational place to look for it...
633                                        V_ConvertBeamCtr_to_pix(fname,detStr,destPath)
634                                else
635                                        // beam center is in pixels, so use the old routine
636                                        V_ConvertBeamCtrPix_to_mm(fname,detStr,destPath)
637                                endif           
638                                                       
639                        // (2.5) Calculate the q-values
640                        // calculating q-values can't be done unless the non-linear corrections are calculated
641                        // so go ahead and put it in this loop.
642                        // TODO :
643                        // -- make sure that everything is present before the calculation
644                        // -- beam center must be properly defined in terms of real distance
645                        // -- distances/zero location/ etc. must be clearly documented for each detector
646                        //      ** this assumes that NonLinearCorrection() has been run to generate data_RealDistX and Y
647                        // ** this routine Makes the waves QTot, qx, qy, qz in each detector folder.
648                        //
649                        V_Detector_CalcQVals(fname,detStr,destPath)
650                       
651                endfor
652
653                if(gIgnoreDetB==1)             
654                        //"B" is separate
655                        detStr = "B"
656                        Wave w = V_getDetectorDataW(fname,detStr)
657                        Wave cal_x = V_getDet_cal_x(fname,detStr)
658                        Wave cal_y = V_getDet_cal_y(fname,detStr)
659                       
660                        V_NonLinearCorrection_B(fname,w,cal_x,cal_y,detStr,destPath)
661       
662        // "B" is always naturally defined in terms of a pixel center. This can be converted to mm,
663        // but the experiment will measure pixel x,y - just like ordela detectors.
664                       
665        //              if(kBCTR_CM)
666        //
667        //                      Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
668        //                      Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
669        //                      WAVE x_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
670        //                      WAVE y_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
671        //                      x_mm[0] = V_getDet_beam_center_x(fname,detStr) * 10             // convert cm to mm
672        //                      y_mm[0] = V_getDet_beam_center_y(fname,detStr) * 10             // convert cm to mm
673        //                     
674        //                      // now I need to convert the beam center in mm to pixels
675        //                      // and have some rational place to look for it...
676        //                      V_ConvertBeamCtr_to_pixB(fname,detStr,destPath)
677        //              else
678                                // beam center is in pixels, so use the old routine
679                                V_ConvertBeamCtrPix_to_mmB(fname,detStr,destPath)
680       
681        //              endif
682                        V_Detector_CalcQVals(fname,detStr,destPath)
683               
684                endif
685               
686        else
687                Print "Non-linear correction NOT DONE"
688        endif
689
690
691        // (3) dead time correction
692        // DONE:
693        // x- test for correct operation
694        // x- loop over all of the detectors
695        // x- B detector is a special case (do separately, then loop over NoB)
696        // x- this DOES alter the data
697        // x- verify the error propagation (not done yet)
698        //
699        Variable countRate
700        NVAR gDoDeadTimeCor = root:Packages:NIST:VSANS:Globals:gDoDeadTimeCor
701        if (gDoDeadTimeCor == 1)
702                Print "Doing DeadTime correction"// for "+ detStr
703                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
704                        detStr = StringFromList(ii, ksDetectorListAll, ";")
705                        Wave w = V_getDetectorDataW(fname,detStr)
706                        Wave w_err = V_getDetectorDataErrW(fname,detStr)
707                        ctTime = V_getCount_time(fname)
708
709                        if(cmpstr(detStr,"B") == 0 )
710                                if(gIgnoreDetB == 0)
711                                        Variable b_dt = V_getDetector_deadtime_B(fname,detStr)
712                                        // do the correction for the back panel
713                                        countRate = sum(w,-inf,inf)/ctTime              //use sum of detector counts
714       
715                                        w = w/(1-countRate*b_dt)
716                                        w_err = w_err/(1-countRate*b_dt)
717                                endif                   
718                        else
719                                // do the corrections for 8 tube panels
720                                Wave w_dt = V_getDetector_deadtime(fname,detStr)
721                                V_DeadTimeCorrectionTubes(w,w_err,w_dt,ctTime)
722                               
723                        endif
724                endfor
725               
726        else
727                Print "Dead Time correction NOT DONE"
728        endif   
729       
730
731        // (4) solid angle correction
732        //  -- this currently calculates the correction factor AND applys it to the data
733        //  -- as a result, the data values are very large since they are divided by a very small
734        //     solid angle per pixel. But all of the count values are now on the basis of
735        //    counts/(solid angle) --- meaning that they can all be binned together for I(q)
736        //    -and- - this is taken into account for absolute scaling (this part is already done)
737        //
738        // the solid angle correction is calculated for ALL detector panels.
739        NVAR gDoSolidAngleCor = root:Packages:NIST:VSANS:Globals:gDoSolidAngleCor
740        if (gDoSolidAngleCor == 1)
741                Print "Doing Solid Angle correction"// for "+ detStr
742                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
743                        detStr = StringFromList(ii, ksDetectorListAll, ";")
744                       
745                        if(cmpstr(detStr,"B") == 0  && gIgnoreDetB == 1)
746                                // do nothing
747                        else
748                                Wave w = V_getDetectorDataW(fname,detStr)
749                                Wave w_err = V_getDetectorDataErrW(fname,detStr)
750                                // any other dimensions to pass in?
751                                V_SolidAngleCorrection(w,w_err,fname,detStr,destPath)
752                        endif
753                       
754                endfor
755        else
756                Print "Solid Angle correction NOT DONE"
757        endif   
758       
759               
760        // (5) angle-dependent tube shadowing
761        // TODO:
762        // -- not sure about this correction yet...
763        //
764        NVAR gDoTubeShadowCor = root:Packages:NIST:VSANS:Globals:gDoTubeShadowCor
765        if (gDoTubeShadowCor == 1)
766                Print "(stub)Tube Shadow correction"
767        else
768                Print "Tube shadowing correction NOT DONE"
769        endif   
770               
771        // (6) angle dependent transmission correction
772        // TODO:
773        // x- this still needs to be filled in
774        // -- still some debate of when/where in the corrections that this is best applied
775        //    - do it here, and it's done whether the output is 1D or 2D
776        //    - do it later (where SAMPLE information is used) since this section is ONLY instrument-specific
777        // -- verify that the calculation is correct
778        // -- verify that the error propagation (in 2D) is correct
779        //
780        NVAR gDoTrans = root:Packages:NIST:VSANS:Globals:gDoTransmissionCor
781        if (gDoTrans == 1)
782                Print "Doing Large-angle transmission correction"// for "+ detStr
783                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
784                        detStr = StringFromList(ii, ksDetectorListAll, ";")
785                       
786                        if(cmpstr(detStr,"B") == 0  && gIgnoreDetB == 1)
787                                // do nothing
788                        else
789                                Wave w = V_getDetectorDataW(fname,detStr)
790                                Wave w_err = V_getDetectorDataErrW(fname,detStr)
791                               
792                                V_LargeAngleTransmissionCorr(w,w_err,fname,detStr,destPath)
793                        endif
794                endfor
795        else
796                Print "Sample Transmission correction NOT DONE"
797        endif   
798       
799        // (7) normalize to default monitor counts
800        // TODO(DONE) x- each detector is rescaled separately, but the rescaling factor is global (only one monitor!)
801        // TODO -- but there are TWO monitors - so how to switch?
802        //  --- AND, there is also /entry/control/monitor_counts !!! Which one is the correct value? Which will NICE write
803        // TODO -- what do I really need to save?
804       
805        NVAR gDoMonitorNormalization = root:Packages:NIST:VSANS:Globals:gDoMonitorNormalization
806        if (gDoMonitorNormalization == 1)
807               
808                Variable monCount,savedMonCount
809                defmon=1e8                      //default monitor counts
810//              monCount = V_getControlMonitorCount(fname)                      // TODO -- this is read in since VCALC fakes this on output
811                monCount = V_getBeamMonNormData(fname)          // TODO -- I think this is the *real* one to read
812                savedMonCount   = monCount
813                scale = defMon/monCount         // scale factor to MULTIPLY data by to rescale to defmon
814
815                // write to newType=fname will put new values in the destination WORK folder
816                V_writeBeamMonNormSaved_count(fname,savedMonCount)                      // save the true count
817                V_writeBeamMonNormData(fname,defMon)            // mon ct is now 10^8
818                                       
819//                      // TODO
820//                      // the low efficiency monitor, expect to use this for white beam mode
821//                              -- need code switch here to determine which monitor to use.
822//
823//                      V_getBeamMonLowData(fname)
824//                      V_getBeamMonLowSaved_count(fname)       
825                       
826                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
827                        detStr = StringFromList(ii, ksDetectorListAll, ";")
828                        Wave w = V_getDetectorDataW(fname,detStr)
829                        Wave w_err = V_getDetectorDataErrW(fname,detStr)
830
831                        // do the calculation right here. It's a simple scaling and not worth sending to another function.     
832                        //scale the data and error to the default monitor counts
833               
834//
835                        w *= scale
836                        w_err *= scale          //assumes total monitor count is so large there is essentially no error
837
838                        // TODO
839                        // -- do I want to update and save the integrated detector count?
840                        Variable integratedCount = V_getDet_IntegratedCount(fname,detStr)
841                        V_writeDet_IntegratedCount(fname,detStr,integratedCount*scale)
842
843                endfor
844        else
845                Print "Monitor Normalization correction NOT DONE"
846        endif
847       
848       
849        // (not done) angle dependent efficiency correction
850        NVAR doEfficiency = root:Packages:NIST:VSANS:Globals:gDoDetectorEffCor
851
852//     
853        //reset the current displaytype to "newtype"
854        String/G root:Packages:NIST:VSANS:Globals:gCurDispType=newType
855       
856        //return to root folder (redundant)
857        SetDataFolder root:
858       
859        Return(0)
860End
861
862//
863//will "ADD" the current contents of the RAW folder to the newType work folder
864//
865// - used when adding multiple runs together
866//(the function V_Raw_to_work(type) makes a fresh workfile)
867//
868//the current display type is updated to newType (global)
869Function V_Add_raw_to_work(newType)
870        String newType
871       
872        String destPath=""
873        // if the desired workfile doesn't exist, let the user know, and just make a new one
874        if(WaveExists($("root:Packages:NIST:VSANS:"+newType + ":entry:instrument:detector_FL:data")) == 0)
875                Print "There is no old work file to add to - a new one will be created"
876                //call V_Raw_to_work(), then return from this function
877                V_Raw_to_Work(newType)
878                Return(0)               //does not generate an error - a single file was converted to work.newtype
879        Endif
880
881
882        // convert the RAW data to a WORK file.
883        // this will do all of the necessary corrections to the data
884        // put this in some separate work folder that can be cleaned out at the end (ADJ)
885        String tmpType="ADJ"
886       
887        //this step removes the read noise from the back so that neither added file will have this constant
888        V_Raw_to_Work(tmpType)         
889                       
890        //now make references to data in newType folder
891        DestPath="root:Packages:NIST:VSANS:"+newType   
892
893
894/////////////////
895//fields that need to be added together
896// entry block
897        // collection_time              V_getCollectionTime(fname)              V_putCollectionTime(fname,val)
898
899// instrument block
900        // beam_monitor_norm
901                // data (this will be 1e8)                              V_getBeamMonNormData(fname)             V_putBeamMonNormData(fname,val)
902                // saved_count (this is the original monitor count)  V_getBeamMonNormSaved_count(fname)         V_putBeamMonNormSaved_count(fname,val)
903
904        // for each detector
905        // data         V_getDetectorDataW(fname,detStr)
906        // integrated_count             V_getDet_IntegratedCount(fname,detStr)   V_putDet_IntegratedCount(fname,detStr,val)
907        // linear_data           V_getDetectorLinearDataW(fname,detStr)
908        // RECALCULATE (or add properly) linear_data_error              V_getDetectorDataErrW(fname,detStr)
909
910
911// control block (these may not actually be used?)
912        // count_time                           V_getCount_time(fname)                                          V_putCount_time(fname,val)
913        // detector_counts              V_getDetector_counts(fname)                             V_putDetector_counts(fname,val)
914        // monitor_counts               V_getControlMonitorCount(fname)         V_putControlMonitorCount(fname,val)
915
916// sample block - nothing
917// reduction block - nothing
918// user block - nothing
919
920// ?? need to add the file name to a list of what was actually added - so it will be saved with I(q)
921//
922////////////////////
923
924//      total_mon = realsread[1]        //saved monitor count
925//      uscale = total_mon/defmon               //unscaling factor
926//      total_det = uscale*realsread[2]         //unscaled detector count
927//      total_trn = uscale*realsread[39]        //unscaled trans det count
928//      total_numruns = integersread[3] //number of runs in workfile
929//      total_rtime = integersread[2]           //total counting time in workfile
930       
931
932        String detStr
933       
934        Variable saved_mon_dest,scale_dest,saved_mon_tmp,scale_tmp
935        Variable collection_time_dest,collection_time_tmp,count_time_dest,count_time_tmp
936        Variable detCount_dest,detCount_tmp,det_integrated_ct_dest,det_integrated_ct_tmp
937        Variable ii,new_scale,defMon
938       
939        defMon=1e8                      //default monitor counts
940
941        // find the scaling factors, one for each folder
942        saved_mon_dest = V_getBeamMonNormSaved_count(newType)
943        scale_dest = saved_mon_dest/defMon              //un-scaling factor
944       
945        saved_mon_tmp = V_getBeamMonNormSaved_count(tmpType)
946        scale_tmp = saved_mon_tmp/defMon                        //un-scaling factor
947
948        new_scale = defMon / (saved_mon_dest+saved_mon_tmp)
949       
950       
951        // get the count time for each (two locations)
952        collection_time_dest = V_getCollectionTime(newType)
953        collection_time_tmp = V_getCollectionTime(tmpType)
954       
955        count_time_dest = V_getCount_time(newType)
956        count_time_tmp = V_getCount_time(tmpType)
957       
958        detCount_dest = V_getDetector_counts(newType)
959        detCount_tmp = V_getDetector_counts(tmpType)
960
961// update the fields that are not in the detector blocks
962// in entry
963        V_putCollectionTime(newType,collection_time_dest+collection_time_tmp)
964
965// in control block
966        V_putCount_time(newType,count_time_dest+count_time_tmp)
967        V_putDetector_counts(newType,detCount_dest+detCount_tmp)
968        V_putControlMonitorCount(newType,saved_mon_dest+saved_mon_tmp)
969
970
971// now loop over all of the detector panels
972        // data
973        // data_err
974        // integrated count
975        // linear_data
976        for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
977                detStr = StringFromList(ii, ksDetectorListAll, ";")
978               
979                Wave data_dest = V_getDetectorDataW(newType,detStr)
980                Wave data_err_dest = V_getDetectorDataErrW(newType,detStr)
981                Wave linear_data_dest = V_getDetectorLinearDataW(newType,detStr)
982                det_integrated_ct_dest = V_getDet_IntegratedCount(newType,detStr)
983
984                Wave data_tmp = V_getDetectorDataW(tmpType,detStr)
985                Wave data_err_tmp = V_getDetectorDataErrW(tmpType,detStr)
986                Wave linear_data_tmp = V_getDetectorLinearDataW(tmpType,detStr)
987                det_integrated_ct_tmp = V_getDet_IntegratedCount(tmpType,detStr)
988               
989                // unscale the data arrays
990                data_dest *= scale_dest
991                data_err_dest *= scale_dest
992                linear_data_dest *= scale_dest
993               
994                data_tmp *= scale_tmp
995                data_err_tmp *= scale_tmp
996                linear_data_tmp *= scale_tmp
997                               
998                // add them together, the dest is a wave so it is automatically changed in the "dest" folder
999                V_putDet_IntegratedCount(tmpType,detStr,det_integrated_ct_dest+det_integrated_ct_tmp)
1000                data_dest += data_tmp
1001                data_err_dest = sqrt(data_err_dest^2 + data_err_tmp^2)          // add in quadrature
1002                linear_data_dest += linear_data_tmp
1003               
1004                // now rescale the data_dest to the monitor counts
1005                data_dest *= new_scale
1006                data_err_dest *= new_scale
1007                linear_data_dest *= new_scale
1008               
1009        endfor
1010
1011       
1012        //Add the added raw filename to the list of files in the workfile
1013        SVAR fileList_dest = $("root:Packages:NIST:VSANS:"+newType+":gFileList")
1014        SVAR fileList_tmp = $("root:Packages:NIST:VSANS:"+tmpType+":gFileList")
1015       
1016        fileList_dest += ";" + fileList_tmp
1017       
1018        //reset the current display type to "newtype"
1019        SVAR gCurDispType = root:Packages:NIST:VSANS:Globals:gCurDispType
1020        gCurDispType = newType
1021       
1022        //return to root folder (redundant)
1023        SetDataFolder root:
1024       
1025        Return(0)
1026End
1027
1028
1029//used for adding DRK (beam shutter CLOSED) data to a workfile
1030//force the monitor count to 1, since it's irrelevant
1031// run data through normal "add" step, then unscale default monitor counts
1032//to get the data back on a simple time basis
1033//
1034Function V_Raw_to_Work_NoNorm(type)
1035        String type
1036       
1037        WAVE reals=$("root:Packages:NIST:RAW:realsread")
1038        reals[1]=1              //true monitor counts, still in raw
1039        V_Raw_to_work(type)
1040        //data is now in "type" folder
1041        WAVE data=$("root:Packages:NIST:"+type+":linear_data")
1042        WAVE data_copy=$("root:Packages:NIST:"+type+":data")
1043        WAVE data_err=$("root:Packages:NIST:"+type+":linear_data_error")
1044        WAVE new_reals=$("root:Packages:NIST:"+type+":realsread")
1045       
1046        Variable norm_mon,tot_mon,scale
1047       
1048        norm_mon = new_reals[0]         //should be 1e8
1049        tot_mon = new_reals[1]          //should be 1
1050        scale= norm_mon/tot_mon
1051       
1052        data /= scale           //unscale the data
1053        data_err /= scale
1054       
1055        // to keep "data" and linear_data in sync
1056        data_copy = data
1057       
1058        return(0)
1059End
1060
1061//used for adding DRK (beam shutter CLOSED) data to a workfile
1062//force the monitor count to 1, since it's irrelevant
1063// run data through normal "add" step, then unscale default monitor counts
1064//to get the data back on a simple time basis
1065//
1066Function V_Add_Raw_to_Work_NoNorm(type)
1067        String type
1068       
1069        WAVE reals=$("root:Packages:NIST:RAW:realsread")
1070        reals[1]=1              //true monitor counts, still in raw
1071        V_Add_Raw_to_work(type)
1072        //data is now in "type" folder
1073        WAVE data=$("root:Packages:NIST:"+type+":linear_data")
1074        WAVE data_copy=$("root:Packages:NIST:"+type+":data")
1075        WAVE data_err=$("root:Packages:NIST:"+type+":linear_data_error")
1076        WAVE new_reals=$("root:Packages:NIST:"+type+":realsread")
1077       
1078        Variable norm_mon,tot_mon,scale
1079       
1080        norm_mon = new_reals[0]         //should be 1e8
1081        tot_mon = new_reals[1]          //should be equal to the number of runs (1 count per run)
1082        scale= norm_mon/tot_mon
1083       
1084        data /= scale           //unscale the data
1085        data_err /= scale
1086       
1087        // to keep "data" and linear_data in sync
1088        data_copy = data
1089       
1090        return(0)
1091End
1092
Note: See TracBrowser for help on using the repository browser.