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

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

Changed the calculation of solid angle to properly reflect the tube geometry. This mode is now the (correct) default calculation. A global switch can be applied for testing to use the old method (with many warning alerts).

File size: 39.3 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==0)             
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        NVAR/Z gDo_OLD_SolidAngleCor = root:Packages:NIST:VSANS:Globals:gDo_OLD_SolidAngleCor
741        // for older experiments, this won't exist, so generate it and default to zero
742        // so the old calculation is not done
743        if(NVAR_Exists(gDo_OLD_SolidAngleCor)==0)
744                Variable/G root:Packages:NIST:VSANS:Globals:gDo_OLD_SolidAngleCor=0
745        endif
746        if (gDoSolidAngleCor == 1)
747                Print "Doing Solid Angle correction"// for "+ detStr
748                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
749                        detStr = StringFromList(ii, ksDetectorListAll, ";")
750                       
751                        if(cmpstr(detStr,"B") == 0  && gIgnoreDetB == 1)
752                                // do nothing if the back is ignored
753                        else
754                                Wave w = V_getDetectorDataW(fname,detStr)
755                                Wave w_err = V_getDetectorDataErrW(fname,detStr)
756                                // any other dimensions to pass in?
757                               
758                                if(gDo_OLD_SolidAngleCor == 0)
759                                        V_SolidAngleCorrection(w,w_err,fname,detStr,destPath)
760                                else
761                                        // for testing ONLY -- the cos^3 correction is incorrect for tubes, and the normal
762                                        // function call above   correctly handles either high-res grid or tubes. This COS3 function
763                                        // will incorrectly treat tubes as a grid       
764                                        //                              Print "TESTING -- using incorrect COS^3 solid angle !"         
765                                        V_SolidAngleCorrection_COS3(w,w_err,fname,detStr,destPath)
766                                endif
767                               
768                        endif
769                       
770                endfor
771                if(gDo_OLD_SolidAngleCor == 1)
772                        DoAlert 0,"TESTING -- using incorrect COS^3 solid angle !"             
773                endif
774
775        else
776                Print "Solid Angle correction NOT DONE"
777        endif   
778       
779               
780        // (5) angle-dependent tube shadowing
781        // TODO:
782        // -- not sure about this correction yet...
783        //
784        NVAR gDoTubeShadowCor = root:Packages:NIST:VSANS:Globals:gDoTubeShadowCor
785        if (gDoTubeShadowCor == 1)
786                Print "(stub)Tube Shadow correction"
787        else
788                Print "Tube shadowing correction NOT DONE"
789        endif   
790               
791        // (6) angle dependent transmission correction
792        // TODO:
793        // x- this still needs to be filled in
794        // -- still some debate of when/where in the corrections that this is best applied
795        //    - do it here, and it's done whether the output is 1D or 2D
796        //    - do it later (where SAMPLE information is used) since this section is ONLY instrument-specific
797        // -- verify that the calculation is correct
798        // -- verify that the error propagation (in 2D) is correct
799        //
800        NVAR gDoTrans = root:Packages:NIST:VSANS:Globals:gDoTransmissionCor
801        if (gDoTrans == 1)
802                Print "Doing Large-angle transmission correction"// for "+ detStr
803                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
804                        detStr = StringFromList(ii, ksDetectorListAll, ";")
805                       
806                        if(cmpstr(detStr,"B") == 0  && gIgnoreDetB == 1)
807                                // do nothing
808                        else
809                                Wave w = V_getDetectorDataW(fname,detStr)
810                                Wave w_err = V_getDetectorDataErrW(fname,detStr)
811                               
812                                V_LargeAngleTransmissionCorr(w,w_err,fname,detStr,destPath)
813                        endif
814                endfor
815        else
816                Print "Sample Transmission correction NOT DONE"
817        endif   
818       
819        // (7) normalize to default monitor counts
820        // TODO(DONE) x- each detector is rescaled separately, but the rescaling factor is global (only one monitor!)
821        // TODO -- but there are TWO monitors - so how to switch?
822        //  --- AND, there is also /entry/control/monitor_counts !!! Which one is the correct value? Which will NICE write
823        // TODO -- what do I really need to save?
824       
825        NVAR gDoMonitorNormalization = root:Packages:NIST:VSANS:Globals:gDoMonitorNormalization
826        if (gDoMonitorNormalization == 1)
827               
828                Variable monCount,savedMonCount
829                defmon=1e8                      //default monitor counts
830//              monCount = V_getControlMonitorCount(fname)                      // TODO -- this is read in since VCALC fakes this on output
831                monCount = V_getBeamMonNormData(fname)          // TODO -- I think this is the *real* one to read
832                savedMonCount   = monCount
833                scale = defMon/monCount         // scale factor to MULTIPLY data by to rescale to defmon
834
835                // write to newType=fname will put new values in the destination WORK folder
836                V_writeBeamMonNormSaved_count(fname,savedMonCount)                      // save the true count
837                V_writeBeamMonNormData(fname,defMon)            // mon ct is now 10^8
838                                       
839//                      // TODO
840//                      // the low efficiency monitor, expect to use this for white beam mode
841//                              -- need code switch here to determine which monitor to use.
842//
843//                      V_getBeamMonLowData(fname)
844//                      V_getBeamMonLowSaved_count(fname)       
845                       
846                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
847                        detStr = StringFromList(ii, ksDetectorListAll, ";")
848                        Wave w = V_getDetectorDataW(fname,detStr)
849                        Wave w_err = V_getDetectorDataErrW(fname,detStr)
850
851                        // do the calculation right here. It's a simple scaling and not worth sending to another function.     
852                        //scale the data and error to the default monitor counts
853               
854//
855                        w *= scale
856                        w_err *= scale          //assumes total monitor count is so large there is essentially no error
857
858                        // TODO
859                        // -- do I want to update and save the integrated detector count?
860                        // I can't trust the integrated count for the back detector (ever) - since the
861                        // data is shifted for registry and some (~ 10% of the detector) is lost
862                        // also, ML panel has the wrong value (Oct-nov 2019) and I don't know why. so sum the data
863                        //
864                        Variable integratedCount = sum(w)
865                        V_writeDet_IntegratedCount(fname,detStr,integratedCount)                //already the scaled value for counts
866//                      Variable integratedCount = V_getDet_IntegratedCount(fname,detStr)
867//                      V_writeDet_IntegratedCount(fname,detStr,integratedCount*scale)
868
869                endfor
870        else
871                Print "Monitor Normalization correction NOT DONE"
872        endif
873       
874       
875        // (not done) angle dependent efficiency correction
876        NVAR doEfficiency = root:Packages:NIST:VSANS:Globals:gDoDetectorEffCor
877
878//     
879        //reset the current displaytype to "newtype"
880        String/G root:Packages:NIST:VSANS:Globals:gCurDispType=newType
881       
882        //return to root folder (redundant)
883        SetDataFolder root:
884       
885        Return(0)
886End
887
888//
889//will "ADD" the current contents of the RAW folder to the newType work folder
890//
891// - used when adding multiple runs together
892//(the function V_Raw_to_work(type) makes a fresh workfile)
893//
894//the current display type is updated to newType (global)
895Function V_Add_raw_to_work(newType)
896        String newType
897       
898        String destPath=""
899        // if the desired workfile doesn't exist, let the user know, and just make a new one
900        if(WaveExists($("root:Packages:NIST:VSANS:"+newType + ":entry:instrument:detector_FL:data")) == 0)
901                Print "There is no old work file to add to - a new one will be created"
902                //call V_Raw_to_work(), then return from this function
903                V_Raw_to_Work(newType)
904                Return(0)               //does not generate an error - a single file was converted to work.newtype
905        Endif
906
907
908        // convert the RAW data to a WORK file.
909        // this will do all of the necessary corrections to the data
910        // put this in some separate work folder that can be cleaned out at the end (ADJ)
911        String tmpType="ADJ"
912       
913        //this step removes the read noise from the back so that neither added file will have this constant
914        V_Raw_to_Work(tmpType)         
915                       
916        //now make references to data in newType folder
917        DestPath="root:Packages:NIST:VSANS:"+newType   
918
919
920/////////////////
921//fields that need to be added together
922// entry block
923        // collection_time              V_getCollectionTime(fname)              V_putCollectionTime(fname,val)
924
925// instrument block
926        // beam_monitor_norm
927                // data (this will be 1e8)                              V_getBeamMonNormData(fname)             V_putBeamMonNormData(fname,val)
928                // saved_count (this is the original monitor count)  V_getBeamMonNormSaved_count(fname)         V_putBeamMonNormSaved_count(fname,val)
929
930        // for each detector
931        // data         V_getDetectorDataW(fname,detStr)
932        // integrated_count             V_getDet_IntegratedCount(fname,detStr)   V_putDet_IntegratedCount(fname,detStr,val)
933        // linear_data           V_getDetectorLinearDataW(fname,detStr)
934        // RECALCULATE (or add properly) linear_data_error              V_getDetectorDataErrW(fname,detStr)
935
936
937// control block (these may not actually be used?)
938        // count_time                           V_getCount_time(fname)                                          V_putCount_time(fname,val)
939        // detector_counts              V_getDetector_counts(fname)                             V_putDetector_counts(fname,val)
940        // monitor_counts               V_getControlMonitorCount(fname)         V_putControlMonitorCount(fname,val)
941
942// sample block - nothing
943// reduction block - nothing
944// user block - nothing
945
946// ?? need to add the file name to a list of what was actually added - so it will be saved with I(q)
947//
948////////////////////
949
950//      total_mon = realsread[1]        //saved monitor count
951//      uscale = total_mon/defmon               //unscaling factor
952//      total_det = uscale*realsread[2]         //unscaled detector count
953//      total_trn = uscale*realsread[39]        //unscaled trans det count
954//      total_numruns = integersread[3] //number of runs in workfile
955//      total_rtime = integersread[2]           //total counting time in workfile
956       
957
958        String detStr
959       
960        Variable saved_mon_dest,scale_dest,saved_mon_tmp,scale_tmp
961        Variable collection_time_dest,collection_time_tmp,count_time_dest,count_time_tmp
962        Variable detCount_dest,detCount_tmp,det_integrated_ct_dest,det_integrated_ct_tmp
963        Variable ii,new_scale,defMon
964       
965        defMon=1e8                      //default monitor counts
966
967        // find the scaling factors, one for each folder
968        saved_mon_dest = V_getBeamMonNormSaved_count(newType)
969        scale_dest = saved_mon_dest/defMon              //un-scaling factor
970       
971        saved_mon_tmp = V_getBeamMonNormSaved_count(tmpType)
972        scale_tmp = saved_mon_tmp/defMon                        //un-scaling factor
973
974        new_scale = defMon / (saved_mon_dest+saved_mon_tmp)
975       
976       
977        // get the count time for each (two locations)
978        collection_time_dest = V_getCollectionTime(newType)
979        collection_time_tmp = V_getCollectionTime(tmpType)
980       
981        count_time_dest = V_getCount_time(newType)
982        count_time_tmp = V_getCount_time(tmpType)
983       
984        detCount_dest = V_getDetector_counts(newType)
985        detCount_tmp = V_getDetector_counts(tmpType)
986
987// update the fields that are not in the detector blocks
988// in entry
989        V_putCollectionTime(newType,collection_time_dest+collection_time_tmp)
990
991// in control block
992        V_putCount_time(newType,count_time_dest+count_time_tmp)
993        V_putDetector_counts(newType,detCount_dest+detCount_tmp)
994        V_putControlMonitorCount(newType,saved_mon_dest+saved_mon_tmp)
995
996// TODO (DONE)
997// the new, unscaled monitor count was written to the control block, but it needs to be
998// written to the BeamMonNormSaved_count field instead, since this is where I read it from.
999// - so this worked in the past for adding two files, but fails on 3+
1000// x- write to the NormSaved_count field...
1001        V_writeBeamMonNormSaved_count(newType,saved_mon_dest+saved_mon_tmp)                     // save the true count
1002
1003
1004// now loop over all of the detector panels
1005        // data
1006        // data_err
1007        // integrated count
1008        // linear_data
1009        for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
1010                detStr = StringFromList(ii, ksDetectorListAll, ";")
1011               
1012                Wave data_dest = V_getDetectorDataW(newType,detStr)
1013                Wave data_err_dest = V_getDetectorDataErrW(newType,detStr)
1014                Wave linear_data_dest = V_getDetectorLinearDataW(newType,detStr)
1015                det_integrated_ct_dest = V_getDet_IntegratedCount(newType,detStr)
1016
1017                Wave data_tmp = V_getDetectorDataW(tmpType,detStr)
1018                Wave data_err_tmp = V_getDetectorDataErrW(tmpType,detStr)
1019                Wave linear_data_tmp = V_getDetectorLinearDataW(tmpType,detStr)
1020                det_integrated_ct_tmp = V_getDet_IntegratedCount(tmpType,detStr)
1021               
1022                // unscale the data arrays
1023                data_dest *= scale_dest
1024                data_err_dest *= scale_dest
1025                linear_data_dest *= scale_dest
1026               
1027                data_tmp *= scale_tmp
1028                data_err_tmp *= scale_tmp
1029                linear_data_tmp *= scale_tmp
1030
1031// TODO SRK-ERROR?
1032// the integrated count may not be correct (ML error Oct/Nov 2019)
1033// and is not correct for the back detector since some pixels were lost due to shifting the image registration
1034//                             
1035                // add them together, the dest is a wave so it is automatically changed in the "dest" folder
1036                V_putDet_IntegratedCount(tmpType,detStr,sum(data_dest)+sum(data_tmp))           // adds the unscaled data sums
1037//              V_putDet_IntegratedCount(tmpType,detStr,det_integrated_ct_dest+det_integrated_ct_tmp)           // wrong for "B", may be wrong for ML
1038                data_dest += data_tmp
1039                data_err_dest = sqrt(data_err_dest^2 + data_err_tmp^2)          // add in quadrature
1040                linear_data_dest += linear_data_tmp
1041               
1042                // now rescale the data_dest to the monitor counts
1043                data_dest *= new_scale
1044                data_err_dest *= new_scale
1045                linear_data_dest *= new_scale
1046               
1047        endfor
1048
1049       
1050        //Add the added raw filename to the list of files in the workfile
1051        SVAR fileList_dest = $("root:Packages:NIST:VSANS:"+newType+":gFileList")
1052        SVAR fileList_tmp = $("root:Packages:NIST:VSANS:"+tmpType+":gFileList")
1053       
1054        fileList_dest += ";" + fileList_tmp
1055       
1056        //reset the current display type to "newtype"
1057        SVAR gCurDispType = root:Packages:NIST:VSANS:Globals:gCurDispType
1058        gCurDispType = newType
1059       
1060        //return to root folder (redundant)
1061        SetDataFolder root:
1062       
1063        Return(0)
1064End
1065
1066
1067//used for adding DRK (beam shutter CLOSED) data to a workfile
1068//force the monitor count to 1, since it's irrelevant
1069// run data through normal "add" step, then unscale default monitor counts
1070//to get the data back on a simple time basis
1071//
1072Function V_Raw_to_Work_NoNorm(type)
1073        String type
1074       
1075        WAVE reals=$("root:Packages:NIST:RAW:realsread")
1076        reals[1]=1              //true monitor counts, still in raw
1077        V_Raw_to_work(type)
1078        //data is now in "type" folder
1079        WAVE data=$("root:Packages:NIST:"+type+":linear_data")
1080        WAVE data_copy=$("root:Packages:NIST:"+type+":data")
1081        WAVE data_err=$("root:Packages:NIST:"+type+":linear_data_error")
1082        WAVE new_reals=$("root:Packages:NIST:"+type+":realsread")
1083       
1084        Variable norm_mon,tot_mon,scale
1085       
1086        norm_mon = new_reals[0]         //should be 1e8
1087        tot_mon = new_reals[1]          //should be 1
1088        scale= norm_mon/tot_mon
1089       
1090        data /= scale           //unscale the data
1091        data_err /= scale
1092       
1093        // to keep "data" and linear_data in sync
1094        data_copy = data
1095       
1096        return(0)
1097End
1098
1099//used for adding DRK (beam shutter CLOSED) data to a workfile
1100//force the monitor count to 1, since it's irrelevant
1101// run data through normal "add" step, then unscale default monitor counts
1102//to get the data back on a simple time basis
1103//
1104Function V_Add_Raw_to_Work_NoNorm(type)
1105        String type
1106       
1107        WAVE reals=$("root:Packages:NIST:RAW:realsread")
1108        reals[1]=1              //true monitor counts, still in raw
1109        V_Add_Raw_to_work(type)
1110        //data is now in "type" folder
1111        WAVE data=$("root:Packages:NIST:"+type+":linear_data")
1112        WAVE data_copy=$("root:Packages:NIST:"+type+":data")
1113        WAVE data_err=$("root:Packages:NIST:"+type+":linear_data_error")
1114        WAVE new_reals=$("root:Packages:NIST:"+type+":realsread")
1115       
1116        Variable norm_mon,tot_mon,scale
1117       
1118        norm_mon = new_reals[0]         //should be 1e8
1119        tot_mon = new_reals[1]          //should be equal to the number of runs (1 count per run)
1120        scale= norm_mon/tot_mon
1121       
1122        data /= scale           //unscale the data
1123        data_err /= scale
1124       
1125        // to keep "data" and linear_data in sync
1126        data_copy = data
1127       
1128        return(0)
1129End
1130
Note: See TracBrowser for help on using the repository browser.