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

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

moved dead time correction to before the solid angle correction, so that the dead time would be correcting counts, not counts per solid angle

added a routine to kill all of the waves and folders possible, if the overall DF kill failed. This is to prevent stray folders and waves from being present if different data files are loaded - since different data blocks are present for say, 3He data, data with temperature logging, etc.
This kill routine is used every time, before raw data is loaded, DIV or MASK loaded, or data is converted to WORK.

changed the "Save I(q)" button on the data display panel to save as ITX format, since the data has not been processed, and data can more easily be used for trimming input.

picking protocols in dialogs now excludes/includes appropriate waves

menus are consolidated

Fixed bug in SANS macros where the DRK[] item in the protocol could be null, and force the read of a DRK file, even if it was not desired.

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