source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/FACILITY_Utils.ipf @ 795

Last change on this file since 795 was 795, checked in by srkline, 11 years ago

Changes to SANS reduction that apply to other Facilities:

These changes are related to the propagation of errors in 2D, on a
per-pixel basis. These changes only affect the errors that are reported in
the QxQy? ASCII file output. The 1D code is unaffected.

If these changes are not implemented, then errors of zero will be substitued as defaults
for these experimental errors.

Upon data loading, an error matrix, linear_data_error is generated and filled with
error values appropriate for Poisson statistics (not simply sqrt(n)).

4 functions in FACILITY_DataReadWrite.ipf have been added, and they are rather
self-explanatory:

In FACILITY_Utils.ipf, the AttenuatorTransmission?() function now returns
an additional parameter, atten_err, which is one standard deviation of the
attenuator transmission value. It returns a default error=0 (which is
correct if no attenuation is used). Facilities can fill this function in
with their own estimates for the uncertainty in the attenutator transmission.

File size: 19.1 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma version=5.0
3#pragma IgorVersion=6.1
4
5// this file contains globals and functions that are specific to a
6// particular facility or data file format
7// branched out 29MAR07 - SRK
8//
9// functions are either labeled with the procedure file that calls them,
10// or noted that they are local to this file
11
12
13// initializes globals that are specific to a particular facility
14// - number of XY pixels
15// - pixexl resolution [cm]
16// - detector deadtime constant [s]
17//
18// called by Initialize.ipf
19//
20Function InitFacilityGlobals()
21
22        //Detector -specific globals
23        Variable/G root:myGlobals:gNPixelsX=128                                 // number of X and Y pixels
24        Variable/G root:myGlobals:gNPixelsY=128
25       
26        // pixel dimensions are now read directly from the file header.
27//      Variable/G root:myGlobals:PixelResDefault = 0.5                 //pixel resolution in cm
28       
29        Variable/G root:myGlobals:DeadtimeDefault = 3.4e-6              //deadtime in seconds
30
31        Variable/G root:myGlobals:apOff = 5.0           // (cm) distance from sample aperture to sample position
32
33End
34
35
36//**********************
37// Resolution calculation - used by the averaging routines
38// to calculate the resolution function at each q-value
39// - the return value is not used
40//
41// equivalent to John's routine on the VAX Q_SIGMA_AVE.FOR
42// Incorporates eqn. 3-15 from J. Appl. Cryst. (1995) v. 28 p105-114
43//
44// - 21 MAR 07 uses projected BS diameter on the detector
45// - APR 07 still need to add resolution with lenses. currently there is no flag in the
46//          raw data header to indicate the presence of lenses.
47//
48// - Aug 07 - added input to switch calculation based on lenses (==1 if in)
49//
50// - called by CircSectAvg.ipf and RectAnnulAvg.ipf
51//
52// passed values are read from RealsRead
53// except DDet and apOff, which are set from globals before passing
54//
55//
56Function/S getResolution(inQ,lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,SigmaQ,QBar,fSubS)
57        Variable inQ, lambda, lambdaWidth, DDet, apOff, S1, S2, L1, L2, BS, del_r,usingLenses
58        Variable &fSubS, &QBar, &SigmaQ         //these are the output quantities at the input Q value
59       
60        //lots of calculation variables
61        Variable a2, q_small, lp, v_lambda, v_b, v_d, vz, yg, v_g
62        Variable r0, delta, inc_gamma, fr, fv, rmd, v_r1, rm, v_r
63
64        //Constants
65        Variable vz_1 = 3.956e5         //velocity [cm/s] of 1 A neutron
66        Variable g = 981.0                              //gravity acceleration [cm/s^2]
67
68        String results
69        results ="Failure"
70
71        S1 *= 0.5*0.1                   //convert to radius and [cm]
72        S2 *= 0.5*0.1
73
74        L1 *= 100.0                     // [cm]
75        L1 -= apOff                             //correct the distance
76
77        L2 *= 100.0
78        L2 += apOff
79        del_r *= 0.1                            //width of annulus, convert mm to [cm]
80       
81        BS *= 0.5*0.1                   //nominal BS diameter passed in, convert to radius and [cm]
82        // 21 MAR 07 SRK - use the projected BS diameter, based on a point sample aperture
83        Variable LB
84        LB = 20.1 + 1.61*BS                     //distance in cm from beamstop to anode plane (empirical)
85        BS = bs + bs*lb/(l2-lb)         //adjusted diameter of shadow from parallax
86       
87        //Start resolution calculation
88        a2 = S1*L2/L1 + S2*(L1+L2)/L1
89        q_small = 2.0*Pi*(BS-a2)*(1.0-lambdaWidth)/(lambda*L2)
90        lp = 1.0/( 1.0/L1 + 1.0/L2)
91
92        v_lambda = lambdaWidth^2/6.0
93       
94//      if(usingLenses==1)                      //SRK 2007
95        if(usingLenses != 0)                    //SRK 2008 allows for the possibility of different numbers of lenses in header
96                v_b = 0.25*(S1*L2/L1)^2 +0.25*(2/3)*(lambdaWidth/lambda)^2*(S2*L2/lp)^2         //correction to 2nd term
97        else
98                v_b = 0.25*(S1*L2/L1)^2 +0.25*(S2*L2/lp)^2              //original form
99        endif
100       
101        v_d = (DDet/2.3548)^2 + del_r^2/12.0
102        vz = vz_1 / lambda
103        yg = 0.5*g*L2*(L1+L2)/vz^2
104        v_g = 2.0*(2.0*yg^2*v_lambda)                                   //factor of 2 correction, B. Hammouda, 2007
105
106        r0 = L2*tan(2.0*asin(lambda*inQ/(4.0*Pi) ))
107        delta = 0.5*(BS - r0)^2/v_d
108
109        if (r0 < BS)
110                inc_gamma=exp(gammln(1.5))*(1-gammp(1.5,delta))
111        else
112                inc_gamma=exp(gammln(1.5))*(1+gammp(1.5,delta))
113        endif
114
115        fSubS = 0.5*(1.0+erf( (r0-BS)/sqrt(2.0*v_d) ) )
116        if (fSubS <= 0.0)
117                fSubS = 1.e-10
118        endif
119        fr = 1.0 + sqrt(v_d)*exp(-1.0*delta) /(r0*fSubS*sqrt(2.0*Pi))
120        fv = inc_gamma/(fSubS*sqrt(Pi)) - r0^2*(fr-1.0)^2/v_d
121
122        rmd = fr*r0
123        v_r1 = v_b + fv*v_d +v_g
124
125        rm = rmd + 0.5*v_r1/rmd
126        v_r = v_r1 - 0.5*(v_r1/rmd)^2
127        if (v_r < 0.0)
128                v_r = 0.0
129        endif
130        QBar = (4.0*Pi/lambda)*sin(0.5*atan(rm/L2))
131        SigmaQ = QBar*sqrt(v_r/rmd^2 +v_lambda)
132
133        results = "success"
134        Return results
135End
136
137
138//Utility function that returns the detector resolution (in cm)
139//Global values are set in the Initialize procedure
140//
141//
142// - called by CircSectAvg.ipf, RectAnnulAvg.ipf, and ProtocolAsPanel.ipf
143//
144// fileStr is passed as TextRead[3]
145// detStr is passed as TextRead[9]
146//
147// *** as of Jan 2008, depricated. Now detector pixel sizes are read from the file header
148// rw[10] = x size (mm); rw[13] = y size (mm)
149//
150// depricated - pixel dimensions are read directly from the file header
151Function xDetectorPixelResolution(fileStr,detStr)
152        String fileStr,detStr
153       
154        Variable DDet
155
156        //your code here
157       
158        return(DDet)
159End
160
161//Utility function that returns the detector deadtime (in seconds)
162//Global values are set in the Initialize procedure
163//
164// - called by WorkFileUtils.ipf
165//
166// fileStr is passed as TextRead[3] and is the filename
167// detStr is passed as TextRead[9] and is an identifier for the detector
168//
169Function DetectorDeadtime(fileStr,detStr)
170        String fileStr,detStr
171       
172        Variable deadtime
173       
174// your code here
175
176        return(deadtime)
177End
178
179
180// item is a filename
181//
182// this function extracts some sort of number from the file
183// presumably some sort of automatically incrementing run number set by the
184// acquisition system
185//
186// this run number should be a unique identifier for the file
187//
188Function GetRunNumFromFile(item)
189        String item
190
191        Variable num=-1         // an invalid return value
192       
193        //your code here
194       
195        return (num)
196End
197
198// item is a filename
199//
200// this function extracts some sort of number from the file
201// presumably some sort of automatically incrementing run number set by the
202// acquisition system
203//
204// this run number should be a unique identifier for the file
205//
206// same as GetRunNumFromFile(0), just with a string return
207//
208// "ABC" returned as an invalid result
209Function/S GetRunNumStrFromFile(item)
210        String item
211       
212        String invalid = "ABC"  //"ABC" is not a valid run number, since it's text
213        String retStr
214        retStr=invalid
215       
216        //your code here
217       
218        return(retStr)
219End
220
221//returns a string containing the full path to the file containing the
222//run number "num". The null string is returned if no valid file can be found.
223//
224//
225// search in the path "catPathName" (hard-wired), will abort if this path does not exist
226//the file returned will be a RAW SANS data file, other types of files are
227//filtered out.
228//
229// called by Buttons.ipf and Transmission.ipf, and locally by parsing routines
230//
231Function/S FindFileFromRunNumber(num)
232        Variable num
233       
234        String fullName="",partialName="",item=""
235       
236        //make sure that path exists
237        PathInfo catPathName
238        String path = S_path
239        if (V_flag == 0)
240                Abort "folder path does not exist - use Pick Path button"
241        Endif
242
243        //your code here       
244
245        return(fullname)
246       
247End
248
249//function to test a file to see if it is a RAW SANS file
250//
251// returns truth 0/1
252//
253// called by many procedures (both external and local)
254//
255Function CheckIfRawData(fname)
256        String fname
257       
258
259//      if(your test here)
260//              //true, is raw data file
261//              Return(1)
262//      else
263//              //some other file
264//              Return(0)
265//      Endif
266
267End
268
269
270// function returns 1 if file is a transmission file, 0 if not
271//
272// called by Transmission.ipf, CatVSTable.ipf, NSORT.ipf
273//
274Function isTransFile(fName)
275        String fname
276       
277//      if(your test here)
278//              //yes, its a transmisison file
279//              Return(1)
280//      else
281//              //some other file
282//              Return(0)
283//      Endif
284End
285
286
287//function to remove all spaces from names when searching for filenames
288//the filename (as saved) will never have interior spaces (TTTTTnnn_AB _Bnnn)
289//but the text field in the header may
290//
291//returns a string identical to the original string, except with the interior spaces removed
292//
293// local function for file name manipulation
294//
295// no change needed here
296Function/S RemoveAllSpaces(str)
297        String str
298       
299        String tempstr = str
300        Variable ii,spc,len             //should never be more than 2 or 3 trailing spaces in a filename
301        ii=0
302        do
303                len = strlen(tempStr)
304                spc = strsearch(tempStr," ",0)          //is the last character a space?
305                if (spc == -1)
306                        break           //no more spaces found, get out
307                endif
308                str = tempstr
309                tempStr = str[0,(spc-1)] + str[(spc+1),(len-1)] //remove the space from the string
310        While(1)        //should never be more than 2 or 3
311       
312        If(strlen(tempStr) < 1)
313                tempStr = ""            //be sure to return a null string if problem found
314        Endif
315       
316        //Print strlen(tempstr)
317       
318        Return(tempStr)
319               
320End
321
322
323//Function attempts to find valid filename from partial name by checking for
324// the existence of the file on disk
325//
326// returns a valid filename (No path prepended) or a null string
327//
328// called by any functions, both external and local
329//
330Function/S FindValidFilename(partialName)
331        String PartialName
332       
333        String retStr=""
334       
335        //your code here
336       
337        return(retStr)
338End
339
340
341//returns a string containing filename (WITHOUT the ;vers)
342//the input string is a full path to the file (Mac-style, still works on Win in IGOR)
343//with the folders separated by colons
344//
345// called by MaskUtils.ipf, ProtocolAsPanel.ipf, WriteQIS.ipf
346//
347// NEEDS NO CHANGES
348//
349Function/S GetFileNameFromPathNoSemi(fullPath)
350        String fullPath
351       
352        Variable offset1,offset2
353        String filename=""
354        //String PartialPath
355        offset1 = 0
356        do
357                offset2 = StrSearch(fullPath, ":", offset1)
358                if (offset2 == -1)                              // no more colons ?
359                        fileName = FullPath[offset1,strlen(FullPath) ]
360                        //PartialPath = FullPath[0, offset1-1]
361                        break
362                endif
363                offset1 = offset2+1
364        while (1)
365       
366        //remove version number from name, if it's there - format should be: filename;N
367        filename =  StringFromList(0,filename,";")              //returns null if error
368       
369        Return filename
370End
371
372//returns a string containing filename (INCLUDING the ;vers)
373//the input string is a full path to the file (Mac-style, still works on Win in IGOR)
374//with the folders separated by colons
375//
376// local, currently unused
377//
378// NEEDS NO CHANGES
379//
380Function/S GetFileNameFromPathKeepSemi(fullPath)
381        String fullPath
382       
383        Variable offset1,offset2
384        String filename
385        //String PartialPath
386        offset1 = 0
387        do
388                offset2 = StrSearch(fullPath, ":", offset1)
389                if (offset2 == -1)                              // no more colons ?
390                        fileName = FullPath[offset1,strlen(FullPath) ]
391                        //PartialPath = FullPath[0, offset1-1]
392                        break
393                endif
394                offset1 = offset2+1
395        while (1)
396       
397        //keep version number from name, if it's there - format should be: filename;N
398       
399        Return filename
400End
401
402//given the full path and filename (fullPath), strips the data path
403//(Mac-style, separated by colons) and returns this path
404//this partial path is the same string that would be returned from PathInfo, for example
405//
406// - allows the user to save to a different path than catPathName
407//
408// called by WriteQIS.ipf
409//
410// NEEDS NO CHANGES
411//
412Function/S GetPathStrFromfullName(fullPath)
413        String fullPath
414       
415        Variable offset1,offset2
416        //String filename
417        String PartialPath
418        offset1 = 0
419        do
420                offset2 = StrSearch(fullPath, ":", offset1)
421                if (offset2 == -1)                              // no more colons ?
422                        //fileName = FullPath[offset1,strlen(FullPath) ]
423                        PartialPath = FullPath[0, offset1-1]
424                        break
425                endif
426                offset1 = offset2+1
427        while (1)
428       
429        Return PartialPath
430End
431
432//given the filename trim or modify the filename to get a new
433//file string that can be used for naming averaged 1-d files
434//
435// called by ProtocolAsPanel.ipf and Tile_2D.ipf
436//
437Function/S GetNameFromHeader(fullName)
438        String fullName
439        String newName = ""
440
441        //your code here
442       
443        Return(newName)
444End
445
446//list (input) is a list, typically returned from IndexedFile()
447//which is semicolon-delimited, and may contain filenames from the VAX
448//that contain version numbers, where the version number appears as a separate list item
449//(and also as a non-existent file)
450//these numbers must be purged from the list, especially for display in a popup
451//or list processing of filenames
452//the function returns the list, cleaned of version numbers (up to 11)
453//raw data files will typically never have a version number other than 1.
454//
455// if there are no version numbers in the list, the input list is returned
456//
457// called by CatVSTable.ipf, NSORT.ipf, Transmission.ipf, WorkFileUtils.ipf
458//
459//
460// NO CHANGE NEEDED
461//
462Function/S RemoveVersNumsFromList(list)
463        String list
464       
465        //get rid of version numbers first (up to 11)
466        Variable ii,num
467        String item
468        num = ItemsInList(list,";")
469        ii=1
470        do
471                item = num2str(ii)
472                list = RemoveFromList(item, list ,";" )
473                ii+=1
474        while(ii<12)
475       
476        return (list)
477End
478
479//input is a list of run numbers, and output is a list of filenames (not the full path)
480//*** input list must be COMMA delimited***
481//output is equivalent to selecting from the CAT table
482//if some or all of the list items are valid filenames, keep them...
483//if an error is encountered, notify of the offending element and return a null list
484//
485//output is COMMA delimited
486//
487// this routine is expecting that the "ask", "none" special cases are handled elsewhere
488//and not passed here
489//
490// called by Marquee.ipf, MultipleReduce.ipf, ProtocolAsPanel.ipf
491//
492// NO CHANGE NEEDED
493//
494Function/S ParseRunNumberList(list)
495        String list
496       
497        String newList="",item="",tempStr=""
498        Variable num,ii,runNum
499       
500        //expand number ranges, if any
501        list = ExpandNumRanges(list)
502       
503        num=itemsinlist(list,",")
504       
505        for(ii=0;ii<num;ii+=1)
506                //get the item
507                item = StringFromList(ii,list,",")
508                //is it already a valid filename?
509                tempStr=FindValidFilename(item) //returns filename if good, null if error
510                if(strlen(tempstr)!=0)
511                        //valid name, add to list
512                        //Print "it's a file"
513                        newList += tempStr + ","
514                else
515                        //not a valid name
516                        //is it a number?
517                        runNum=str2num(item)
518                        //print runnum
519                        if(numtype(runNum) != 0)
520                                //not a number -  maybe an error                       
521                                DoAlert 0,"List item "+item+" is not a valid run number or filename. Please enter a valid number or filename."
522                                return("")
523                        else
524                                //a run number or an error
525                                tempStr = GetFileNameFromPathNoSemi( FindFileFromRunNumber(runNum) )
526                                if(strlen(tempstr)==0)
527                                        //file not found, error
528                                        DoAlert 0,"List item "+item+" is not a valid run number. Please enter a valid number."
529                                        return("")
530                                else
531                                        newList += tempStr + ","
532                                endif
533                        endif
534                endif
535        endfor          //loop over all items in list
536       
537        return(newList)
538End
539
540//takes a comma delimited list that MAY contain number range, and
541//expands any range of run numbers into a comma-delimited list...
542//and returns the new list - if not a range, return unchanged
543//
544// local function
545//
546// NO CHANGE NEEDED
547//
548Function/S ExpandNumRanges(list)
549        String list
550       
551        String newList="",dash="-",item,str
552        Variable num,ii,hasDash
553       
554        num=itemsinlist(list,",")
555//      print num
556        for(ii=0;ii<num;ii+=1)
557                //get the item
558                item = StringFromList(ii,list,",")
559                //does it contain a dash?
560                hasDash = strsearch(item,dash,0)                //-1 if no dash found
561                if(hasDash == -1)
562                        //not a range, keep it in the list
563                        newList += item + ","
564                else
565                        //has a dash (so it's a range), expand (or add null)
566                        newList += ListFromDash(item)           
567                endif
568        endfor
569       
570        return newList
571End
572
573//be sure to add a trailing comma to the return string...
574//
575// local function
576//
577// NO CHANGE NEEDED
578//
579Function/S ListFromDash(item)
580        String item
581       
582        String numList="",loStr="",hiStr=""
583        Variable lo,hi,ii
584       
585        loStr=StringFromList(0,item,"-")        //treat the range as a list
586        hiStr=StringFromList(1,item,"-")
587        lo=str2num(loStr)
588        hi=str2num(hiStr)
589        if( (numtype(lo) != 0) || (numtype(hi) !=0 ) || (lo > hi) )
590                numList=""
591                return numList
592        endif
593        for(ii=lo;ii<=hi;ii+=1)
594                numList += num2str(ii) + ","
595        endfor
596       
597        Return numList
598End
599
600
601//returns the proper attenuation factor based on the instrument
602//
603// filestr is passed from TextRead[3] = the default directory, used to identify the instrument
604// lam is passed from RealsRead[26]
605// AttenNo is passed from ReaslRead[3]
606//
607// Attenuation factor as defined here is <= 1
608//
609// Facilities can pass ("",1,attenuationFactor) and have this function simply
610// spit back the attenuationFactor (that was read into rw[3])
611//
612// called by Correct.ipf, ProtocolAsPanel.ipf, Transmission.ipf
613// atten_err is one std. deviation, passed back by reference
614Function AttenuationFactor(fileStr,lam,attenNo,atten_err)
615        String fileStr
616        Variable lam,attenNo, &atten_err
617       
618        Variable attenFactor=1
619       
620        // your code here
621
622        return(attenFactor)
623End
624
625//function called by the popups to get a file list of data that can be sorted
626// this procedure simply removes the raw data files from the string - there
627//can be lots of other junk present, but this is very fast...
628//
629// could also use the alternate procedure of keeping only file with the proper extension
630//
631// another possibility is to get a listing of the text files, but is unreliable on
632// Windows, where the data file must be .txt (and possibly OSX)
633//
634// called by FIT_Ops.ipf, NSORT.ipf, PlotUtils.ipf
635//
636// modify for specific facilities by changing the "*.SA1*","*.SA2*","*.SA3*" stringmatch
637// items which are specific to NCNR
638//
639Function/S ReducedDataFileList(ctrlName)
640        String ctrlName
641
642        String list="",newList="",item=""
643        Variable num,ii
644       
645        //check for the path
646        PathInfo catPathName
647        if(V_Flag==0)
648                DoAlert 0, "Data path does not exist - pick the data path from the button on the main panel"
649                Return("")
650        Endif
651       
652        list = IndexedFile(catpathName,-1,"????")
653        num=ItemsInList(list,";")
654        //print "num = ",num
655        for(ii=(num-1);ii>=0;ii-=1)
656                item = StringFromList(ii, list  ,";")
657                //simply remove all that are not raw data files (SA1 SA2 SA3)
658                if( !stringmatch(item,"*.SA1*") && !stringmatch(item,"*.SA2*") && !stringmatch(item,"*.SA3*") )
659                        if( !stringmatch(item,".*") && !stringmatch(item,"*.pxp") && !stringmatch(item,"*.DIV"))                //eliminate mac "hidden" files, pxp, and div files
660                                newlist += item + ";"
661                        endif
662                endif
663        endfor
664        //remove VAX version numbers
665        newList = RemoveVersNumsFromList(newList)
666        //sort
667        newList = SortList(newList,";",0)
668
669        return newlist
670End
671
672// returns a list of raw data files in the catPathName directory on disk
673// - list is SEMICOLON-delimited
674//
675// called by PatchFiles.ipf, Tile_2D.ipf
676//
677Function/S GetRawDataFileList()
678       
679        //make sure that path exists
680        PathInfo catPathName
681        if (V_flag == 0)
682                Abort "Folder path does not exist - use Pick Path button on Main Panel"
683        Endif
684       
685        String list=""
686       
687        // your code here
688       
689        return(list)
690End
691
692//**********************
693// 2D resolution function calculation - in terms of X and Y
694//
695// based on notes from David Mildner, 2008
696//
697// the final NCNR version is located in NCNR_Utils.ipf
698//
699Function/S get2DResolution(inQ,phi,lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,r_dist,SigmaQX,SigmaQY,fSubS)
700        Variable inQ, phi,lambda, lambdaWidth, DDet, apOff, S1, S2, L1, L2, BS, del_r,usingLenses,r_dist
701        Variable &SigmaQX,&SigmaQY,&fSubS               //these are the output quantities at the input Q value
702       
703        return("Function Empty")
704End
705
706
707// Return the filename that represents the previous or next file.
708// Input is current filename and increment.
709// Increment should be -1 or 1
710// -1 => previous file
711// 1 => next file
712Function/S GetPrevNextRawFile(curfilename, prevnext)
713        String curfilename
714        Variable prevnext
715
716        String filename
717       
718        //get the run number
719        Variable num = GetRunNumFromFile(curfilename)
720               
721        //find the next specified file by number
722        fileName = FindFileFromRunNumber(num+prevnext)
723
724        if(cmpstr(fileName,"")==0)
725                //null return, do nothing
726                fileName = FindFileFromRunNumber(num)
727        Endif
728
729//      print "in FU "+filename
730
731        Return filename
732End
Note: See TracBrowser for help on using the repository browser.