source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/USANS/BT5_Loader.ipf @ 1122

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

changes to the USANS procedures to accommodate the re-ordering of the data columns in the raw USANS data that is to be output from NICE, versus what was previously written out by ICP.

A global flag switches between the two reading modes. On startup, the preferences panel is automatically opened to the USANS tab so that the user can immediately set the flag correctly. Currently it defaults to "checked" to read the New data format.

  • Property eol-style set to native
  • Property svn:executable set to *
File size: 12.1 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma Version=2.20
3#pragma IgorVersion=6.1
4
5/////////////////////////
6// 101001 Vers. 1
7//
8// functions to load and parse the ASCII ICP files generated at BT5 USANS
9// Nearly all of the important information about the data is dragged around
10// with the DetCts wave as a wave note
11//
12// - thes wave note is a string of KEY:value items
13//
14//      str = "FILE:"+filen+";"
15//      str += "TIMEPT:"+num2str(counttime)+";"
16//      str += "PEAKANG:;"              //no value yet
17//      str += "STATUS:;"               //no value yet
18//      str += "LABEL:"+fileLabel+";"
19//      str += "PEAKVAL:;"              //no value yet
20//      str += "TWIDE:0;"               //no value yet
21//
22/////////////////////////
23
24
25//given the filename, loads a single data file into a "TYPE" folder, into a
26//series of 1D waves. counting time is associated with the DetCts wave as a note
27//"blank" note entries are set up here as well, so later they can be polled/updated
28Function LoadBT5File(fname,type)
29        String fname,type
30
31        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
32       
33        Variable num=500,err=0,refnum
34        Make/O/D/N=(num) $(USANSFolder+":"+type+":Angle")
35        Make/O/D/N=(num) $(USANSFolder+":"+type+":DetCts")
36        Make/O/D/N=(num) $(USANSFolder+":"+type+":ErrDetCts")
37        Make/O/D/N=(num) $(USANSFolder+":"+type+":MonCts")
38        Make/O/D/N=(num) $(USANSFolder+":"+type+":TransCts")
39        Wave Angle = $(USANSFolder+":"+type+":Angle")
40        Wave DetCts = $(USANSFolder+":"+type+":DetCts")
41        Wave ErrDetCts = $(USANSFolder+":"+type+":ErrDetCts")
42        Wave MonCts = $(USANSFolder+":"+type+":MonCts")
43        Wave TransCts = $(USANSFolder+":"+type+":TransCts")
44       
45        Open/R refNum as fname          //if fname is "", a dialog will be presented
46        if(refnum==0)
47                return(1)               //user cancelled
48        endif
49        //read in the ASCII data line-by-line
50        Variable numLinesLoaded = 0,firstchar,countTime
51        Variable v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16,ii,valuesRead
52        String buffer ="",s1,s2,s3,s4,s5,s6,s7,s8,s9,s10
53        String filen="",fileLabel="",filedt=""
54       
55        Variable MainDeadTime,TransDeadTime
56       
57        //parse the first line
58        FReadLine refnum,buffer
59        sscanf buffer, "%s%s%s%s%s%s%g%g%s%g%s",s1,s2,s3,s4,s5,s6,v1,v2,s7,v3,s8
60        ii=strlen(s1)
61        filen=s1[1,(ii-2)]                      //remove the ' from beginning and end
62        //Print "filen = ",filen,strlen(filen)
63        //v2 is the monitor prefactor. multiply monitor time by prefactor. AJJ 5 March 07
64        countTime=v1*v2
65        //Print "time (sec) = ",countTime
66       
67        //Deadtime correction
68        //Discovered significant deadtime on detectors in Oct 2010
69        //JGB and AJJ changed pre-amps during shutdown Oct 27 - Nov 7 2010
70        //Need different deadtime before and after 8th November 2010
71        filedt = s2+" "+s3+" "+s4+" "+s5
72        ii=strlen(filedt)
73        filedt = filedt[1,(ii-2)]
74               
75//      print filedt
76//      print BT5Date2Secs(filedt)
77//      print date2secs(2010,11,7)
78       
79        // as of Feb 6 2019, USANS is using NICE for data collection rather than ICP
80        // as a consequence, the data file is somewhat different. Specifically, the
81        // order of the detectors in the comma-delimited list is different, with the
82        // 5 main detectors all listed in order, rather than separated.
83        // --- I can either flag on the date, or set a global to switch
84        // I also do not know if there will be any other changes in the data file format
85       
86        Variable useNewDataFormat
87
88// test by date
89//      Variable thisFileSecs,switchSecs
90//      switchSecs = date2secs(2019,2,6)
91//      thisFileSecs = BT5Date2Secs(filedt)             // could use BT5DateTime2Secs() to include HR:MIN
92//      if(thisFileSecs >= switchSecs)
93//              useNewDataFormat = 1
94//      else
95//              useNewDataFormat = 0
96//      endif
97
98// or use the global
99        NVAR gVal = root:Packages:NIST:gUseNICEDataFormat       
100        useNewDataFormat = gVal
101       
102       
103        USANS_DetectorDeadtime(filedt,MainDeadTime,TransDeadTime)
104       
105        //skip line 2
106        FReadLine refnum,buffer
107        //the next line is the sample label, use it all, minus the terminator
108        FReadLine refnum,filelabel
109       
110        //skip the next 10 lines
111        For(ii=0;ii<10;ii+=1)
112                FReadLine refnum,buffer
113        EndFor
114       
115        //read the data until EOF - assuming always a pair or lines
116        do
117                FReadLine refNum, buffer
118                if (strlen(buffer) == 0)
119                        break                                                   // End of file
120                endif
121                firstChar = char2num(buffer[0])
122                if ( (firstChar==10) || (firstChar==13) )
123                        break                                                   // Hit blank line. End of data in the file.
124                endif
125                sscanf buffer,"%g%g%g%g%g",v1,v2,v3,v4,v5               // 5 values here now
126                angle[numlinesloaded] = v1              //[0] is the ANGLE
127               
128                FReadLine refNum,buffer //assume a 2nd line is there, w/16 values
129                sscanf buffer,"%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g",v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16
130                //valuesRead = V_flag
131                //print valuesRead
132                if(useNewDataFormat == 1)
133                // new order is MonCt, det1,det2,det3,det4,det5, Trans
134                        MonCts[numlinesloaded] = v1             //monitor
135                        v2 = v2/(1.0-v2*MainDeadTime/countTime)   // Deadtime correction
136                        v3 = v3/(1.0-v3*MainDeadTime/countTime)   // Deadtime correction
137                        v4 = v6/(1.0-v4*MainDeadTime/countTime)   // Deadtime correction
138                        v5 = v5/(1.0-v5*MainDeadTime/countTime)   // Deadtime correction
139                        v6 = v6/(1.0-v6*MainDeadTime/countTime)   // Deadtime correction
140                        DetCts[numlinesloaded] = v2 + v3 + v4 + v5 + v6         //5 detectors
141                        TransCts[numlinesloaded] = v7/(1.0-v7*TransDeadTime/countTime)          //trans detector+deadtime correction
142                        ErrDetCts[numlinesloaded] = sqrt(detCts[numlinesloaded])
143                        //values 8-16 are always zero
144                else
145                // this is the original format from ICP
146                        MonCts[numlinesloaded] = v1             //monitor
147                        v2 = v2/(1.0-v2*MainDeadTime/countTime)   // Deadtime correction
148                        v3 = v3/(1.0-v3*MainDeadTime/countTime)   // Deadtime correction
149                        v5 = v5/(1.0-v5*MainDeadTime/countTime)   // Deadtime correction
150                        v6 = v6/(1.0-v6*MainDeadTime/countTime)   // Deadtime correction
151                        v7 = v7/(1.0-v7*MainDeadTime/countTime)   // Deadtime correction
152                        DetCts[numlinesloaded] = v2 + v3 + v5 + v6 + v7         //5 detectors
153                        TransCts[numlinesloaded] = v4/(1.0-v4*TransDeadTime/countTime)          //trans detector+deadtime correction
154                        ErrDetCts[numlinesloaded] = sqrt(detCts[numlinesloaded])
155                        //values 8-16 are always zero
156                endif
157               
158                numlinesloaded += 1             //2 lines in file read, one "real" line of data
159        while(1)
160               
161        Close refNum            // Close the file.
162        //Print "numlines = ",numlinesloaded
163        //trim the waves to the correct number of points
164        Redimension/N=(numlinesloaded) Angle,MonCts,DetCts,TransCts,ErrDetCts
165       
166        //remove LF from end of filelabel
167        filelabel=fileLabel[0,(strlen(fileLabel)-2)]
168       
169        //set the wave note for the DetCts
170        String str=""
171        str = "FILE:"+filen+";"
172        str += "TIMEPT:"+num2str(counttime)+";"
173        str += "PEAKANG:;"              //no value yet
174        str += "STATUS:;"               //no value yet
175        str += "LABEL:"+fileLabel+";"
176        str += "PEAKVAL:;"              //no value yet
177        str += "TWIDE:0;"               //no value yet
178        Note DetCts,str
179       
180        String/G fileLbl = filelabel
181        return err                      // Zero signifies no error.     
182End
183
184
185//will search for peak in root:type:DetCts,  find and returns the angle
186// of the peak poistion (= q=0)
187//returns -9999 if error, appends value to wavenote
188// !can't return -1,0,1, since these may be the peak angle!
189//
190Function FindZeroAngle(type)
191        String type
192       
193        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder       
194       
195        Variable pkNotFound,pkPt,pkAngle,pkVal,temp
196        Wave angle = $(USANSFolder+":"+type+":Angle")
197        Wave detCts = $(USANSFolder+":"+type+":DetCts")
198
199
200        WaveStats/Q detcts
201        temp=V_Maxloc           //in points
202        FindPeak/P/Q/M=10/R=[(temp-10),(temp+10)] DetCts                //+/- 10 pts from maximum
203       
204        //FindPeak/P/Q/M=(1000) detCts          // /M=min ht -- peak must be at least 100x higher than 1st pt
205        pkNotFound=V_Flag               //V_Flag==1 if no pk found
206        pkPt = V_PeakLoc
207        pkVal = V_PeakVal               //for calc of T_rock
208        if(pkNotFound)
209                //DoAlert 0, "Peak not found"
210                Return (-9999)          //fatal error
211        Endif
212        pkAngle = Angle[pkPt]
213        //Print "Peak Angle = ",pkAngle
214        //update the note
215        String str=""
216        str=note(DetCts)
217        str = ReplaceNumberByKey("PEAKANG",str,pkangle,":",";")
218        str = ReplaceNumberByKey("PEAKVAL",str,pkVal,":",";")
219        Note/K DetCts
220        Note detCts,str
221       
222        return(pkAngle)
223End
224
225// given the peakAngle (q=0) and the "type" folder of scattering angle,
226// convert the angle to Q-values [degrees to (1/A)]
227// makes a new Qvals wave, duplicating the Angle wave, which is assumed to exist
228// Uses a conversion constant supplied by John Barker, and is hard-wired in
229//
230Function ConvertAngle2Qvals(type,pkAngle)
231        String type
232        Variable pkAngle
233
234        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
235       
236        Wave angle = $(USANSFolder+":"+type+":Angle")
237        Variable num=numpnts(angle)
238//      Variable deg2QConv=5.55e-5              //JGB -- 2/24/01
239        NVAR deg2QConv=root:Packages:NIST:USANS:Globals:MainPanel:deg2QConv
240
241       
242        Make/O/N=(num) $(USANSFolder+":"+type+":Qvals")
243        Wave qvals = $(USANSFolder+":"+type+":Qvals")   
244        Qvals = deg2QConv*(angle[p] - pkAngle)
245       
246        return(0)       //no error
247End
248
249//updates the wavenote with the average Trans det cts for calculation of T_Wide
250//  finds the average number of counts on the transmission detector at angles
251// greater than 2 deg. (per John, practical experience where the trans data levels off)
252//
253// error value of 1 is returned and wavenote not updated if level is not found
254//
255// now normalizes to the monitor counts
256//
257// 17 APR 2013 SRK - coverted the 2 degree "cutoff" to use q-values instead. Then the
258//  location is generic and can be used for either NCNR or KUSANS
259//
260Function FindTWideCts(type)
261        String type
262       
263        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
264       
265        Variable levNotFound,levPt,Cts,num,ii,sumMonCts
266       
267        Wave angle = $(USANSFolder+":"+type+":Angle")
268        Wave detCts = $(USANSFolder+":"+type+":DetCts")
269        Wave TransCts = $(USANSFolder+":"+type+":TransCts")
270        Wave monCts = $(USANSFolder+":"+type+":MonCts")
271        Wave Qvals = $(USANSFolder+":"+type+":Qvals")
272       
273        num=numpnts(TransCts)
274
275//      FindLevel/Q/P angle,2           //use angles greater than 2 deg
276//      Print "Using angle, pt = ",V_levelX
277       
278        FindLevel/Q/P Qvals,1e-4                //use angles greater than 2 deg = 2*5.55e-5 = 1e-4 (1/A)
279        Print "Using Qval, pt = ",V_levelX
280
281        levNotFound=V_Flag              //V_Flag==1 if no pk found
282       
283        if(levNotFound)
284                //post a warning, and use just the last 4 points...
285//              DoAlert 0,"You don't have "+type+" data past 2 degrees - so Twide may not be reliable"
286                DoAlert 0,"You don't have "+type+" data past 1e-4 A-1 (~2 degrees) - so Twide may not be reliable"
287                levPt = num-4
288        else
289                levPt = trunc(V_LevelX)         // in points, force to integer
290        endif
291       
292        //average the trans cts from levPt to the end of the dataset
293        Cts=0
294        sumMonCts=0
295       
296        for(ii=levPt;ii<num;ii+=1)
297                Cts += transCts[ii]
298                sumMonCts += monCts[ii]
299        endfor
300        // get the average over that number of data points
301        sumMonCts /= (num-levPt-1)
302        Cts /= (num-levPt-1)
303       
304//      Print "cts = ",cts
305//      Print "sumMoncts = ",sumMoncts
306       
307        //normalize to the average monitor counts
308        Cts /= sumMonCts
309       
310//      Print "normalized counts = ",cts
311
312
313        //update the note
314        Wave DetCts = $(USANSFolder+":"+type+":DetCts")
315        String str,strVal
316        str=note(DetCts)
317        str = ReplaceNumberByKey("TWIDE",str,Cts,":",";")
318        Note/K DetCts
319        Note detCts,str
320
321        return(0)
322End
323
324
325Function BT5DateTime2Secs(datestring)
326        String datestring
327       
328        Variable bt5secs
329       
330        String monthnums = "Jan:1;Feb:2;Mar:3;Apr:4;May:5;Jun:6;Jul:7;Aug:8;Sep:9;Oct:10;Nov:11;Dec:12"
331               
332        Variable bt5month = str2num(StringByKey(stringfromlist(0,datestring, " "),monthnums))
333        Variable bt5day = str2num(stringfromlist(1,datestring," "))
334        Variable bt5year = str2num(stringfromlist(2,datestring," "))
335        Variable bt5hours = str2num(stringfromlist(0,stringfromlist(3,datestring," "),":"))
336        Variable bt5mins = str2num(stringfromlist(1,stringfromlist(3,datestring," "),":"))
337       
338        bt5secs = date2secs(bt5year,bt5month,bt5day) + 3600*bt5hours + 60*bt5mins
339
340        return bt5secs
341End
342
343Function BT5Date2Secs(datestring)
344        String datestring
345       
346        Variable bt5secs
347       
348        String monthnums = "Jan:1;Feb:2;Mar:3;Apr:4;May:5;Jun:6;Jul:7;Aug:8;Sep:9;Oct:10;Nov:11;Dec:12"
349               
350        Variable bt5month = str2num(StringByKey(stringfromlist(0,datestring, " "),monthnums))
351        Variable bt5day = str2num(stringfromlist(1,datestring," "))
352        Variable bt5year = str2num(stringfromlist(2,datestring," "))
353       
354        bt5secs = date2secs(bt5year,bt5month,bt5day)
355
356        return bt5secs
357End
Note: See TracBrowser for help on using the repository browser.