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

Last change on this file since 1077 was 1077, checked in by srkline, 5 years ago

more attempts to correct the slit binning - to get the absolute scaling correct - still in progress.

File size: 34.9 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                // everything on the top level
111                V_DuplicateDataFolder($(fromDF+":entry"),fromStr,toStr,0,"",0)  //no recursion here
112                // control
113                V_DuplicateDataFolder($(fromDF+":entry:control"),fromStr,toStr,0,"",1)  //yes recursion here
114                // instrument
115                V_DuplicateDataFolder($(fromDF+":entry:instrument"),fromStr,toStr,0,"",1)       //yes recursion here
116                // reduction
117                V_DuplicateDataFolder($(fromDF+":entry:reduction"),fromStr,toStr,0,"",1)        //yes recursion here
118                // sample
119                V_DuplicateDataFolder($(fromDF+":entry:sample"),fromStr,toStr,0,"",1)   //yes recursion here
120                // user
121                V_DuplicateDataFolder($(fromDF+":entry:user"),fromStr,toStr,0,"",1)     //yes recursion here
122               
123        endif   
124       
125        return(0)
126end
127
128
129
130Function V_KillWavesInFolder(folderStr)
131        String folderStr
132       
133        if(DataFolderExists(folderStr) && strlen(folderStr) != 0)
134                SetDataFolder $folderStr
135                KillWaves/A/Z
136        endif
137       
138        SetDataFolder root:
139        return(0)
140end
141
142////////
143// see the help entry for IndexedDir for help on (possibly) how to do this faster
144// -- see  Function ScanDirectories(pathName, printDirNames)
145//
146
147
148// from IgorExchange On July 17th, 2011 jtigor
149// started from "Recursively List Data Folder Contents"
150// Posted July 15th, 2011 by hrodstein
151//
152//
153//
154Proc V_CopyWorkFolderProc(dataFolderStr, fromStr, toStr, level, sNBName, recurse)
155        String dataFolderStr="root:Packages:NIST:VSANS:RAW"
156        String fromStr = "RAW"
157        String toStr="SAM"
158        Variable level=0
159        String sNBName="DataFolderTree"
160        Variable recurse = 1
161       
162        V_DuplicateDataFolder($dataFolderStr, fromStr, toStr, level, sNBName, recurse)
163
164
165end
166
167// ListDataFolder(dfr, level)
168// Recursively lists objects in data folder.
169// Pass data folder path for dfr and 0 for level.
170// pass level == 0 for the first call
171//  sNBName = "" prints nothing. any name will generate a notebook
172//
173// recurse == 0 will do only the specified folder, anything else will recurse all levels
174// toStr is the string name of the top-level folder only, not the full path
175//
176//
177Function V_DuplicateDataFolder(dfr, fromStr, toStr, level, sNBName,recurse)
178        DFREF dfr
179        String fromStr
180        String toStr
181        Variable level                  // Pass 0 to start
182        String sNBName
183        Variable recurse
184 
185        String name
186        String dfName
187        String sString
188       
189        String toDF = ""
190 
191        if (level == 0)         // this is the data folder, generate if needed in the destination
192                name = GetDataFolder(1, dfr)
193//              sPrintf sString, "%s (data folder)\r", name
194                toDF = ReplaceString(fromStr,name,toStr,1)              // case-sensitive replace
195                sprintf sString, "NewDataFolder/O %s\r",toDF
196                NewDataFolder/O $(RemoveEnding(toDF,":"))                       // remove trailing semicolon if it's there
197               
198                V_WriteBrowserInfo_test(sString, 1, sNBName)
199        endif
200 
201        dfName = GetDataFolder(1, dfr)
202        toDF = ReplaceString(fromStr,dfName,toStr,1)            // case-sensitive replace
203        Variable i
204 
205        String indentStr = "\t"
206        for(i=0; i<level; i+=1)
207                indentStr += "\t"
208        endfor
209 
210        Variable numWaves = CountObjectsDFR(dfr, 1)
211        for(i=0; i<numWaves; i+=1)
212                name = GetIndexedObjNameDFR(dfr, 1, i)
213                //
214                // wave type does not matter now. Duplicate does not care
215                //
216                sPrintf sString, "Duplicate/O  %s,  %s\r",dfName+name,toDF+name
217                Duplicate/O $(dfName+name),$(toDF+name)
218               
219                V_WriteBrowserInfo_test(sString, 2, sNBName)
220        endfor 
221 
222        Variable numNumericVariables = CountObjectsDFR(dfr, 2) 
223        for(i=0; i<numNumericVariables; i+=1)
224                name = GetIndexedObjNameDFR(dfr, 2, i)
225                sPrintf sString, "%s%s (numeric variable)\r", indentStr, name
226                V_WriteBrowserInfo_test(sString, 3, sNBName)
227        endfor 
228 
229        Variable numStringVariables = CountObjectsDFR(dfr, 3)   
230        for(i=0; i<numStringVariables; i+=1)
231                name = GetIndexedObjNameDFR(dfr, 3, i)
232                sPrintf sString, "%s%s (string variable)\r", indentStr, name
233                V_WriteBrowserInfo_test(sString, 4, sNBName)
234        endfor 
235
236        if(recurse)
237                Variable numDataFolders = CountObjectsDFR(dfr, 4)       
238                for(i=0; i<numDataFolders; i+=1)
239                        name = GetIndexedObjNameDFR(dfr, 4, i)
240//                      sPrintf sString, "%s%s (data folder)\r", indentStr, name
241                         dfName = GetDataFolder(1, dfr)
242                         
243                        toDF = ReplaceString(fromStr,dfName,toStr,1)            // case-sensitive replace
244                        sprintf sString, "NewDataFolder/O %s\r",toDF+name
245                        NewDataFolder/O $(toDF+name)
246                       
247                       
248                        V_WriteBrowserInfo_test(sString, 1, sNBName)
249                        DFREF childDFR = dfr:$(name)
250                        V_DuplicateDataFolder(childDFR, fromStr, toStr, level+1, sNBName, recurse)
251                endfor 
252        endif
253         
254
255End
256
257
258// ListDataFolder(dfr, level)
259// Recursively lists objects in data folder.
260// Pass data folder path for dfr and 0 for level.
261// pass level == 0 for the first call
262//  sNBName = "" prints nothing. any name will generate a notebook
263//
264// recurse == 0 will do only the specified folder, anything else will recurse all levels
265// toStr is the string name of the top-level folder only, not the full path
266//
267//
268Function V_KillWavesFullTree(dfr, fromStr, level, sNBName,recurse)
269        DFREF dfr
270        String fromStr
271//      String toStr
272        Variable level                  // Pass 0 to start
273        String sNBName
274        Variable recurse
275 
276        String name
277        String dfName
278        String sString
279       
280        String toDF = ""
281 
282        if (level == 0)         // this is the data folder, generate if needed in the destination
283                name = GetDataFolder(1, dfr)
284                sPrintf sString, "%s (data folder)\r", name
285//              toDF = ReplaceString(fromStr,name,toStr,1)              // case-sensitive replace
286//              sprintf sString, "NewDataFolder/O %s\r",toDF
287//              NewDataFolder/O $(RemoveEnding(toDF,":"))                       // remove trailing semicolon if it's there
288               
289                V_WriteBrowserInfo_test(sString, 1, sNBName)
290        endif
291 
292        dfName = GetDataFolder(1, dfr)
293//      toDF = ReplaceString(fromStr,dfName,toStr,1)            // case-sensitive replace
294        Variable i
295 
296        String indentStr = "\t"
297        for(i=0; i<level; i+=1)
298                indentStr += "\t"
299        endfor
300 
301        Variable numWaves = CountObjectsDFR(dfr, 1)
302        for(i=0; i<numWaves; i+=1)
303                name = GetIndexedObjNameDFR(dfr, 1, i)
304                //
305                // wave type does not matter now. Kill does not care
306                //
307                sPrintf sString, "Killing  %s\r",dfName+name
308                KillWaves/Z $(dfName+name)
309               
310                V_WriteBrowserInfo_test(sString, 2, sNBName)
311        endfor 
312 
313 // now kill the data folder if possible
314        KillDataFolder/Z $dfName
315       
316       
317        Variable numNumericVariables = CountObjectsDFR(dfr, 2) 
318        for(i=0; i<numNumericVariables; i+=1)
319                name = GetIndexedObjNameDFR(dfr, 2, i)
320                sPrintf sString, "%s%s (numeric variable)\r", indentStr, name
321                V_WriteBrowserInfo_test(sString, 3, sNBName)
322        endfor 
323 
324        Variable numStringVariables = CountObjectsDFR(dfr, 3)   
325        for(i=0; i<numStringVariables; i+=1)
326                name = GetIndexedObjNameDFR(dfr, 3, i)
327                sPrintf sString, "%s%s (string variable)\r", indentStr, name
328                V_WriteBrowserInfo_test(sString, 4, sNBName)
329        endfor 
330
331        if(recurse)
332                Variable numDataFolders = CountObjectsDFR(dfr, 4)       
333                for(i=0; i<numDataFolders; i+=1)
334                        name = GetIndexedObjNameDFR(dfr, 4, i)
335                        sPrintf sString, "%s%s (data folder)\r", indentStr, name
336                         dfName = GetDataFolder(1, dfr)
337                         
338//                      toDF = ReplaceString(fromStr,dfName,toStr,1)            // case-sensitive replace
339//                      sprintf sString, "NewDataFolder/O %s\r",toDF+name
340//                      NewDataFolder/O $(toDF+name)
341                       
342                       
343                        V_WriteBrowserInfo_test(sString, 1, sNBName)
344                        DFREF childDFR = dfr:$(name)
345                        V_KillWavesFullTree(childDFR, fromStr, level+1, sNBName, recurse)
346                endfor 
347        endif
348         
349
350End
351 
352
353
354Function V_WriteBrowserInfo_test(sString, vType, sNBName)
355        String sString
356        Variable vType
357        String sNBName
358 
359        if(strlen(sNBName) == 0)
360//              print sString
361                return 0
362        endif
363        DoWindow $sNBName
364        if(V_flag != 1)
365                NewNoteBook/F=0 /N=$sNBName /V=1 as sNBName
366        else
367                DoWindow/F $sNBName
368        endif
369        Notebook $sNBName selection={endOfFile, endOfFile}
370        if(vType == 1)          // a data folder
371//              Notebook $sNBName fstyle=1
372                Notebook $sNBName text=sString
373//              Notebook $sNBName fstyle=-1
374        else
375                Notebook $sNBName text=sString 
376        endif
377 
378End
379
380///////////////////////////////
381
382
383//
384// given the folder, duplicate the data -> linear_data and generate the error
385//
386// x- do I want to use different names here? If it turns out that I don't need to drag a copy of
387//    the data around as "linear_data", then I can eliminate that, and rename the error wave
388// x- be sure the data is either properly written as 2D in the file, or converted to 2D before
389//    duplicating here
390// x- ? do I recast to DP here? No- V_MakeDataWaves_DP() is called directly prior to this call, so data
391//     coming in is already DP
392Function V_MakeDataError(folderStr)
393        String folderStr
394       
395        SetDataFolder $folderStr
396        Wave data=data
397        Duplicate/O data linear_data                    // at this point, the data is still the raw data, and is linear_data
398       
399        // proper error for counting statistics, good for low count values too
400        // rather than just sqrt(n)
401        // see N. Gehrels, Astrophys. J., 303 (1986) 336-346, equation (7)
402        // for S = 1 in eq (7), this corresponds to one sigma error bars
403        Duplicate/O linear_data linear_data_error
404        linear_data_error = 1 + sqrt(linear_data + 0.75)                               
405        //
406       
407        SetDataFolder root:
408        return(0)
409End
410
411
412/////////////////////
413
414
415
416//testing procedure
417// TODO -- can't duplicate this with another proceudre, but if I change the name of the variable
418//   "newType" to "type", then when Raw_to_work() gets to CopyHDFToWorkFolder(), the KillDataFolder/Z
419//   line fails (but reports no error), then DuplicateDataFolder fails, and reports an error. Trying
420//   to simplify this condition, I can't reproduce the bug for WM...
421Proc V_Convert_to_Workfile(newtype, doadd)
422        String newtype,doadd
423        Prompt newtype,"WORK data type",popup,"SAM;EMP;BGD;ADJ;"
424        Prompt doadd,"Add to current WORK contents?",popup,"No;Yes;"
425       
426        //macro will take whatever is in RAW folder and "ADD" it to the folder specified
427        //in the popup menu
428       
429        //"add" = yes/no, don't add to previous runs
430        //switch here - two separate functions to avoid (my) confusion
431        Variable err// = Raw_to_work(newtype)
432        if(cmpstr(doadd,"No")==0)
433                //don't add to prev work contents, copy RAW contents to work and convert
434                err = V_Raw_to_work(newtype)
435        else
436                //yes, add RAW to the current work folder contents
437                Abort "Adding RAW data files is currently unsupported"
438                err = V_Add_raw_to_work(newtype)
439        endif
440       
441        String newTitle = "WORK_"+newtype
442        DoWindow/F VSANS_Data
443        DoWindow/T VSANS_Data, newTitle
444        KillStrings/Z newTitle
445       
446        //need to update the display with "data" from the correct dataFolder
447        V_UpdateDisplayInformation(newtype)
448       
449End
450
451
452//
453// THIS IS THE MAJOR ROUTINE TO APPLY DATA CORRECTIONS
454//
455//will copy the current contents of the RAW folder to the newType work folder
456//and do the geometric corrections and normalization to monitor counts
457//(the function Add_Raw_to_work(type) adds multiple runs together - and is LOW priority)
458//
459//the current display type is updated to newType (global)
460//
461Function V_Raw_to_work(newType)
462        String newType
463       
464        Variable deadTime,defmon,total_mon,total_det,total_trn,total_numruns,total_rtime
465        Variable ii,jj,itim,cntrate,dscale,scale,uscale
466        String destPath
467       
468        String fname = newType
469        String detStr
470        Variable ctTime
471
472        //initialize values before normalization
473        total_mon=0
474        total_det=0
475        total_trn=0
476        total_numruns=0
477        total_rtime=0
478       
479        //Not adding multiple runs, so wipe out the old contents of the work folder and
480        // replace with the contents of raw
481
482        destPath = "root:Packages:NIST:VSANS:" + newType
483       
484        //copy from current dir (RAW) to work, defined by newType
485        V_CopyHDFToWorkFolder("RAW",newType)
486       
487        // now work with the waves from the destination folder.
488       
489        // apply corrections ---
490        // switches to control what is done, don't do the transmission correction for the BGD measurement
491        // start with the DIV correction, before conversion to mm
492        // then do all of the other corrections, order doesn't matter.
493        // rescaling to default monitor counts however, must be LAST.
494
495// each correction must loop over each detector. tedious.
496
497        // (0) Redimension the data waves in the destination folder
498        //     so that they are DP, not integer
499        // TODO
500        // -- currently only redimensioning the data and linear_data_error - What else???
501        // -- ?? some of this is done at load time for RAW data. shouldn't be an issue to re-do the redimension
502        for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
503                detStr = StringFromList(ii, ksDetectorListAll, ";")
504                Wave w = V_getDetectorDataW(fname,detStr)
505                Wave w_err = V_getDetectorDataErrW(fname,detStr)
506                Redimension/D w,w_err
507        endfor
508       
509       
510        // (1) DIV correction
511        // do this in terms of pixels.
512        // TODO : This must also exist at the time the first work folder is generated.
513        //   So it must be in the user folder at the start of the experiment, and defined.
514        NVAR gDoDIVCor = root:Packages:NIST:VSANS:Globals:gDoDIVCor
515        if (gDoDIVCor == 1)
516                // need extra check here for file existence
517                // if not in DIV folder, load.
518                // if unable to load, skip correction and report error (Alert?) (Ask to Load?)
519                Print "Doing DIV correction"// for "+ detStr
520                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
521                        detStr = StringFromList(ii, ksDetectorListAll, ";")
522                        Wave w = V_getDetectorDataW(fname,detStr)
523                        Wave w_err = V_getDetectorDataErrW(fname,detStr)
524                       
525                        V_DIVCorrection(w,w_err,detStr,newType)
526                endfor
527        else
528                Print "DIV correction NOT DONE"         // not an error since correction was unchecked
529        endif
530       
531        // (2) non-linear correction   
532        // TODO:
533        // x-  the "B" detector is calculated in its own routines
534        // -- document what is generated here:
535        //    **in each detector folder: data_realDistX and data_realDistY (2D waves of the mm? position of the pixel)
536        // x- these spatial calculations ARE DONE as the RAW data is loaded. It allows the RAW
537        //    data to be properly displayed, but without all of the (complete) set of detector corrections
538        // * the corrected distances are calculated into arrays, but nothing is done with them yet
539        // * there is enough information now to calculate the q-arrays, so it is done now
540        // - other corrections may modify the data, this calculation does NOT modify the data
541        NVAR gDoNonLinearCor = root:Packages:NIST:VSANS:Globals:gDoNonLinearCor
542        // generate a distance matrix for each of the detectors
543        if (gDoNonLinearCor == 1)
544                Print "Doing Non-linear correction"// for "+ detStr
545                for(ii=0;ii<ItemsInList(ksDetectorListNoB);ii+=1)
546                        detStr = StringFromList(ii, ksDetectorListNoB, ";")
547                        Wave w = V_getDetectorDataW(fname,detStr)
548//                      Wave w_err = V_getDetectorDataErrW(fname,detStr)
549                        Wave w_calib = V_getDetTube_spatialCalib(fname,detStr)
550                        Variable tube_width = V_getDet_tubeWidth(fname,detStr)
551                        V_NonLinearCorrection(fname,w,w_calib,tube_width,detStr,destPath)
552
553                        //(2.4) Convert the beam center values from pixels to mm
554                        // TODO -- there needs to be a permanent location for these values??
555                        //
556                                // TODO
557                                // -- the beam center value in mm needs to be present - it is used in calculation of Qvalues
558                                // -- but having both the same is wrong...
559                                // -- the pixel value is needed for display of the panels
560                                if(kBCTR_CM)
561                                        //V_ConvertBeamCtr_to_mm(folder,detStr,destPath)
562                                        //
563       
564                                        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
565                                        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
566                                        WAVE x_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
567                                        WAVE y_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
568                                        x_mm[0] = V_getDet_beam_center_x(fname,detStr) * 10             // convert cm to mm
569                                        y_mm[0] = V_getDet_beam_center_y(fname,detStr) * 10             // convert cm to mm
570                                       
571                                        // TODO:::
572                                // now I need to convert the beam center in mm to pixels
573                                // and have some rational place to look for it...
574                                        V_ConvertBeamCtr_to_pix(fname,detStr,destPath)
575                                else
576                                        // beam center is in pixels, so use the old routine
577                                        V_ConvertBeamCtr_to_mm(fname,detStr,destPath)
578                                endif           
579                                                       
580                        // (2.5) Calculate the q-values
581                        // calculating q-values can't be done unless the non-linear corrections are calculated
582                        // so go ahead and put it in this loop.
583                        // TODO :
584                        // -- make sure that everything is present before the calculation
585                        // -- beam center must be properly defined in terms of real distance
586                        // -- distances/zero location/ etc. must be clearly documented for each detector
587                        //      ** this assumes that NonLinearCorrection() has been run to generate data_RealDistX and Y
588                        // ** this routine Makes the waves QTot, qx, qy, qz in each detector folder.
589                        //
590                        V_Detector_CalcQVals(fname,detStr,destPath)
591                       
592                endfor
593               
594                //"B" is separate
595                Wave w = V_getDetectorDataW(fname,"B")
596                Wave cal_x = V_getDet_cal_x(fname,"B")
597                Wave cal_y = V_getDet_cal_y(fname,"B")
598                V_NonLinearCorrection_B(fname,w,cal_x,cal_y,"B",destPath)
599                V_ConvertBeamCtr_to_mmB(fname,"B",destPath)
600                V_Detector_CalcQVals(fname,"B",destPath)
601               
602        else
603                Print "Non-linear correction NOT DONE"
604        endif
605
606
607        // (3) dead time correction
608        // TODO:
609        // -- test for correct operation
610        // x- loop over all of the detectors
611        // x- B detector is a special case (do separately, then loop over NoB)
612        // -- this DOES alter the data
613        // -- verify the error propagation (not done yet)
614        //
615        Variable countRate
616        NVAR gDoDeadTimeCor = root:Packages:NIST:VSANS:Globals:gDoDeadTimeCor
617        if (gDoDeadTimeCor == 1)
618                Print "Doing DeadTime correction"// for "+ detStr
619                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
620                        detStr = StringFromList(ii, ksDetectorListAll, ";")
621                        Wave w = V_getDetectorDataW(fname,detStr)
622                        Wave w_err = V_getDetectorDataErrW(fname,detStr)
623                        ctTime = V_getCount_time(fname)
624
625                        if(cmpstr(detStr,"B") == 0)
626                                Variable b_dt = V_getDetector_deadtime_B(fname,detStr)
627                                // do the correction for the back panel
628                                countRate = sum(w,-inf,inf)/ctTime              //use sum of detector counts
629
630                                w = w/(1-countRate*b_dt)
631                                w_err = w_err/(1-countRate*b_dt)
632                                                               
633                        else
634                                // do the corrections for 8 tube panels
635                                Wave w_dt = V_getDetector_deadtime(fname,detStr)
636                                V_DeadTimeCorrectionTubes(w,w_err,w_dt,ctTime)
637                               
638                        endif
639                endfor
640               
641        else
642                Print "Dead Time correction NOT DONE"
643        endif   
644       
645
646        // (4) solid angle correction
647        //  -- this currently calculates the correction factor AND applys it to the data
648        //  -- as a result, the data values are very large since they are divided by a very small
649        //     solid angle per pixel. But all of the count values are now on the basis of
650        //    counts/(solid angle) --- meaning that they can all be binned together for I(q)
651        //    -and- - this is taken into account for absolute scaling (this part is already done)
652        //
653        // the solid angle correction is calculated for ALL detector panels.
654        NVAR gDoSolidAngleCor = root:Packages:NIST:VSANS:Globals:gDoSolidAngleCor
655        if (gDoSolidAngleCor == 1)
656                Print "Doing Solid Angle correction"// for "+ detStr
657                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
658                        detStr = StringFromList(ii, ksDetectorListAll, ";")
659                        Wave w = V_getDetectorDataW(fname,detStr)
660                        Wave w_err = V_getDetectorDataErrW(fname,detStr)
661                        // any other dimensions to pass in?
662                        V_SolidAngleCorrection(w,w_err,fname,detStr,destPath)
663                       
664                endfor
665        else
666                Print "Solid Angle correction NOT DONE"
667        endif   
668       
669               
670        // (5) angle-dependent tube shadowing
671        // TODO:
672        // -- not sure about this correction yet...
673        //
674        NVAR gDoTubeShadowCor = root:Packages:NIST:VSANS:Globals:gDoTubeShadowCor
675        if (gDoTubeShadowCor == 1)
676                Print "(stub)Tube Shadow correction"
677        else
678                Print "Tube shadowing correction NOT DONE"
679        endif   
680               
681        // (6) angle dependent transmission correction
682        // TODO:
683        // x- this still needs to be filled in
684        // -- still some debate of when/where in the corrections that this is best applied
685        //    - do it here, and it's done whether the output is 1D or 2D
686        //    - do it later (where SAMPLE information is used) since this section is ONLY instrument-specific
687        // -- verify that the calculation is correct
688        // -- verify that the error propagation (in 2D) is correct
689        //
690        NVAR gDoTrans = root:Packages:NIST:VSANS:Globals:gDoTransmissionCor
691        if (gDoTrans == 1)
692                Print "Doing Large-angle transmission correction"// for "+ detStr
693                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
694                        detStr = StringFromList(ii, ksDetectorListAll, ";")
695                        Wave w = V_getDetectorDataW(fname,detStr)
696                        Wave w_err = V_getDetectorDataErrW(fname,detStr)
697                       
698                        V_LargeAngleTransmissionCorr(w,w_err,fname,detStr,destPath)
699                       
700                endfor
701        else
702                Print "Sample Transmission correction NOT DONE"
703        endif   
704       
705        // (7) normalize to default monitor counts
706        // TODO(DONE) x- each detector is rescaled separately, but the rescaling factor is global (only one monitor!)
707        // TODO -- but there are TWO monitors - so how to switch?
708        //  --- AND, there is also /entry/control/monitor_counts !!! Which one is the correct value? Which will NICE write
709        // TODO -- what do I really need to save?
710       
711        NVAR gDoMonitorNormalization = root:Packages:NIST:VSANS:Globals:gDoMonitorNormalization
712        if (gDoMonitorNormalization == 1)
713               
714                Variable monCount,savedMonCount
715                defmon=1e8                      //default monitor counts
716                monCount = V_getMonitorCount(fname)                     // TODO -- this is read in since VCALC fakes this on output
717//              monCount = V_getBeamMonNormData(fname)          // TODO -- I think this is the *real* one to read
718                savedMonCount   = monCount
719                scale = defMon/monCount         // scale factor to MULTIPLY data by to rescale to defmon
720
721                // write to newType=fname will put new values in the destination WORK folder
722                V_writeBeamMonNormSaved_count(fname,savedMonCount)                      // save the true count
723                V_writeBeamMonNormData(fname,defMon)            // mon ct is now 10^8
724                                       
725//                      // TODO
726//                      // the low efficiency monitor, expect to use this for white beam mode
727//                              -- need code switch here to determine which monitor to use.
728//
729//                      V_getBeamMonLowData(fname)
730//                      V_getBeamMonLowSaved_count(fname)       
731                       
732                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
733                        detStr = StringFromList(ii, ksDetectorListAll, ";")
734                        Wave w = V_getDetectorDataW(fname,detStr)
735                        Wave w_err = V_getDetectorDataErrW(fname,detStr)
736
737                        // do the calculation right here. It's a simple scaling and not worth sending to another function.     
738                        //scale the data and error to the default monitor counts
739               
740//
741                        w *= scale
742                        w_err *= scale          //assumes total monitor count is so large there is essentially no error
743
744                        // TODO
745                        // -- do I want to update and save the integrated detector count?
746                        Variable integratedCount = V_getDet_IntegratedCount(fname,detStr)
747                        V_writeDet_IntegratedCount(fname,detStr,integratedCount*scale)
748
749                endfor
750        else
751                Print "Monitor Normalization correction NOT DONE"
752        endif
753       
754       
755        // (not done) angle dependent efficiency correction
756        NVAR doEfficiency = root:Packages:NIST:VSANS:Globals:gDoDetectorEffCor
757
758//
759///// these lines are if files are added together, not done (yet) for VSANS
760//
761//update totals to put in the work header (at the end of the function)
762//      total_mon += realsread[0]
763//
764//      total_det += dscale*realsread[2]
765//
766//      total_trn += realsread[39]
767//      total_rtime += integersread[2]
768//      total_numruns +=1
769//     
770
771//     
772        //reset the current displaytype to "newtype"
773        String/G root:Packages:NIST:VSANS:Globals:gCurDispType=newType
774       
775        //return to root folder (redundant)
776        SetDataFolder root:
777       
778        Return(0)
779End
780
781
782//will "ADD" the current contents of the RAW folder to the newType work folder
783//and will ADD the RAW contents to the existing content of the newType folder
784// - used when adding multiple runs together
785//(the function Raw_to_work(type) makes a fresh workfile)
786//
787//the current display type is updated to newType (global)
788Function V_Add_raw_to_work(newType)
789        String newType
790       
791        // NEW OCT 2014
792        // this corrects for adding raw data files with different attenuation   
793        // does nothing if the attenuation of RAW and destination are the same
794        NVAR doAdjustRAW_Atten = root:Packages:NIST:gDoAdjustRAW_Atten
795        if(doAdjustRAW_Atten)
796                V_Adjust_RAW_Attenuation(newType)
797        endif
798       
799        String destPath=""
800       
801        // if the desired workfile doesn't exist, let the user know, and just make a new one
802        if(WaveExists($("root:Packages:NIST:"+newType + ":data")) == 0)
803                Print "There is no old work file to add to - a new one will be created"
804                //call Raw_to_work(), then return from this function
805                V_Raw_to_Work(newType)
806                Return(0)               //does not generate an error - a single file was converted to work.newtype
807        Endif
808       
809        NVAR pixelsX = root:myGlobals:gNPixelsX
810        NVAR pixelsY = root:myGlobals:gNPixelsY
811       
812        //now make references to data in newType folder
813        DestPath="root:Packages:NIST:"+newType 
814        WAVE data=$(destPath +":linear_data")                   // these wave references point to the EXISTING work data
815        WAVE data_copy=$(destPath +":data")                     // these wave references point to the EXISTING work data
816        WAVE dest_data_err=$(destPath +":linear_data_error")                    // these wave references point to the EXISTING work data
817        WAVE/T textread=$(destPath + ":textread")
818        WAVE integersread=$(destPath + ":integersread")
819        WAVE realsread=$(destPath + ":realsread")
820       
821        Variable deadTime,defmon,total_mon,total_det,total_trn,total_numruns,total_rtime
822        Variable ii,jj,itim,cntrate,dscale,scale,uscale,wrk_beamx,wrk_beamy,xshift,yshift
823
824
825        defmon=1e8                      //default monitor counts
826       
827        //Yes, add to previous run(s) in work, that does exist
828        //use the actual monitor count run.savmon rather than the normalized monitor count
829        //in run.moncnt and unscale the work data
830       
831        total_mon = realsread[1]        //saved monitor count
832        uscale = total_mon/defmon               //unscaling factor
833        total_det = uscale*realsread[2]         //unscaled detector count
834        total_trn = uscale*realsread[39]        //unscaled trans det count
835        total_numruns = integersread[3] //number of runs in workfile
836        total_rtime = integersread[2]           //total counting time in workfile
837        //retrieve workfile beamcenter
838        wrk_beamx = realsread[16]
839        wrk_beamy = realsread[17]
840        //unscale the workfile data in "newType"
841        //
842        //check for log-scaling and adjust if necessary
843        // should not be needed now - using display flag instead
844//      ConvertFolderToLinearScale(newType)
845        //
846        //then unscale the data array
847        data *= uscale
848        dest_data_err *= uscale
849       
850        //DetCorr() has not been applied to the data in RAW , do it now in a local reference to the raw data
851        WAVE raw_data = $"root:Packages:NIST:RAW:linear_data"
852        WAVE raw_data_err = $"root:Packages:NIST:RAW:linear_data_error"
853        WAVE raw_reals =  $"root:Packages:NIST:RAW:realsread"
854        WAVE/T raw_text = $"root:Packages:NIST:RAW:textread"
855        WAVE raw_ints = $"root:Packages:NIST:RAW:integersread"
856       
857        //check for log-scaling of the raw data - make sure it's linear
858        // should not be needed now - using display flag instead
859//      ConvertFolderToLinearScale("RAW")
860       
861        // switches to control what is done, don't do the transmission correction for the BGD measurement
862        NVAR doEfficiency = root:Packages:NIST:gDoDetectorEffCor
863        NVAR gDoTrans = root:Packages:NIST:gDoTransmissionCor
864        Variable doTrans = gDoTrans
865        if(cmpstr("BGD",newtype) == 0)
866                doTrans = 0             //skip the trans correction for the BGD file but don't change the value of the global
867        endif   
868       
869        V_DetCorr(raw_data,raw_data_err,raw_reals,doEfficiency,doTrans) //applies correction to raw_data, and overwrites it
870       
871        //deadtime corrections to raw data
872        // TODO - do the tube correction for dead time now
873        deadTime = 1//DetectorDeadtime(raw_text[3],raw_text[9],dateAndTimeStr=raw_text[1],dtime=raw_reals[48])          //pick the correct detector deadtime, switch on date too
874        itim = raw_ints[2]
875        cntrate = sum(raw_data,-inf,inf)/itim           //080802 use data sum, rather than scaler value
876        dscale = 1/(1-deadTime*cntrate)
877
878#if (exists("ILL_D22")==6)
879        Variable tubeSum
880        // for D22 detector might need to use cntrate/128 as it is the tube response
881        for(ii=0;ii<pixelsX;ii+=1)
882                //sum the counts in each tube
883                tubeSum = 0
884                for(jj=0;jj<pixelsY;jj+=1)
885                        tubeSum += data[jj][ii]
886                endfor
887                // countrate in tube ii
888                cntrate = tubeSum/itim
889                // deadtime scaling in tube ii
890                dscale = 1/(1-deadTime*cntrate)
891                // multiply data[ii][] by the dead time
892                raw_data[][ii] *= dscale
893                raw_data_err[][ii] *= dscale
894        endfor
895#else
896        // dead time correction on all other RAW data, including NCNR
897        raw_data *= dscale
898        raw_data_err *= dscale
899#endif
900
901        //update totals by adding RAW values to the local ones (write to work header at end of function)
902        total_mon += raw_reals[0]
903
904        total_det += dscale*raw_reals[2]
905
906        total_trn += raw_reals[39]
907        total_rtime += raw_ints[2]
908        total_numruns +=1
909       
910        //do the beamcenter shifting if there is a mismatch
911        //and then add the two data sets together, changing "data" since it is the workfile data
912        xshift = raw_reals[16] - wrk_beamx
913        yshift = raw_reals[17] - wrk_beamy
914       
915        If((xshift != 0) || (yshift != 0))
916                DoAlert 1,"Do you want to ignore the beam center mismatch?"
917                if(V_flag==1)
918                        xshift=0
919                        yshift=0
920                endif
921        endif
922       
923        If((xshift == 0) && (yshift == 0))              //no shift, just add them
924                data += raw_data                //deadtime correction has already been done to the raw data
925                dest_data_err = sqrt(dest_data_err^2 + raw_data_err^2)                  // error of the sum
926        Endif
927       
928        //scale the data to the default montor counts
929        scale = defmon/total_mon
930        data *= scale
931        dest_data_err *= scale
932       
933        // keep "data" and linear_data in sync in the destination folder
934        data_copy = data
935       
936        //all is done, except for the bookkeeping of updating the header info in the work folder
937        textread[1] = date() + " " + time()             //date + time stamp
938        integersread[3] = total_numruns                                         //numruns = more than one
939        realsread[1] = total_mon                        //save the true monitor count
940        realsread[0] = defmon                                   //monitor ct = defmon
941        integersread[2] = total_rtime                   // total counting time
942        realsread[2] = scale*total_det                  //scaled detector counts
943        realsread[39] = scale*total_trn                 //scaled transmission counts
944       
945        //Add the added raw filename to the list of files in the workfile
946        String newfile = ";" + raw_text[0]
947        SVAR oldList = $(destPath + ":fileList")
948        String/G $(destPath + ":fileList") = oldList + newfile
949       
950        //reset the current display type to "newtype"
951        SVAR gCurDispType = root:Packages:NIST:VSANS:Globals:gCurDispType
952        gCurDispType = newType
953       
954        //return to root folder (redundant)
955        SetDataFolder root:
956       
957        Return(0)
958End
959
960
961//used for adding DRK (beam shutter CLOSED) data to a workfile
962//force the monitor count to 1, since it's irrelevant
963// run data through normal "add" step, then unscale default monitor counts
964//to get the data back on a simple time basis
965//
966Function V_Raw_to_Work_NoNorm(type)
967        String type
968       
969        WAVE reals=$("root:Packages:NIST:RAW:realsread")
970        reals[1]=1              //true monitor counts, still in raw
971        V_Raw_to_work(type)
972        //data is now in "type" folder
973        WAVE data=$("root:Packages:NIST:"+type+":linear_data")
974        WAVE data_copy=$("root:Packages:NIST:"+type+":data")
975        WAVE data_err=$("root:Packages:NIST:"+type+":linear_data_error")
976        WAVE new_reals=$("root:Packages:NIST:"+type+":realsread")
977       
978        Variable norm_mon,tot_mon,scale
979       
980        norm_mon = new_reals[0]         //should be 1e8
981        tot_mon = new_reals[1]          //should be 1
982        scale= norm_mon/tot_mon
983       
984        data /= scale           //unscale the data
985        data_err /= scale
986       
987        // to keep "data" and linear_data in sync
988        data_copy = data
989       
990        return(0)
991End
992
993//used for adding DRK (beam shutter CLOSED) data to a workfile
994//force the monitor count to 1, since it's irrelevant
995// run data through normal "add" step, then unscale default monitor counts
996//to get the data back on a simple time basis
997//
998Function V_Add_Raw_to_Work_NoNorm(type)
999        String type
1000       
1001        WAVE reals=$("root:Packages:NIST:RAW:realsread")
1002        reals[1]=1              //true monitor counts, still in raw
1003        V_Add_Raw_to_work(type)
1004        //data is now in "type" folder
1005        WAVE data=$("root:Packages:NIST:"+type+":linear_data")
1006        WAVE data_copy=$("root:Packages:NIST:"+type+":data")
1007        WAVE data_err=$("root:Packages:NIST:"+type+":linear_data_error")
1008        WAVE new_reals=$("root:Packages:NIST:"+type+":realsread")
1009       
1010        Variable norm_mon,tot_mon,scale
1011       
1012        norm_mon = new_reals[0]         //should be 1e8
1013        tot_mon = new_reals[1]          //should be equal to the number of runs (1 count per run)
1014        scale= norm_mon/tot_mon
1015       
1016        data /= scale           //unscale the data
1017        data_err /= scale
1018       
1019        // to keep "data" and linear_data in sync
1020        data_copy = data
1021       
1022        return(0)
1023End
1024
Note: See TracBrowser for help on using the repository browser.