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

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

updating the version number in VSANS, even though there were no changes there.

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