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

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

bug fixes for the back detector dimensions and beam center (pixels)

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_getMonitorCount(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.