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

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

two significant changes:

1- in the averaging, before the data is finally written out, duplicate q-values (within 0.1%) are averaged so that duplicate q-values do not appear in the data file. This has a very bad effect on the calculation of the smearing matrix (dq=0).

2-for the data from the back (high res) detector, two steps have been added to the processing at the point where the data is converted to WORK. First, a constant read noise value of 200 cts/pixel is subtracted (found from average of multiple runs with beam off) Second, a 3x3 median filter is applied to the whole image to eliminate the stray bright pixels. Some are saturatd, some are simply outliers. Very effective.

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