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

Last change on this file since 1117 was 1117, checked in by srkline, 4 years ago

extensive changes to accomodate 1x1 binning of the HighRes? detector. It is implemented as a global flag. Currently only 4x4 and 1x1 are allowed. 1x1 has never been tested in reality, only simulated data - so my assumed dimensions may not be correct. look for TODOHIGHRES in the file for places that may need to be updated for different file dimensions. Testing of the simulated data is proving to be excruciatingly slow, but passable for a test. Speed optimization will be needed if this is the final solution. Memory management will also be an issue since every "copy" of the highRes matrix is enormous. Carry as few of these around as possible in the future to keep the experiment size to something reasonable.

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