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

Last change on this file since 1103 was 907, checked in by srkline, 10 years ago

Changed the calculation of the wide angle transmission to look for q-values greater than 1e-4 (1/A) rather than > 2 degrees so that it would be more generic, and can be used for KUSANS. The 2 degree point is where the transmission counts levels off (empirically) and I simply converted this to Q.

  • Property eol-style set to native
  • Property svn:executable set to *
File size: 10.4 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        USANS_DetectorDeadtime(filedt,MainDeadTime,TransDeadTime)
80       
81        //skip line 2
82        FReadLine refnum,buffer
83        //the next line is the sample label, use it all, minus the terminator
84        FReadLine refnum,filelabel
85       
86        //skip the next 10 lines
87        For(ii=0;ii<10;ii+=1)
88                FReadLine refnum,buffer
89        EndFor
90       
91        //read the data until EOF - assuming always a pair or lines
92        do
93                FReadLine refNum, buffer
94                if (strlen(buffer) == 0)
95                        break                                                   // End of file
96                endif
97                firstChar = char2num(buffer[0])
98                if ( (firstChar==10) || (firstChar==13) )
99                        break                                                   // Hit blank line. End of data in the file.
100                endif
101                sscanf buffer,"%g%g%g%g%g",v1,v2,v3,v4,v5               // 5 values here now
102                angle[numlinesloaded] = v1              //[0] is the ANGLE
103               
104                FReadLine refNum,buffer //assume a 2nd line is there, w/16 values
105                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
106                //valuesRead = V_flag
107                //print valuesRead
108                MonCts[numlinesloaded] = v1             //monitor
109                v2 = v2/(1.0-V2*MainDeadTime/countTime)   // Deadtime correction
110                v3 = v3/(1.0-V3*MainDeadTime/countTime)   // Deadtime correction
111                v5 = v5/(1.0-V5*MainDeadTime/countTime)   // Deadtime correction
112                v6 = v6/(1.0-V6*MainDeadTime/countTime)   // Deadtime correction
113                v7 = v7/(1.0-V7*MainDeadTime/countTime)   // Deadtime correction
114                DetCts[numlinesloaded] = v2 + v3 + v5 + v6 + v7         //5 detectors
115                TransCts[numlinesloaded] = v4/(1.0-V4*TransDeadTime/countTime)          //trans detector+deadtime correction
116                ErrDetCts[numlinesloaded] = sqrt(detCts[numlinesloaded])
117                //values 8-16 are always zero
118                numlinesloaded += 1             //2 lines in file read, one "real" line of data
119        while(1)
120               
121        Close refNum            // Close the file.
122        //Print "numlines = ",numlinesloaded
123        //trim the waves to the correct number of points
124        Redimension/N=(numlinesloaded) Angle,MonCts,DetCts,TransCts,ErrDetCts
125       
126        //remove LF from end of filelabel
127        filelabel=fileLabel[0,(strlen(fileLabel)-2)]
128       
129        //set the wave note for the DetCts
130        String str=""
131        str = "FILE:"+filen+";"
132        str += "TIMEPT:"+num2str(counttime)+";"
133        str += "PEAKANG:;"              //no value yet
134        str += "STATUS:;"               //no value yet
135        str += "LABEL:"+fileLabel+";"
136        str += "PEAKVAL:;"              //no value yet
137        str += "TWIDE:0;"               //no value yet
138        Note DetCts,str
139       
140        String/G fileLbl = filelabel
141        return err                      // Zero signifies no error.     
142End
143
144
145//will search for peak in root:type:DetCts,  find and returns the angle
146// of the peak poistion (= q=0)
147//returns -9999 if error, appends value to wavenote
148// !can't return -1,0,1, since these may be the peak angle!
149//
150Function FindZeroAngle(type)
151        String type
152       
153        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder       
154       
155        Variable pkNotFound,pkPt,pkAngle,pkVal,temp
156        Wave angle = $(USANSFolder+":"+type+":Angle")
157        Wave detCts = $(USANSFolder+":"+type+":DetCts")
158
159
160        WaveStats/Q detcts
161        temp=V_Maxloc           //in points
162        FindPeak/P/Q/M=10/R=[(temp-10),(temp+10)] DetCts                //+/- 10 pts from maximum
163       
164        //FindPeak/P/Q/M=(1000) detCts          // /M=min ht -- peak must be at least 100x higher than 1st pt
165        pkNotFound=V_Flag               //V_Flag==1 if no pk found
166        pkPt = V_PeakLoc
167        pkVal = V_PeakVal               //for calc of T_rock
168        if(pkNotFound)
169                //DoAlert 0, "Peak not found"
170                Return (-9999)          //fatal error
171        Endif
172        pkAngle = Angle[pkPt]
173        //Print "Peak Angle = ",pkAngle
174        //update the note
175        String str=""
176        str=note(DetCts)
177        str = ReplaceNumberByKey("PEAKANG",str,pkangle,":",";")
178        str = ReplaceNumberByKey("PEAKVAL",str,pkVal,":",";")
179        Note/K DetCts
180        Note detCts,str
181       
182        return(pkAngle)
183End
184
185// given the peakAngle (q=0) and the "type" folder of scattering angle,
186// convert the angle to Q-values [degrees to (1/A)]
187// makes a new Qvals wave, duplicating the Angle wave, which is assumed to exist
188// Uses a conversion constant supplied by John Barker, and is hard-wired in
189//
190Function ConvertAngle2Qvals(type,pkAngle)
191        String type
192        Variable pkAngle
193
194        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
195       
196        Wave angle = $(USANSFolder+":"+type+":Angle")
197        Variable num=numpnts(angle)
198//      Variable deg2QConv=5.55e-5              //JGB -- 2/24/01
199        NVAR deg2QConv=root:Packages:NIST:USANS:Globals:MainPanel:deg2QConv
200
201       
202        Make/O/N=(num) $(USANSFolder+":"+type+":Qvals")
203        Wave qvals = $(USANSFolder+":"+type+":Qvals")   
204        Qvals = deg2QConv*(angle[p] - pkAngle)
205       
206        return(0)       //no error
207End
208
209//updates the wavenote with the average Trans det cts for calculation of T_Wide
210//  finds the average number of counts on the transmission detector at angles
211// greater than 2 deg. (per John, practical experience where the trans data levels off)
212//
213// error value of 1 is returned and wavenote not updated if level is not found
214//
215// now normalizes to the monitor counts
216//
217// 17 APR 2013 SRK - coverted the 2 degree "cutoff" to use q-values instead. Then the
218//  location is generic and can be used for either NCNR or KUSANS
219//
220Function FindTWideCts(type)
221        String type
222       
223        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
224       
225        Variable levNotFound,levPt,Cts,num,ii,sumMonCts
226       
227        Wave angle = $(USANSFolder+":"+type+":Angle")
228        Wave detCts = $(USANSFolder+":"+type+":DetCts")
229        Wave TransCts = $(USANSFolder+":"+type+":TransCts")
230        Wave monCts = $(USANSFolder+":"+type+":MonCts")
231        Wave Qvals = $(USANSFolder+":"+type+":Qvals")
232       
233        num=numpnts(TransCts)
234
235//      FindLevel/Q/P angle,2           //use angles greater than 2 deg
236//      Print "Using angle, pt = ",V_levelX
237       
238        FindLevel/Q/P Qvals,1e-4                //use angles greater than 2 deg = 2*5.55e-5 = 1e-4 (1/A)
239        Print "Using Qval, pt = ",V_levelX
240
241        levNotFound=V_Flag              //V_Flag==1 if no pk found
242       
243        if(levNotFound)
244                //post a warning, and use just the last 4 points...
245//              DoAlert 0,"You don't have "+type+" data past 2 degrees - so Twide may not be reliable"
246                DoAlert 0,"You don't have "+type+" data past 1e-4 A-1 (~2 degrees) - so Twide may not be reliable"
247                levPt = num-4
248        else
249                levPt = trunc(V_LevelX)         // in points, force to integer
250        endif
251       
252        //average the trans cts from levPt to the end of the dataset
253        Cts=0
254        sumMonCts=0
255       
256        for(ii=levPt;ii<num;ii+=1)
257                Cts += transCts[ii]
258                sumMonCts += monCts[ii]
259        endfor
260        // get the average over that number of data points
261        sumMonCts /= (num-levPt-1)
262        Cts /= (num-levPt-1)
263       
264//      Print "cts = ",cts
265//      Print "sumMoncts = ",sumMoncts
266       
267        //normalize to the average monitor counts
268        Cts /= sumMonCts
269       
270//      Print "normalized counts = ",cts
271
272
273        //update the note
274        Wave DetCts = $(USANSFolder+":"+type+":DetCts")
275        String str,strVal
276        str=note(DetCts)
277        str = ReplaceNumberByKey("TWIDE",str,Cts,":",";")
278        Note/K DetCts
279        Note detCts,str
280
281        return(0)
282End
283
284
285Function BT5DateTime2Secs(datestring)
286        String datestring
287       
288        Variable bt5secs
289       
290        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"
291               
292        Variable bt5month = str2num(StringByKey(stringfromlist(0,datestring, " "),monthnums))
293        Variable bt5day = str2num(stringfromlist(1,datestring," "))
294        Variable bt5year = str2num(stringfromlist(2,datestring," "))
295        Variable bt5hours = str2num(stringfromlist(0,stringfromlist(3,datestring," "),":"))
296        Variable bt5mins = str2num(stringfromlist(1,stringfromlist(3,datestring," "),":"))
297       
298        bt5secs = date2secs(bt5year,bt5month,bt5day) + 3600*bt5hours + 60*bt5mins
299
300        return bt5secs
301End
302
303Function BT5Date2Secs(datestring)
304        String datestring
305       
306        Variable bt5secs
307       
308        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"
309               
310        Variable bt5month = str2num(StringByKey(stringfromlist(0,datestring, " "),monthnums))
311        Variable bt5day = str2num(stringfromlist(1,datestring," "))
312        Variable bt5year = str2num(stringfromlist(2,datestring," "))
313       
314        bt5secs = date2secs(bt5year,bt5month,bt5day)
315
316        return bt5secs
317End
Note: See TracBrowser for help on using the repository browser.