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

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

Added functions to check that files in a protocol, including the sample data are all from the same configuration. In SANS, the only flag was a beam center mismatch. In VSANS, a more complete check of the configuration in necessary.

Changed all instances of calls to the read and normalize the monitor count to use the "norm"al monitor, not the value listed in the Control block.

File size: 35.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 - anything else???
501        // -- ?? some of this is done at load time for RAW data. Not a problem 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 each 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                detStr = "B"
596                Wave w = V_getDetectorDataW(fname,detStr)
597                Wave cal_x = V_getDet_cal_x(fname,detStr)
598                Wave cal_y = V_getDet_cal_y(fname,detStr)
599               
600                V_NonLinearCorrection_B(fname,w,cal_x,cal_y,detStr,destPath)
601
602// "B" is always naturally defined in terms of a pixel center. This can be converted to mm,
603// but the experiment will measure pixel x,y - just like ordela detectors.
604               
605//              if(kBCTR_CM)
606//
607//                      Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
608//                      Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
609//                      WAVE x_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
610//                      WAVE y_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
611//                      x_mm[0] = V_getDet_beam_center_x(fname,detStr) * 10             // convert cm to mm
612//                      y_mm[0] = V_getDet_beam_center_y(fname,detStr) * 10             // convert cm to mm
613//                     
614//                      // now I need to convert the beam center in mm to pixels
615//                      // and have some rational place to look for it...
616//                      V_ConvertBeamCtr_to_pixB(fname,detStr,destPath)
617//              else
618                        // beam center is in pixels, so use the old routine
619                        V_ConvertBeamCtr_to_mmB(fname,detStr,destPath)
620
621//              endif
622                V_Detector_CalcQVals(fname,detStr,destPath)
623               
624        else
625                Print "Non-linear correction NOT DONE"
626        endif
627
628
629        // (3) dead time correction
630        // DONE:
631        // x- test for correct operation
632        // x- loop over all of the detectors
633        // x- B detector is a special case (do separately, then loop over NoB)
634        // x- this DOES alter the data
635        // x- verify the error propagation (not done yet)
636        //
637        Variable countRate
638        NVAR gDoDeadTimeCor = root:Packages:NIST:VSANS:Globals:gDoDeadTimeCor
639        if (gDoDeadTimeCor == 1)
640                Print "Doing DeadTime correction"// for "+ detStr
641                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
642                        detStr = StringFromList(ii, ksDetectorListAll, ";")
643                        Wave w = V_getDetectorDataW(fname,detStr)
644                        Wave w_err = V_getDetectorDataErrW(fname,detStr)
645                        ctTime = V_getCount_time(fname)
646
647                        if(cmpstr(detStr,"B") == 0)
648                                Variable b_dt = V_getDetector_deadtime_B(fname,detStr)
649                                // do the correction for the back panel
650                                countRate = sum(w,-inf,inf)/ctTime              //use sum of detector counts
651
652                                w = w/(1-countRate*b_dt)
653                                w_err = w_err/(1-countRate*b_dt)
654                                                               
655                        else
656                                // do the corrections for 8 tube panels
657                                Wave w_dt = V_getDetector_deadtime(fname,detStr)
658                                V_DeadTimeCorrectionTubes(w,w_err,w_dt,ctTime)
659                               
660                        endif
661                endfor
662               
663        else
664                Print "Dead Time correction NOT DONE"
665        endif   
666       
667
668        // (4) solid angle correction
669        //  -- this currently calculates the correction factor AND applys it to the data
670        //  -- as a result, the data values are very large since they are divided by a very small
671        //     solid angle per pixel. But all of the count values are now on the basis of
672        //    counts/(solid angle) --- meaning that they can all be binned together for I(q)
673        //    -and- - this is taken into account for absolute scaling (this part is already done)
674        //
675        // the solid angle correction is calculated for ALL detector panels.
676        NVAR gDoSolidAngleCor = root:Packages:NIST:VSANS:Globals:gDoSolidAngleCor
677        if (gDoSolidAngleCor == 1)
678                Print "Doing Solid Angle correction"// for "+ detStr
679                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
680                        detStr = StringFromList(ii, ksDetectorListAll, ";")
681                        Wave w = V_getDetectorDataW(fname,detStr)
682                        Wave w_err = V_getDetectorDataErrW(fname,detStr)
683                        // any other dimensions to pass in?
684                        V_SolidAngleCorrection(w,w_err,fname,detStr,destPath)
685                       
686                endfor
687        else
688                Print "Solid Angle correction NOT DONE"
689        endif   
690       
691               
692        // (5) angle-dependent tube shadowing
693        // TODO:
694        // -- not sure about this correction yet...
695        //
696        NVAR gDoTubeShadowCor = root:Packages:NIST:VSANS:Globals:gDoTubeShadowCor
697        if (gDoTubeShadowCor == 1)
698                Print "(stub)Tube Shadow correction"
699        else
700                Print "Tube shadowing correction NOT DONE"
701        endif   
702               
703        // (6) angle dependent transmission correction
704        // TODO:
705        // x- this still needs to be filled in
706        // -- still some debate of when/where in the corrections that this is best applied
707        //    - do it here, and it's done whether the output is 1D or 2D
708        //    - do it later (where SAMPLE information is used) since this section is ONLY instrument-specific
709        // -- verify that the calculation is correct
710        // -- verify that the error propagation (in 2D) is correct
711        //
712        NVAR gDoTrans = root:Packages:NIST:VSANS:Globals:gDoTransmissionCor
713        if (gDoTrans == 1)
714                Print "Doing Large-angle transmission correction"// for "+ detStr
715                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
716                        detStr = StringFromList(ii, ksDetectorListAll, ";")
717                        Wave w = V_getDetectorDataW(fname,detStr)
718                        Wave w_err = V_getDetectorDataErrW(fname,detStr)
719                       
720                        V_LargeAngleTransmissionCorr(w,w_err,fname,detStr,destPath)
721                       
722                endfor
723        else
724                Print "Sample Transmission correction NOT DONE"
725        endif   
726       
727        // (7) normalize to default monitor counts
728        // TODO(DONE) x- each detector is rescaled separately, but the rescaling factor is global (only one monitor!)
729        // TODO -- but there are TWO monitors - so how to switch?
730        //  --- AND, there is also /entry/control/monitor_counts !!! Which one is the correct value? Which will NICE write
731        // TODO -- what do I really need to save?
732       
733        NVAR gDoMonitorNormalization = root:Packages:NIST:VSANS:Globals:gDoMonitorNormalization
734        if (gDoMonitorNormalization == 1)
735               
736                Variable monCount,savedMonCount
737                defmon=1e8                      //default monitor counts
738//              monCount = V_getControlMonitorCount(fname)                      // TODO -- this is read in since VCALC fakes this on output
739                monCount = V_getBeamMonNormData(fname)          // TODO -- I think this is the *real* one to read
740                savedMonCount   = monCount
741                scale = defMon/monCount         // scale factor to MULTIPLY data by to rescale to defmon
742
743                // write to newType=fname will put new values in the destination WORK folder
744                V_writeBeamMonNormSaved_count(fname,savedMonCount)                      // save the true count
745                V_writeBeamMonNormData(fname,defMon)            // mon ct is now 10^8
746                                       
747//                      // TODO
748//                      // the low efficiency monitor, expect to use this for white beam mode
749//                              -- need code switch here to determine which monitor to use.
750//
751//                      V_getBeamMonLowData(fname)
752//                      V_getBeamMonLowSaved_count(fname)       
753                       
754                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
755                        detStr = StringFromList(ii, ksDetectorListAll, ";")
756                        Wave w = V_getDetectorDataW(fname,detStr)
757                        Wave w_err = V_getDetectorDataErrW(fname,detStr)
758
759                        // do the calculation right here. It's a simple scaling and not worth sending to another function.     
760                        //scale the data and error to the default monitor counts
761               
762//
763                        w *= scale
764                        w_err *= scale          //assumes total monitor count is so large there is essentially no error
765
766                        // TODO
767                        // -- do I want to update and save the integrated detector count?
768                        Variable integratedCount = V_getDet_IntegratedCount(fname,detStr)
769                        V_writeDet_IntegratedCount(fname,detStr,integratedCount*scale)
770
771                endfor
772        else
773                Print "Monitor Normalization correction NOT DONE"
774        endif
775       
776       
777        // (not done) angle dependent efficiency correction
778        NVAR doEfficiency = root:Packages:NIST:VSANS:Globals:gDoDetectorEffCor
779
780//
781///// these lines are if files are added together, not done (yet) for VSANS
782//
783//update totals to put in the work header (at the end of the function)
784//      total_mon += realsread[0]
785//
786//      total_det += dscale*realsread[2]
787//
788//      total_trn += realsread[39]
789//      total_rtime += integersread[2]
790//      total_numruns +=1
791//     
792
793//     
794        //reset the current displaytype to "newtype"
795        String/G root:Packages:NIST:VSANS:Globals:gCurDispType=newType
796       
797        //return to root folder (redundant)
798        SetDataFolder root:
799       
800        Return(0)
801End
802
803
804//will "ADD" the current contents of the RAW folder to the newType work folder
805//and will ADD the RAW contents to the existing content of the newType folder
806// - used when adding multiple runs together
807//(the function Raw_to_work(type) makes a fresh workfile)
808//
809//the current display type is updated to newType (global)
810Function V_Add_raw_to_work(newType)
811        String newType
812       
813        // NEW OCT 2014
814        // this corrects for adding raw data files with different attenuation   
815        // does nothing if the attenuation of RAW and destination are the same
816        NVAR doAdjustRAW_Atten = root:Packages:NIST:gDoAdjustRAW_Atten
817        if(doAdjustRAW_Atten)
818                V_Adjust_RAW_Attenuation(newType)
819        endif
820       
821        String destPath=""
822       
823        // if the desired workfile doesn't exist, let the user know, and just make a new one
824        if(WaveExists($("root:Packages:NIST:"+newType + ":data")) == 0)
825                Print "There is no old work file to add to - a new one will be created"
826                //call Raw_to_work(), then return from this function
827                V_Raw_to_Work(newType)
828                Return(0)               //does not generate an error - a single file was converted to work.newtype
829        Endif
830       
831        NVAR pixelsX = root:myGlobals:gNPixelsX
832        NVAR pixelsY = root:myGlobals:gNPixelsY
833       
834        //now make references to data in newType folder
835        DestPath="root:Packages:NIST:"+newType 
836        WAVE data=$(destPath +":linear_data")                   // these wave references point to the EXISTING work data
837        WAVE data_copy=$(destPath +":data")                     // these wave references point to the EXISTING work data
838        WAVE dest_data_err=$(destPath +":linear_data_error")                    // these wave references point to the EXISTING work data
839        WAVE/T textread=$(destPath + ":textread")
840        WAVE integersread=$(destPath + ":integersread")
841        WAVE realsread=$(destPath + ":realsread")
842       
843        Variable deadTime,defmon,total_mon,total_det,total_trn,total_numruns,total_rtime
844        Variable ii,jj,itim,cntrate,dscale,scale,uscale,wrk_beamx,wrk_beamy,xshift,yshift
845
846
847        defmon=1e8                      //default monitor counts
848       
849        //Yes, add to previous run(s) in work, that does exist
850        //use the actual monitor count run.savmon rather than the normalized monitor count
851        //in run.moncnt and unscale the work data
852       
853        total_mon = realsread[1]        //saved monitor count
854        uscale = total_mon/defmon               //unscaling factor
855        total_det = uscale*realsread[2]         //unscaled detector count
856        total_trn = uscale*realsread[39]        //unscaled trans det count
857        total_numruns = integersread[3] //number of runs in workfile
858        total_rtime = integersread[2]           //total counting time in workfile
859        //retrieve workfile beamcenter
860        wrk_beamx = realsread[16]
861        wrk_beamy = realsread[17]
862        //unscale the workfile data in "newType"
863        //
864        //check for log-scaling and adjust if necessary
865        // should not be needed now - using display flag instead
866//      ConvertFolderToLinearScale(newType)
867        //
868        //then unscale the data array
869        data *= uscale
870        dest_data_err *= uscale
871       
872        //DetCorr() has not been applied to the data in RAW , do it now in a local reference to the raw data
873        WAVE raw_data = $"root:Packages:NIST:RAW:linear_data"
874        WAVE raw_data_err = $"root:Packages:NIST:RAW:linear_data_error"
875        WAVE raw_reals =  $"root:Packages:NIST:RAW:realsread"
876        WAVE/T raw_text = $"root:Packages:NIST:RAW:textread"
877        WAVE raw_ints = $"root:Packages:NIST:RAW:integersread"
878       
879        //check for log-scaling of the raw data - make sure it's linear
880        // should not be needed now - using display flag instead
881//      ConvertFolderToLinearScale("RAW")
882       
883        // switches to control what is done, don't do the transmission correction for the BGD measurement
884        NVAR doEfficiency = root:Packages:NIST:gDoDetectorEffCor
885        NVAR gDoTrans = root:Packages:NIST:gDoTransmissionCor
886        Variable doTrans = gDoTrans
887        if(cmpstr("BGD",newtype) == 0)
888                doTrans = 0             //skip the trans correction for the BGD file but don't change the value of the global
889        endif   
890       
891        V_DetCorr(raw_data,raw_data_err,raw_reals,doEfficiency,doTrans) //applies correction to raw_data, and overwrites it
892       
893        //deadtime corrections to raw data
894        // TODO - do the tube correction for dead time now
895        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
896        itim = raw_ints[2]
897        cntrate = sum(raw_data,-inf,inf)/itim           //080802 use data sum, rather than scaler value
898        dscale = 1/(1-deadTime*cntrate)
899
900#if (exists("ILL_D22")==6)
901        Variable tubeSum
902        // for D22 detector might need to use cntrate/128 as it is the tube response
903        for(ii=0;ii<pixelsX;ii+=1)
904                //sum the counts in each tube
905                tubeSum = 0
906                for(jj=0;jj<pixelsY;jj+=1)
907                        tubeSum += data[jj][ii]
908                endfor
909                // countrate in tube ii
910                cntrate = tubeSum/itim
911                // deadtime scaling in tube ii
912                dscale = 1/(1-deadTime*cntrate)
913                // multiply data[ii][] by the dead time
914                raw_data[][ii] *= dscale
915                raw_data_err[][ii] *= dscale
916        endfor
917#else
918        // dead time correction on all other RAW data, including NCNR
919        raw_data *= dscale
920        raw_data_err *= dscale
921#endif
922
923        //update totals by adding RAW values to the local ones (write to work header at end of function)
924        total_mon += raw_reals[0]
925
926        total_det += dscale*raw_reals[2]
927
928        total_trn += raw_reals[39]
929        total_rtime += raw_ints[2]
930        total_numruns +=1
931       
932        //do the beamcenter shifting if there is a mismatch
933        //and then add the two data sets together, changing "data" since it is the workfile data
934        xshift = raw_reals[16] - wrk_beamx
935        yshift = raw_reals[17] - wrk_beamy
936       
937        If((xshift != 0) || (yshift != 0))
938                DoAlert 1,"Do you want to ignore the beam center mismatch?"
939                if(V_flag==1)
940                        xshift=0
941                        yshift=0
942                endif
943        endif
944       
945        If((xshift == 0) && (yshift == 0))              //no shift, just add them
946                data += raw_data                //deadtime correction has already been done to the raw data
947                dest_data_err = sqrt(dest_data_err^2 + raw_data_err^2)                  // error of the sum
948        Endif
949       
950        //scale the data to the default montor counts
951        scale = defmon/total_mon
952        data *= scale
953        dest_data_err *= scale
954       
955        // keep "data" and linear_data in sync in the destination folder
956        data_copy = data
957       
958        //all is done, except for the bookkeeping of updating the header info in the work folder
959        textread[1] = date() + " " + time()             //date + time stamp
960        integersread[3] = total_numruns                                         //numruns = more than one
961        realsread[1] = total_mon                        //save the true monitor count
962        realsread[0] = defmon                                   //monitor ct = defmon
963        integersread[2] = total_rtime                   // total counting time
964        realsread[2] = scale*total_det                  //scaled detector counts
965        realsread[39] = scale*total_trn                 //scaled transmission counts
966       
967        //Add the added raw filename to the list of files in the workfile
968        String newfile = ";" + raw_text[0]
969        SVAR oldList = $(destPath + ":fileList")
970        String/G $(destPath + ":fileList") = oldList + newfile
971       
972        //reset the current display type to "newtype"
973        SVAR gCurDispType = root:Packages:NIST:VSANS:Globals:gCurDispType
974        gCurDispType = newType
975       
976        //return to root folder (redundant)
977        SetDataFolder root:
978       
979        Return(0)
980End
981
982
983//used for adding DRK (beam shutter CLOSED) data to a workfile
984//force the monitor count to 1, since it's irrelevant
985// run data through normal "add" step, then unscale default monitor counts
986//to get the data back on a simple time basis
987//
988Function V_Raw_to_Work_NoNorm(type)
989        String type
990       
991        WAVE reals=$("root:Packages:NIST:RAW:realsread")
992        reals[1]=1              //true monitor counts, still in raw
993        V_Raw_to_work(type)
994        //data is now in "type" folder
995        WAVE data=$("root:Packages:NIST:"+type+":linear_data")
996        WAVE data_copy=$("root:Packages:NIST:"+type+":data")
997        WAVE data_err=$("root:Packages:NIST:"+type+":linear_data_error")
998        WAVE new_reals=$("root:Packages:NIST:"+type+":realsread")
999       
1000        Variable norm_mon,tot_mon,scale
1001       
1002        norm_mon = new_reals[0]         //should be 1e8
1003        tot_mon = new_reals[1]          //should be 1
1004        scale= norm_mon/tot_mon
1005       
1006        data /= scale           //unscale the data
1007        data_err /= scale
1008       
1009        // to keep "data" and linear_data in sync
1010        data_copy = data
1011       
1012        return(0)
1013End
1014
1015//used for adding DRK (beam shutter CLOSED) data to a workfile
1016//force the monitor count to 1, since it's irrelevant
1017// run data through normal "add" step, then unscale default monitor counts
1018//to get the data back on a simple time basis
1019//
1020Function V_Add_Raw_to_Work_NoNorm(type)
1021        String type
1022       
1023        WAVE reals=$("root:Packages:NIST:RAW:realsread")
1024        reals[1]=1              //true monitor counts, still in raw
1025        V_Add_Raw_to_work(type)
1026        //data is now in "type" folder
1027        WAVE data=$("root:Packages:NIST:"+type+":linear_data")
1028        WAVE data_copy=$("root:Packages:NIST:"+type+":data")
1029        WAVE data_err=$("root:Packages:NIST:"+type+":linear_data_error")
1030        WAVE new_reals=$("root:Packages:NIST:"+type+":realsread")
1031       
1032        Variable norm_mon,tot_mon,scale
1033       
1034        norm_mon = new_reals[0]         //should be 1e8
1035        tot_mon = new_reals[1]          //should be equal to the number of runs (1 count per run)
1036        scale= norm_mon/tot_mon
1037       
1038        data /= scale           //unscale the data
1039        data_err /= scale
1040       
1041        // to keep "data" and linear_data in sync
1042        data_copy = data
1043       
1044        return(0)
1045End
1046
Note: See TracBrowser for help on using the repository browser.