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

Last change on this file since 999 was 999, checked in by srkline, 6 years ago

changes to a few analysis models to make these Igor 7-ready

adding mask editing utilities

many changes to event mode for easier processing of split lists

updated event mode help file

+ lots more!

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