source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/ANSTO_DataReadWrite.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: 54.9 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma version=5.0
3#pragma IgorVersion = 6.1 //required to read and write files created with HDF 1.8 library
4
5//**************************
6//
7// Vers. 1.2 092101
8// Vers. 5.0 29MAR07 - branched from main reduction to split out facility
9//                     specific calls
10//
11// functions for reading raw data files from the VAX
12// - RAW data files are read into the RAW folder - integer data from the detector
13//   is decompressed and given the proper orientation
14// - header information is placed into real,integer, or text waves in the order they appear
15//   in the file header
16//
17// Work data (DIV File) is read into the DIV folder
18//
19//*****************************
20
21//simple, main entry procedure that will load a RAW sans data file (not a work file)
22//into the RAW dataFolder. It is up to the calling procedure to display the file
23//
24// called by MainPanel.ipf and ProtocolAsPanel.ipf
25//
26Function LoadRawSANSData(msgStr)
27        String msgStr
28
29        String filename=""
30
31        //each routine is responsible for checking the current (displayed) data folder
32        //selecting it, and returning to root when done
33        PathInfo/S catPathName          //should set the next dialog to the proper path...
34        //get the filename, then read it in
35        filename = PromptForPath(msgStr)                //in SANS_Utils.ipf
36        //check for cancel from dialog
37        if(strlen(filename)==0)
38                //user cancelled, abort
39                SetDataFolder root:
40                DoAlert 0, "No file selected, action aborted"
41                return(1)
42        Endif
43
44        ReadHeaderAndData(filename)     //this is the full Path+file
45       
46        Return(0)
47End
48
49
50//function that does the guts of reading the binary data file
51//fname is the full path:name;vers required to open the file
52//
53// The final root:Packages:NIST:RAW:data wave is the real
54//neutron counts and can be directly operated on
55//
56// header information is put into three waves: integersRead, realsRead, and textRead
57// logicals in the header are currently skipped, since they are no use in the
58// data reduction process.
59//
60// The output is the three R/T/I waves that are filled at least with minimal values
61// and the detector data loaded into an array named "data"
62//
63// see documentation for the expected information in each element of the R/T/I waves
64// and the minimum set of information. These waves can be increased in length so that
65// more information can be accessed as needed (propagating changes...)
66//
67// called by multiple .ipfs (when the file name is known)
68//
69//
70// THIS FUNCTION DOES NOT NEED TO BE MODIFIED. ONLY THE DATA ACCESSORS NEED TO BE CONSTRUCTED
71//
72Function ReadHeaderAndData(fname)
73        String fname
74        //this function is for reading in RAW data only, so it will always put data in RAW folder
75        String curPath = "root:Packages:NIST:RAW"
76        SetDataFolder curPath           //use the full path, so it will always work
77        variable/g root:Packages:NIST:RAW:gIsLogScale = 0               //initial state is linear, keep this in RAW folder
78       
79        Variable refNum,integer,realval
80        String sansfname,textstr
81       
82        Make/D/O/N=23 $"root:Packages:NIST:RAW:IntegersRead"
83        Make/D/O/N=52 $"root:Packages:NIST:RAW:RealsRead"
84        Make/O/T/N=11 $"root:Packages:NIST:RAW:TextRead"
85        Make/O/N=7 $"root:Packages:NIST:RAW:LogicalsRead"
86
87        Wave intw=$"root:Packages:NIST:RAW:IntegersRead"
88        Wave realw=$"root:Packages:NIST:RAW:RealsRead"
89        Wave/T textw=$"root:Packages:NIST:RAW:TextRead"
90        Wave logw=$"root:Packages:NIST:RAW:LogicalsRead"
91       
92        // FILL IN 3 ARRAYS WITH HEADER INFORMATION FOR LATER USE
93        // THESE ARE JUST THE MINIMALLY NECESSARY VALUES
94       
95        // filename as stored in the file header
96        textw[0]= ParseFilePath(0, fname, ":", 1, 0)   
97       
98        // date and time of collection
99        textw[1]= getFileCreationDate(fname)
100       
101        // run type identifier (this is a reader for RAW data)
102        textw[2]= "RAW"
103       
104        // user account identifier (currently used only for NCNR-specific operations)
105        textw[3]= ""
106
107        // sample label
108        textw[6]= getSampleLabel(fname)
109       
110        // identifier of detector type, useful for setting detector constants
111        //(currently used only for NCNR-specific operations)
112        textw[9]= ""
113
114        //total counting time in seconds
115        intw[2] = getCountTime(fname)
116       
117//      realw[1] = 0
118//      realw[6] = 0
119//      realw[7] = 0
120//      realw[8] = 0
121//      realw[9] = 0
122//      realw[19] = 0
123//      realw[22] = 0
124
125        // total monitor count
126        realw[0] = getMonitorCount(fname)
127
128       
129        // attenuator number (NCNR-specific, your stub returns 0)
130        // may also be used to hold attenuator transmission (< 1)
131        realw[3] = getAttenNumber(fname)
132       
133        // sample transmission
134        realw[4] = getSampleTrans(fname)
135       
136        //sample thickness (cm)
137        realw[5] = getSampleThickness(fname)
138       
139        // following 6 values are for non-linear spatial corrections to a detector (RC timing)
140        // these values are set for a detctor that has a linear correspondence between
141        // pixel number and real-space distance
142        // 10 and 13 are the X and Y pixel dimensions, respectively (in mm!)
143        //(11,12 and 13,14 are set to values for a linear response, as from a new Ordela detector)
144        realw[10] = getDetectorPixelXSize(fname)
145        realw[11] = 10000
146        realw[12] = 0
147        realw[13] = getDetectorPixelYSize(fname)
148        realw[14] = 10000
149        realw[15] = 0
150       
151        // beam center X,Y on the detector (in units of pixel coordinates (1,N))
152        realw[16] = getBeamXPos(fname)
153        realw[17] = getBeamYPos(fname)
154       
155        // sample to detector distance (meters)
156        realw[18] = getSDD(fname)
157
158        // detector physical width (right now assumes square...) (in cm)
159        //realw[20] = 65
160        realw[20] = getPhysicalWidth(fname)
161       
162        // beam stop diameter (assumes circular) (in mm)
163        realw[21] = getBSDiameter(fname)
164       
165        // source aperture diameter (mm)
166        realw[23] = getSourceApertureDiam(fname)
167       
168        // sample aperture diameter (mm)
169        realw[24] = getSampleApertureDiam(fname)
170       
171        // source aperture to sample aperture distance
172        realw[25] = getSourceToSampleDist(fname)
173       
174        // wavelength (A)
175        realw[26] = getWavelength(fname)
176       
177        // wavelength spread (FWHM)
178        realw[27] = getWavelengthSpread(fname)
179       
180        // beam stop X-position (motor reading, approximate cm from zero position)
181        // currently NCNR-specific use to identify transmission measurements
182        // you return 0
183        realw[37] = 0
184
185// the actual data array, dimensions are set as globals in
186// InitFacilityGlobals()
187        NVAR XPix = root:myGlobals:gNPixelsX
188        NVAR YPix = root:myGlobals:gNPixelsX
189       
190        Make/D/O/N=(XPix,YPix) $"root:Packages:NIST:RAW:data"
191        WAVE data=$"root:Packages:NIST:RAW:data"
192
193        // fill the data array with the detector values
194        getDetectorData(fname,data)
195       
196        // total detector count
197        //nha 21/5/10 moved here because it requires the detector data to already be written
198        //Result of issue with 0 counts being written for a while in metadata.
199        realw[2] = getDetCount(fname)
200       
201        //keep a string with the filename in the RAW folder
202        String/G root:Packages:NIST:RAW:fileList = textw[0]
203       
204        Duplicate/O data linear_data                    //data read in is on linear scale, copy it now
205
206        // proper error for counting statistics, good for low count values too
207        // rather than just sqrt(n)
208        // see N. Gehrels, Astrophys. J., 303 (1986) 336-346, equation (7)
209        // for S = 1 in eq (7), this corresponds to one sigma error bars
210        Duplicate/O linear_data linear_data_error
211        linear_data_error = 1 + sqrt(linear_data + 0.75)                               
212        //
213       
214        Return 0
215
216End
217
218//****************
219//main entry procedure for reading a "WORK.DIV" file
220//displays a quick image of the  file, to check that it's correct
221//data is deposited in root:Packages:NIST:DIV data folder
222//
223// local, just for testing
224//
225Proc ReadWork_DIV()
226       
227        String fname = PromptForPath("Select detector sensitivity file")
228        ReadHeaderAndWork("DIV",fname)          //puts what is read in work.div
229       
230        String waveStr = "root:Packages:NIST:DIV:data"
231        NewImage/F/K=1/S=2 $waveStr
232        ModifyImage '' ctab= {*,*,YellowHot,0}
233       
234        String/G root:Packages:NIST:DIV:fileList = "WORK.DIV"
235       
236        SetDataFolder root:             //(redundant)
237End
238
239
240
241// Detector sensitivity files have the same header structure as RAW SANS data
242// as NCNR, just with a different data array (real, rather than integer data)
243//
244// So for your facility, make this reader specific to the format of whatever
245// file you use for a pixel-by-pixel division of the raw detector data
246// to correct for non-uniform sensitivities of the detector. This is typically a
247// measurement of water, plexiglas, or another uniform scattering sample.
248//
249// The data must be normalized to a mean value of 1
250//
251// called from ProtocolAsPanel.ipf
252//
253// type is "DIV" on input
254Function ReadHeaderAndWork(type,fname)
255        String type,fname
256       
257        //type is the desired folder to read the workfile to
258        //this data will NOT be automatically displayed
259        // gDataDisplayType is unchanged
260
261        String cur_folder = type
262        String curPath = "root:Packages:NIST:"+cur_folder
263        SetDataFolder curPath           //use the full path, so it will always work
264       
265        Variable refNum,integer,realval
266        String sansfname,textstr
267        Variable/G $(curPath + ":gIsLogScale") = 0              //initial state is linear, keep this in DIV folder
268       
269        Make/O/D/N=23 $(curPath + ":IntegersRead")
270        Make/O/D/N=52 $(curPath + ":RealsRead")
271        Make/O/T/N=11 $(curPath + ":TextRead")
272       
273        WAVE intw=$(curPath + ":IntegersRead")
274        WAVE realw=$(curPath + ":RealsRead")
275        WAVE/T textw=$(curPath + ":TextRead")
276       
277        // the actual data array, dimensions are set as globals in
278        // InitFacilityGlobals()
279        NVAR XPix = root:myGlobals:gNPixelsX
280        NVAR YPix = root:myGlobals:gNPixelsX
281
282        Make/O/D/N=(XPix,YPix) $(curPath + ":data")
283        WAVE data = $(curPath + ":data")
284       
285        // (1)
286        // fill in your reader for the header here so that intw, realw, and textW are filled in
287        // ? possibly a duplication of the raw data reader
288               
289        //(2)
290        // then fill in a reader for the data array that will DIVIDE your data
291        // to get the corrected values.
292        String dfName=""
293        variable err
294        err = hdfRead(fname, dfName)
295        Wave tempData = $dfName+":data:div"
296        data = tempData
297       
298        //funky edge correction bodgy ???
299        //copy second column to first column
300        data[][0] = data[p][1]
301        //copy second last column to last column
302        data[][191] = data[p][190]
303        //copy second row to first row
304        data[0][] = data[1][q]
305        //copy second last row to last row
306        data[191][] = data[190][q]
307        //keep a string with the filename in the DIV folder
308        String/G $(curPath + ":fileList") = textw[0]
309       
310        Duplicate/O data linear_data                    //data read in is on linear scale, copy it now
311
312        //return the data folder to root
313        SetDataFolder root:
314       
315        Return(0)
316End
317
318Function WriteHeaderAndWork(type)
319        String type
320       
321        // your writer here
322        NVAR XPix = root:myGlobals:gNPixelsX
323        NVAR YPix = root:myGlobals:gNPixelsX   
324       
325        Wave wData=$("root:Packages:NIST:"+type+":data")
326       
327//      Variable refnum,ii=0,hdrBytes=516,a,b,offset
328        String fname=""
329//      Duplicate/O wData,tempData
330       
331        //changed for Quokka
332//      Redimension/S/N=(XPix*YPix) tempData
333//      tempData *= 4
334       
335        PathInfo/S catPathName
336        fname = DoSaveFileDialog("Save data as")          //won't actually open the file
337        If(cmpstr(fname,"")==0)
338                //user cancel, don't write out a file
339          Close/A
340          Abort "no data file was written"
341        Endif
342       
343        variable fileID
344        HDF5CreateFile /O fileID as fname
345       
346        //nha ??? Should make this wave in our own DataFolder to avoid clashing names.
347        //Make /N=(1,1) wTransmission
348        String groupName = "/reduce"
349        String varName = "div"
350        // your code returning value
351        variable err
352        err = hdfWrite(fname, groupName, varName, wData)
353
354        // KillWaves wData, tempData
355        return(0)
356End
357
358
359
360/////   ASC FORMAT READER  //////
361/////   FOR WORKFILE MATH PANEL //////
362//
363//function to read in the ASC output of SANS reduction
364// currently the file has 20 header lines, followed by a single column
365// of 16384 values, Data is written by row, starting with Y=1 and X=(1->128)
366//
367//returns 0 if read was ok
368//returns 1 if there was an error
369//
370// called by WorkFileUtils.ipf
371//
372//
373// If the ASC data was generated by the NCNR data writer, then
374// NOTHING NEEDS TO BE CHANGED HERE
375//
376Function ReadASCData(fname,destPath)
377        String fname, destPath
378        //this function is for reading in ASCII data so put data in user-specified folder
379        SetDataFolder "root:Packages:NIST:"+destPath
380
381        NVAR pixelsX = root:myGlobals:gNPixelsX
382        NVAR pixelsY = root:myGlobals:gNPixelsY
383        Variable refNum=0,ii,p1,p2,tot,num=pixelsX,numHdrLines=20
384        String str=""
385        //data is initially linear scale
386        Variable/G :gIsLogScale=0
387        Make/O/T/N=(numHdrLines) hdrLines
388        Make/O/D/N=(pixelsX*pixelsY) data                       //,linear_data
389       
390        //full filename and path is now passed in...
391        //actually open the file
392//      SetDataFolder destPath
393        Open/R/Z refNum as fname                // /Z flag means I must handle open errors
394        if(refnum==0)           //FNF error, get out
395                DoAlert 0,"Could not find file: "+fname
396                Close/A
397                SetDataFolder root:
398                return(1)
399        endif
400        if(V_flag!=0)
401                DoAlert 0,"File open error: V_flag="+num2Str(V_Flag)
402                Close/A
403                SetDataFolder root:
404                return(1)
405        Endif
406        //
407        for(ii=0;ii<numHdrLines;ii+=1)          //read (or skip) 18 header lines
408                FReadLine refnum,str
409                hdrLines[ii]=str
410        endfor
411        //     
412        Close refnum
413       
414//      SetDataFolder destPath
415        LoadWave/Q/G/D/N=temp fName
416        Wave/Z temp0=temp0
417        data=temp0
418        Redimension/N=(pixelsX,pixelsY) data            //,linear_data
419
420        Duplicate/O data linear_data_error
421        linear_data_error = 1 + sqrt(data + 0.75)
422       
423        //just in case there are odd inputs to this, like negative intensities
424        WaveStats/Q linear_data_error
425        linear_data_error = numtype(linear_data_error[p]) == 0 ? linear_data_error[p] : V_avg
426        linear_data_error = linear_data_error[p] != 0 ? linear_data_error[p] : V_avg
427       
428        //linear_data = data
429       
430        KillWaves/Z temp0
431       
432        //return the data folder to root
433        SetDataFolder root:
434       
435        Return(0)
436End
437
438// fills the "default" fake header so that the SANS Reduction machinery does not have to be altered
439// pay attention to what is/not to be trusted due to "fake" information.
440// uses what it can from the header lines from the ASC file (hdrLines wave)
441//
442// destFolder is of the form "myGlobals:WorkMath:AAA"
443//
444//
445// called by WorkFileUtils.ipf
446//
447// If the ASC data was generated by the NCNR data writer, then
448// NOTHING NEEDS TO BE CHANGED HERE
449//
450Function FillFakeHeader_ASC(destFolder)
451        String destFolder
452        Make/O/D/N=23 $("root:Packages:NIST:"+destFolder+":IntegersRead")
453        Make/O/D/N=52 $("root:Packages:NIST:"+destFolder+":RealsRead")
454        Make/O/T/N=11 $("root:Packages:NIST:"+destFolder+":TextRead")
455       
456        Wave intw=$("root:Packages:NIST:"+destFolder+":IntegersRead")
457        Wave realw=$("root:Packages:NIST:"+destFolder+":RealsRead")
458        Wave/T textw=$("root:Packages:NIST:"+destFolder+":TextRead")
459       
460        //Put in appropriate "fake" values
461        //parse values as needed from headerLines
462        Wave/T hdr=$("root:Packages:NIST:"+destFolder+":hdrLines")
463        Variable monCt,lam,offset,sdd,trans,thick
464        Variable xCtr,yCtr,a1,a2,a1a2Dist,dlam,bsDiam
465        String detTyp=""
466        String tempStr="",formatStr="",junkStr=""
467        formatStr = "%g %g %g %g %g %g"
468        tempStr=hdr[3]
469        sscanf tempStr, formatStr, monCt,lam,offset,sdd,trans,thick
470//      Print monCt,lam,offset,sdd,trans,thick,avStr,step
471        formatStr = "%g %g %g %g %g %g %g %s"
472        tempStr=hdr[5]
473        sscanf tempStr,formatStr,xCtr,yCtr,a1,a2,a1a2Dist,dlam,bsDiam,detTyp
474//      Print xCtr,yCtr,a1,a2,a1a2Dist,dlam,bsDiam,detTyp
475       
476        realw[16]=xCtr          //xCtr(pixels)
477        realw[17]=yCtr  //yCtr (pixels)
478        realw[18]=sdd           //SDD (m)
479        realw[26]=lam           //wavelength (A)
480        //
481        // necessary values
482        realw[10]=5                     //detector calibration constants, needed for averaging
483        realw[11]=10000
484        realw[12]=0
485        realw[13]=5
486        realw[14]=10000
487        realw[15]=0
488        //
489        // used in the resolution calculation, ONLY here to keep the routine from crashing
490        realw[20]=65            //det size
491        realw[27]=dlam  //delta lambda
492        realw[21]=bsDiam        //BS size
493        realw[23]=a1            //A1
494        realw[24]=a2    //A2
495        realw[25]=a1a2Dist      //A1A2 distance
496        realw[4]=trans          //trans
497        realw[3]=0              //atten
498        realw[5]=thick          //thick
499        //
500        //
501        realw[0]=monCt          //def mon cts
502
503        // fake values to get valid deadtime and detector constants
504        //
505        textw[9]=detTyp+"  "            //6 characters 4+2 spaces
506        textw[3]="[NGxSANS00]"  //11 chars, NGx will return default values for atten trans, deadtime...
507       
508        //set the string values
509        formatStr="FILE: %s CREATED: %s"
510        sscanf hdr[0],formatStr,tempStr,junkStr
511//      Print tempStr
512//      Print junkStr
513        String/G $("root:Packages:NIST:"+destFolder+":fileList") = tempStr
514        textw[0] = tempStr              //filename
515        textw[1] = junkStr              //run date-time
516       
517        //file label = hdr[1]
518        tempStr = hdr[1]
519        tempStr = tempStr[0,strlen(tempStr)-2]          //clean off the last LF
520//      Print tempStr
521        textW[6] = tempStr      //sample label
522       
523        return(0)
524End
525
526
527///////// ACCESSORS FOR WRITING DATA TO HEADER  /////////
528/////
529
530// Write* routines all must:
531// (1) open file "fname". fname is a full file path and name to the file on disk
532// (2) write the specified value to the header at the correct location in the file
533// (3) close the file
534
535// new, April 2011 for error propagation. fill these in with the facility-
536// specific versions, if desired.
537Function WriteTransmissionErrorToHeader(fname,transErr)
538        String fname
539        Variable transErr
540       
541
542        return(0)
543End
544
545Function WriteBoxCountsErrorToHeader(fname,rel_err)
546        String fname
547        Variable rel_err
548       
549        return(0)
550End
551
552Function getSampleTransError(fname)
553        String fname
554       
555        return(0)
556end
557
558Function getBoxCountsError(fname)
559        String fname
560       
561        return(0)
562end
563
564
565// end April 2011 additions
566
567//whole transmission is NCNR-specific right now
568// leave this stub empty
569Function WriteWholeTransToHeader(fname,wholeTrans)
570        String fname
571        Variable wholeTrans
572       
573        String groupName = "/reduce"
574        variable err
575       
576        Wave wCounts
577        Make /N=(1,1) wWholeTrans
578               
579        wWholeTrans[0] = wholeTrans
580       
581        String varName = "wholeTransmission"   
582        err = hdfWrite(fname, groupName, varName,wWholeTrans)
583
584        KillWaves wWholeTrans
585       
586        //err not handled here
587        return(0)       
588End
589
590//box sum counts is a real value
591// used for transmission calculation module
592Function WriteBoxCountsToHeader(fname,counts)
593        String fname
594        Variable counts
595       
596        // do nothing if not using NCNR Transmission module
597       
598        String groupName = "/reduce"
599        variable err
600       
601        Wave wCounts
602        Make /N=(1,1) wCounts
603               
604        wCounts[0] = counts
605       
606        String varName = "boxCounts"   
607        err = hdfWrite(fname, groupName, varName,wCounts)
608
609        KillWaves wCounts
610       
611        //err not handled here
612        return(0)       
613End
614
615//beam stop X-position
616// used for transmission module to manually tag transmission files
617Function WriteBSXPosToHeader(fname,xpos)
618        String fname
619        Variable xpos
620       
621        // do nothing if not using NCNR Transmission module
622       
623        return(0)
624End
625
626
627
628
629
630//attenuator number (not its transmission)
631// if your beam attenuation is indexed in some way, use that number here
632// if not, write a 1 to the file here as a default
633//
634Function WriteAttenNumberToHeader(fname,attenNumber)
635        String fname
636        Variable attenNumber
637       
638        // your writer here
639        Wave wAttenNumber
640        //nha ??? Should make this wave in our own DataFolder to avoid clashing names.
641        Make /N=(1,1) wAttenNumber
642        String groupName = "/instrument/collimator"
643        String varName = "att"
644        //convert number to a rotation angle
645        attenNumber = attenNumber * 30
646        wAttenNumber[0] = attenNumber //
647        // your code returning value
648        variable err
649        err = hdfWrite(fname, groupName, varName, wAttenNumber)
650        KillWaves wAttenNumber
651        //err not handled here
652               
653        return(0)
654
655End
656
657
658// total monitor count during data collection
659Function WriteMonitorCountToHeader(fname,num)
660        String fname
661        Variable num
662       
663        // your code here
664       
665        return(0)
666End
667
668//total detector count
669Function WriteDetectorCountToHeader(fname,num)
670        String fname
671        Variable num
672       
673        // your code here
674       
675        return(0)
676End
677
678//transmission detector count
679// (currently unused in data reduction)
680Function WriteTransDetCountToHeader(fname,num)
681        String fname
682        Variable num
683       
684        // do nothing for now
685       
686        return(0)
687End
688
689
690
691//temperature of the sample (C)
692Function WriteTemperatureToHeader(fname,num)
693        String fname
694        Variable num
695       
696        //  your code here
697       
698        return(0)
699End
700
701//magnetic field (Oe)
702Function WriteMagnFieldToHeader(fname,num)
703        String fname
704        Variable num
705       
706        // your code here
707       
708        return(0)
709End
710
711//lateral detector offset (centimeters)
712Function WriteDetectorOffsetToHeader(fname,DetectorOffset)
713        String fname
714        Variable DetectorOffset
715       
716        // your writer here
717        Wave wDetectorOffset
718        //nha ??? Should make this wave in our own DataFolder to avoid clashing names.
719        Make /N=(1,1) wDetectorOffset
720        String groupName = "/instrument/detector"
721        String varName = "detector_x"
722        //convert from cm (NIST standard) to mm (NeXus standard)
723        DetectorOffset = DetectorOffset * 10
724        wDetectorOffset[0] = DetectorOffset //
725        // your code returning value
726        variable err
727        err = hdfWrite(fname, groupName, varName, wDetectorOffset)
728        KillWaves wDetectorOffset
729        //err not handled here
730       
731        return(0)
732End
733
734
735
736
737// physical dimension of detector x-pixel (mm)
738Function WriteDetPixelXToHeader(fname,num)
739        String fname
740        Variable num
741       
742        //your code here
743       
744        return(0)
745end
746
747// physical dimension of detector y-pixel (mm)
748Function WriteDetPixelYToHeader(fname,num)
749        String fname
750        Variable num
751       
752        //your code here
753       
754        return(0)
755end
756
757// sample label string
758Function WriteSamLabelToHeader(fname,str)
759        String fname,str
760       
761        // your code here
762
763        return(0)
764End
765
766// total counting time (seconds)
767Function WriteCountTimeToHeader(fname,num)
768        String fname
769        Variable num
770       
771        // your code here
772       
773        return(0)
774End
775
776
777
778//////// ACCESSORS FOR READING DATA FROM THE HEADER  //////////////
779//
780// read specific bits of information from the header
781// each of these operations MUST take care of open/close on their own
782//
783// fname is the full path and filname to the file on disk
784// return values are either strings or real values as appropriate
785//
786//////
787
788
789// function that reads in the 2D detector data and fills the array
790// data[nx][ny] with the data values
791// fname is the full name and path to the data file for opening and closing
792//
793//
794Function getDetectorData(fname,data)
795        String fname
796        Wave data
797        NVAR XPix = root:myGlobals:gNPixelsX
798       
799        // your reader here
800        variable err
801        string dfName = ""
802        err = hdfRead(fname, dfName)
803        if(err)
804                return 0
805        endif
806
807        Wave hmm_xy = $(dfName+":data:hmm_xy")
808       
809        //redimension /I /N = (dimsize(hmm_xy, 2), dimsize(hmm_xy, 1)), data
810        //nha. Count arrays need to be floating point, since the data will be divided, normalised etc.
811        redimension /N = (dimsize(hmm_xy, 2), dimsize(hmm_xy, 1)), data
812        data[][] = hmm_xy[0][q][p]
813       
814        //nha workaround. for wrongly dimensioned Quokka data 191x192
815        variable x_dim = dimsize(data,0)
816        if (x_dim!=XPix)
817                //redimension to add an extra row(s) to the data
818                redimension /I /N = (XPix,-1) data
819                //copy row 190 to row 191
820                data[191][] = data[190][q]
821        endif   
822        // end workaround
823       
824        //funky edge correction bodgy ???
825        //copy second column to first column
826        data[][0] = data[p][1]
827        //copy second last column to last column
828        data[][191] = data[p][190]
829        //copy second row to first row
830        data[0][] = data[1][q]
831        //copy second last row to last row
832        data[191][] = data[190][q]
833                       
834        KillWaves hmm_xy
835       
836       
837        return(0)
838End
839
840// file suffix (NCNR data file name specific)
841// return filename as suffix
842Function/S getSuffix(fname)
843        String fname
844       
845        return(ParseFilePath(0, fname, ":", 1, 0))
846End
847
848// associated file suffix (for transmission)
849// NCNR Transmission calculation specific
850// return null string
851Function/S getAssociatedFileSuffix(fname)
852        String fname
853       
854        return(getFileAssoc(fname))
855End
856
857// sample label
858Function/S getSampleLabel(fname)
859        String fname
860        String str
861       
862        // your code, returning str
863        variable err
864        string dfName = ""
865        err = hdfRead(fname, dfName)
866        //err not handled here
867       
868        Wave/T wSampleName = $dfname+":sample:name"
869        str = wSampleName[0]
870        KillWaves wSampleName
871       
872        return(str)
873End
874
875// file creation date
876Function/S getFileCreationDate(fname)
877        String fname
878        String str
879       
880        // your code, returning str
881        variable err
882        string dfName = ""
883        err = hdfRead(fname, dfName)
884        //err not handled here
885       
886        Wave/T wStartTime = $dfName+":start_time"
887        str = wStartTime[0]
888        KillWaves wStartTime
889       
890        return(str)
891End
892
893
894//monitor count
895Function getMonitorCount(fname)
896//not patched
897        String fname
898        Variable value
899       
900        // your code returning value
901        variable err
902        string dfName = ""
903        err = hdfRead(fname, dfName)
904        //err not handled here
905       
906        Wave wCounts = $dfName+":monitor:bm1_counts"
907        value = wCounts[0]
908        KillWaves wCounts
909       
910        return(value)
911end
912
913//saved monitor count
914// NCNR specific for pre-processed data, never used
915// return 0
916Function getSavMon(fname)
917        String fname
918       
919        Variable value
920       
921        // your code returning value
922        //??? to do. Is this required if 'never used'? nha
923       
924        return(0)
925end
926
927//total detector count
928Function getDetCount(fname)
929//not patched, but could be
930        String fname
931        Variable value
932       
933        // your code returning value
934        variable err
935        string dfName = ""
936        err = hdfRead(fname, dfName)
937        //err not handled here
938
939        //      Wave wCounts = $(dfName+":data:total_counts")
940        // changed 22/12/09 nha
941                if(WaveExists($(dfName+":data:total_counts")))
942                        Wave wCounts = $(dfName+":data:total_counts")
943                elseif(WaveExists($(dfName+":instrument:detector:total_counts")))
944                Wave wCounts = $(dfName+":instrument:detector:total_counts")
945               else
946                print "Can't find detector total_counts in " + fname
947               endif
948       
949        value = wCounts[0]
950       
951        //nha 21/5/10 temporary glitch wrote detector count to file as 0       
952                if (value<1)
953                        NVAR XPix = root:myGlobals:gNPixelsX
954                        NVAR YPix = root:myGlobals:gNPixelsX
955                        Make/D/O/N=(XPix,YPix) $"root:Packages:NIST:RAW:data"
956                        WAVE data=$"root:Packages:NIST:RAW:data"
957                        getDetectorData(fname,data)
958                        value = sum(data)
959                endif
960       
961        KillWaves wCounts
962       
963        return(value)
964end
965
966//Attenuator number, return 1 if your attenuators are not
967// indexed by number
968Function getAttenNumber(fname)
969        String fname
970        Variable value
971        Variable att, err
972        string dfName = ""
973               
974        err = hdfRead(fname, dfName)
975        //err not handled here
976
977        if(WaveExists($(dfName+":instrument:collimator:att")))
978                        Wave wAttrotdeg = $(dfName+":instrument:collimator:att")
979        elseif(WaveExists($(dfName+":instrument:parameters:derived_parameters:AttRotDeg")))
980                        Wave wAttrotdeg = $(dfName+":instrument:parameters:derived_parameters:AttRotDeg")
981        else
982                        print "Can't find attenuator in " + fname
983        endif   
984       
985        value = wAttrotdeg[0]
986        att = round(value)/30
987        KillWaves wAttrotdeg
988        return(att)
989end
990
991
992//transmission
993Function getSampleTrans(fname)
994        String fname
995       
996        Variable value
997       
998        // your code returning value
999        variable err
1000        string dfName = ""
1001        err = hdfRead(fname, dfName)
1002        //err not handled here
1003
1004        if(WaveExists($(dfName+":reduce:Transmission"))) //canonical location
1005                        Wave wTransmission = $(dfName+":reduce:Transmission")
1006        elseif(WaveExists($(dfName+":instrument:parameters:Transmission")))
1007                        Wave wTransmission = $(dfName+":instrument:parameters:Transmission")
1008        else
1009                        print "Can't find Transmission in " + fname
1010        endif
1011        value = wTransmission[0]
1012        KillWaves wTransmission
1013               
1014        return(value)
1015end
1016
1017//sample transmission (0<T<=1)
1018Function WriteTransmissionToHeader(fname,trans)
1019        String fname
1020        Variable trans
1021       
1022        // your writer here
1023        Wave wTransmission
1024        //nha ??? Should make this wave in our own DataFolder to avoid clashing names.
1025        Make /N=(1,1) wTransmission
1026        String groupName = "/reduce"
1027        String varName = "Transmission"
1028        wTransmission[0] = trans //
1029        // your code returning value
1030        variable err
1031        err = hdfWrite(fname, groupName, varName, wTransmission)
1032        KillWaves wTransmission
1033        //err not handled here
1034               
1035        return(0)
1036End
1037
1038//box counts from stored transmission calculation
1039// return 0 if not using NCNR transmission module
1040Function getBoxCounts(fname)
1041        String fname
1042        Variable value
1043       
1044        // do nothing if not using NCNR Transmission module
1045        variable err
1046        string dfName = ""
1047        err = hdfRead(fname, dfName)
1048        //err not handled here
1049
1050        Wave wBoxCounts = $(dfName+":reduce:boxCounts")
1051        if (waveexists(wBoxCounts) == 0)
1052                //boxcounts not yet set in  reduce group
1053                //return 0
1054                value = 0
1055        else
1056                value = wBoxCounts[0]
1057        endif
1058
1059        KillWaves/Z wBoxCounts
1060       
1061        return(value)
1062end
1063
1064//whole detector transmission
1065// return 0 if not using NCNR transmission module
1066Function getSampleTransWholeDetector(fname)
1067        String fname
1068        Variable value
1069       
1070        // your code returning value
1071        // ??? don't know what to put here. nha
1072        value=0
1073        return(value)
1074end
1075
1076//SampleThickness in centimeters
1077Function getSampleThickness(fname)
1078        String fname
1079        Variable value
1080
1081        // your code returning value
1082       
1083        variable err
1084        string dfName = ""
1085        err = hdfRead(fname, dfName)
1086        //err not handled here
1087
1088        if(WaveExists($(dfName+":sample:SampleThickness"))) //canonical location - a bit ugly and verbose, but that's just my opinion
1089                Wave wThickness = $(dfName+":sample:SampleThickness")
1090        elseif(WaveExists($(dfName+":sample:thickness")))
1091                Wave wThickness = $(dfName+":sample:thickness")
1092        else
1093                print "Can't find Sample Thickness in " + fname
1094        endif
1095                       
1096        value = wThickness[0]/10
1097        //value = 1 //??? temporary fix. nha
1098        KillWaves wThickness
1099       
1100        return(value)
1101end
1102
1103//sample thickness in cm
1104Function WriteThicknessToHeader(fname,thickness)
1105        String fname
1106        Variable thickness
1107       
1108        // your writer here
1109        Wave wThickness
1110        //nha ??? Should make this wave in our own DataFolder to avoid clashing names.
1111        Make /N=(1,1) wThickness
1112        String groupName = "/sample"
1113        String varName = "SampleThickness"
1114        wThickness[0] = thickness*10 //
1115        // your code returning value
1116        variable err
1117        err = hdfWrite(fname, groupName, varName, wThickness) //does not exist ???
1118        KillWaves wThickness
1119        //err not handled here
1120               
1121        return(0)
1122End
1123
1124//Sample Rotation Angle (degrees)
1125Function getSampleRotationAngle(fname)
1126        String fname
1127        Variable value
1128       
1129        // your code returning value
1130        variable err
1131        string dfName = ""
1132        err = hdfRead(fname, dfName)
1133        //err not handled here
1134
1135        Wave wSample_rotation_angle = $(dfName+":sample:sample_theta") //is this correct
1136        value = wSample_rotation_angle[0]
1137        KillWaves wSample_rotation_angle
1138               
1139        return(value)
1140end
1141
1142//temperature (C)
1143Function getTemperature(fname)
1144        String fname
1145       
1146        Variable value
1147       
1148        // your code returning value
1149       
1150        return(value)
1151end
1152
1153//field strength (Oe)
1154Function getFieldStrength(fname)
1155        String fname
1156       
1157        Variable value
1158       
1159        // your code returning value
1160       
1161        return(value)
1162end
1163
1164//center of beam xPos in pixel coordinates
1165Function getBeamXPos(fname)
1166        String fname
1167        Variable value
1168        // your code returning value
1169        variable err
1170        string dfName = ""
1171        err = hdfRead(fname, dfName)
1172        //err not handled here
1173
1174        if(WaveExists($(dfName+":instrument:reduce:BeamCenterX"))) //canonical location
1175                        Wave wBeamXPos = $(dfName+":instrument:reduce:BeamCenterX")
1176        elseif(WaveExists($(dfName+":instrument:parameters:BeamCenterX")))
1177                        Wave wBeamXPos = $(dfName+":instrument:parameters:BeamCenterX")
1178        else
1179                        print "Can't find BeamCenterX in" $fname
1180        endif
1181        value = wBeamXPos[0]   
1182        KillWaves wBeamXPos
1183       
1184        return(value)   
1185
1186end
1187
1188//beam center X pixel location (detector coordinates)
1189Function WriteBeamCenterXToHeader(fname,beamCenterX)
1190        String fname
1191        Variable beamCenterX
1192       
1193        // your writer here
1194        Wave wBeamCenterX
1195        //nha ??? Should make this wave in our own DataFolder to avoid clashing names.
1196        Make /N=(1,1) wBeamCenterX
1197        String groupName = "/instrument/reduce"
1198        String varName = "BeamCenterX"
1199        wBeamCenterX[0] = beamCenterX //
1200        // your code returning value
1201        variable err
1202        err = hdfWrite(fname, groupName, varName, wBeamCenterX)
1203        KillWaves wBeamCenterX
1204        //err not handled here
1205       
1206        return(0)
1207End
1208
1209//center of beam Y pos in pixel coordinates
1210Function getBeamYPos(fname)
1211        String fname
1212        Variable value
1213        // your code returning value
1214        variable err
1215        string dfName = ""
1216        err = hdfRead(fname, dfName)
1217        //err not handled here
1218
1219        if(WaveExists($(dfName+":instrument:reduce:BeamCenterZ"))) //canonical location
1220                        Wave wBeamYPos = $(dfName+":instrument:reduce:BeamCenterZ")
1221        elseif(WaveExists($(dfName+":instrument:parameters:BeamCenterZ")))
1222                        Wave wBeamYPos = $(dfName+":instrument:parameters:BeamCenterZ")
1223        else
1224                        print "Can't find BeamCenterZ in" $fname
1225        endif
1226        value = wBeamYPos[0]   
1227        KillWaves wBeamYPos
1228               
1229        return(value)
1230end
1231
1232//beam center Y pixel location (detector coordinates)
1233Function WriteBeamCenterYToHeader(fname,beamCenterY)
1234        String fname
1235        Variable beamCenterY
1236       
1237        // your writer here
1238        Wave wBeamCenterY
1239        //nha ??? Should make this wave in our own DataFolder to avoid clashing names.
1240        Make /N=(1,1) wBeamCenterY
1241        String groupName = "/instrument/reduce"
1242        String varName = "BeamCenterZ"
1243        wBeamCenterY[0] = beamCenterY //
1244        // your code returning value
1245        variable err
1246        err = hdfWrite(fname, groupName, varName, wBeamCenterY)
1247        KillWaves wBeamCenterY
1248        //err not handled here
1249
1250        return(0)
1251End
1252
1253
1254//sample to detector distance (meters)
1255Function getSDD(fname)
1256        String fname
1257        Variable value
1258        // your code returning value
1259        variable err
1260        string dfName = ""
1261        err = hdfRead(fname, dfName)
1262        //err not handled here
1263
1264        //workaround for bad HDF5 dataset
1265        if(WaveExists($(dfName+":instrument:parameters:L2"))) //canonical location
1266                Wave wSourceToDetectorDist = $(dfName+":instrument:parameters:L2")
1267        elseif(WaveExists($(dfName+":instrument:parameters:L2mm")))
1268                Wave wSourceToDetectorDist = $(dfName+":instrument:parameters:L2mm")
1269        elseif(WaveExists($(dfName+":instrument:parameters:derived_parameters:L2mm"))) 
1270                Wave wSourceToDetectorDist = $(dfName+":instrument:parameters:derived_parameters:L2mm")
1271        else
1272                print "Can't find L2 in " + fname
1273        endif
1274       
1275        value = wSourceToDetectorDist[0]/1000   
1276        KillWaves wSourceToDetectorDist
1277               
1278        return(value)
1279end
1280
1281//sample to detector distance (meters)
1282Function WriteSDDToHeader(fname,sdd)
1283        String fname
1284        Variable sdd
1285       
1286// your writer here
1287        Wave wSDD
1288        //nha ??? Should make this wave in our own DataFolder to avoid clashing names.
1289        Make /N=(1,1) wSDD
1290        String groupName = "/instrument/parameters"
1291        String varName = "L2"
1292        wSDD[0] = sdd * 1000 //
1293        // your code returning value
1294        variable err
1295        err = hdfWrite(fname, groupName, varName, wSDD)
1296        KillWaves wSDD
1297        //err not handled here
1298       
1299        return(0)
1300End
1301
1302//lateral detector offset (centimeters)
1303Function getDetectorOffset(fname)
1304        String fname
1305        Variable value
1306        // your code returning value
1307        variable err
1308        string dfName = ""
1309        err = hdfRead(fname, dfName)
1310        //err not handled here
1311
1312        Wave wDetectorOffset = $(dfName+":instrument:detector:detector_x") //is this correct
1313        value = wDetectorOffset[0]/10
1314        KillWaves wDetectorOffset
1315       
1316        return(value)
1317end
1318
1319//Beamstop diameter (millimeters)
1320Function getBSDiameter(fname)
1321        String fname
1322        Variable value
1323        // your code returning value
1324        variable err
1325        string dfName = ""
1326        err = hdfRead(fname, dfName)
1327        //err not handled here
1328
1329        if(WaveExists($(dfName+":instrument:parameters:BSdiam"))) //canonical location
1330                Wave wBSdiameter = $(dfName+":instrument:parameters:BSdiam")
1331        elseif(WaveExists($(dfName+":instrument:parameters:BSXmm")))
1332                Wave wBSdiameter = $(dfName+":instrument:parameters:BSXmm")
1333        else
1334                print "Can't find Beamstop Diameter in " + fname
1335        endif
1336        value = wBSdiameter[0]
1337        KillWaves wBSdiameter
1338       
1339        return(value)   
1340end
1341
1342//beam stop diameter (millimeters)
1343Function WriteBeamStopDiamToHeader(fname,bs)
1344        String fname
1345        Variable bs     
1346       
1347        // your writer here
1348        Wave wBS
1349        //nha ??? Should make this wave in our own DataFolder to avoid clashing names.
1350        Make /N=(1,1) wBS
1351        String groupName = "/instrument/parameters"
1352        String varName = "BSdiam"
1353        wBS[0] = bs //
1354        // your code returning value
1355        variable err
1356        err = hdfWrite(fname, groupName, varName, wBS)
1357        KillWaves wBS
1358        //err not handled here
1359        return(0)
1360End
1361
1362//source aperture diameter (millimeters)
1363Function getSourceApertureDiam(fname)
1364        String fname
1365        Variable value
1366        // your code returning value
1367        variable err
1368        string dfName = ""
1369        err = hdfRead(fname, dfName)
1370        //err not handled here
1371
1372        if(WaveExists($(dfName+":instrument:parameters:EApX")))
1373                Wave wSourceApertureDiam = $(dfName+":instrument:parameters:EApX") // canonical location
1374        elseif(WaveExists($(dfName+":instrument:parameters:derived_parameters:EApXmm")))
1375                Wave wSourceApertureDiam = $(dfName+":instrument:parameters:derived_parameters:EApXmm")
1376        else
1377                print "Can't find Source Aperture Diameter in " + fname
1378        endif   
1379        value = wSourceApertureDiam[0]
1380        KillWaves wSourceApertureDiam
1381       
1382        return(value)
1383end
1384
1385//Source Aperture diameter (millimeters)
1386Function WriteSourceApDiamToHeader(fname,source)
1387        String fname
1388        Variable source
1389       
1390        // your writer here
1391        Wave wsource
1392        //nha ??? Should make this wave in our own DataFolder to avoid clashing names.
1393        Make /N=(1,1) wsource
1394        String groupName = "/instrument/parameters"
1395        String varName = "EApX"
1396        wsource[0] = source //
1397        // your code returning value
1398        variable err
1399        err = hdfWrite(fname, groupName, varName, wsource)
1400        KillWaves wsource
1401        //err not handled here
1402        return(0)
1403End
1404
1405//sample aperture diameter (millimeters)
1406Function getSampleApertureDiam(fname)
1407        String fname
1408        Variable value
1409        // your code returning value
1410        variable err
1411        string dfName = ""
1412        err = hdfRead(fname, dfName)
1413        //err not handled here
1414
1415        if(WaveExists($(dfName+":sample:diameter"))) //canonical location
1416                Wave wSampleApertureDiam = $(dfName+":sample:diameter")
1417        elseif(WaveExists($(dfName+":instrument:parameters:autoSampleAp:diameter"))) //canonical location
1418                Wave wSampleApertureDiam = $(dfName+":instrument:parameters:autoSampleAp:diameter")
1419        elseif (WaveExists($(dfName+":instrument:sample_aperture:geometry:shape:SApXmm")))
1420                Wave wSampleApertureDiam = $(dfName+":instrument:sample_aperture:geometry:shape:SApXmm")
1421        else
1422                print "Can't find Sample Aperture Diameter in " + fname
1423        endif   
1424        value = wSampleApertureDiam[0]
1425        KillWaves wSampleApertureDiam
1426
1427        return(value)   
1428end
1429
1430//Sample Aperture diameter (millimeters)
1431Function WriteSampleApDiamToHeader(fname,source)
1432        String fname
1433        Variable source
1434       
1435        // your writer here
1436        Wave wsource
1437        //nha ??? Should make this wave in our own DataFolder to avoid clashing names.
1438        Make /N=(1,1) wsource
1439        String groupName = "/sample"
1440        String varName = "diameter"
1441        wsource[0] = source //
1442        // your code returning value
1443        variable err
1444        err = hdfWrite(fname, groupName, varName, wsource)
1445        KillWaves wsource
1446        //err not handled here
1447        return(0)
1448End
1449
1450//source AP to Sample AP distance (meters)
1451Function getSourceToSampleDist(fname)
1452        String fname
1453       
1454        Variable value
1455       
1456        // your code returning value
1457        variable err
1458        string dfName = ""
1459               
1460        err = hdfRead(fname, dfName)
1461        //err not handled here
1462       
1463        if(WaveExists($(dfName+":instrument:parameters:L1"))) //canonical location
1464                Wave wSourceToSampleDist = $(dfName+":instrument:parameters:L1")       
1465        elseif(WaveExists($(dfName+":instrument:parameters:L1mm")))
1466                Wave wSourceToSampleDist = $(dfName+":instrument:parameters:L1mm")
1467        elseif(WaveExists($(dfName+":instrument:parameters:derived_parameters:L1mm")))
1468                Wave wSourceToSampleDist = $(dfName+":instrument:parameters:derived_parameters:L1mm")
1469        else
1470                print "Can't find L1 in " + fname
1471        endif
1472       
1473        value = wSourceToSampleDist[0]/1000
1474        KillWaves wSourceToSampleDist   
1475               
1476        return(value)
1477end
1478
1479//Source aperture to sample aperture distance (meters)
1480Function WriteSrcToSamDistToHeader(fname,SSD)
1481        String fname
1482        Variable SSD
1483       
1484// your writer here
1485        Wave wSSD
1486        //nha ??? Should make this wave in our own DataFolder to avoid clashing names.
1487        Make /N=(1,1) wSSD
1488        String groupName = "/instrument/parameters"
1489        String varName = "L1"
1490        wSSD[0] = SSD * 1000 //input in metres, converted to mm for saving to file.
1491        // your code returning value
1492        variable err
1493        err = hdfWrite(fname, groupName, varName, wSSD)
1494        KillWaves wSSD
1495        //err not handled here
1496       
1497        return(0)
1498End
1499
1500//wavelength (Angstroms)
1501Function getWavelength(fname)
1502        String fname
1503        Variable value
1504        // your code returning value
1505        variable err
1506        string dfName = ""
1507        err = hdfRead(fname, dfName)
1508        //err not handled here
1509
1510        //      Wave wWavelength = $(dfName+":data:LambdaA")
1511        //change 22/12/09 nha
1512        // all these locations to be deprecated
1513        if(WaveExists($(dfName+":instrument:velocity_selector:Lambda")))  // canonical location
1514                Wave wWavelength = $(dfName+":instrument:velocity_selector:Lambda")
1515        elseif(WaveExists($(dfName+":data:Lambda")))
1516                Wave wWavelength = $(dfName+":data:Lambda")
1517        elseif(WaveExists($(dfName+":data:LambdaA")))
1518                Wave wWavelength = $(dfName+":data:LambdaA")
1519        elseif(WaveExists($(dfName+":instrument:velocity_selector:LambdaA")))
1520                Wave wWavelength = $(dfName+":instrument:velocity_selector:LambdaA")
1521        else
1522                print "Can't find Lambda in " + fname
1523        endif
1524        value = wWavelength[0] 
1525        KillWaves wWavelength
1526       
1527        return(value)
1528end
1529
1530//wavelength (Angstroms)
1531Function WriteWavelengthToHeader(fname,wavelength)
1532        String fname
1533        Variable wavelength
1534       
1535        // your writer here
1536        Wave wWavelength
1537        //nha ??? Should make this wave in our own DataFolder to avoid clashing names.
1538        Make /N=(1,1) wWavelength
1539        String groupName = "/instrument/velocity_selector"
1540        String varName = "Lambda"
1541        wWavelength[0] = wavelength //
1542        // your code returning value
1543        variable err
1544        err = hdfWrite(fname, groupName, varName, wWavelength)
1545       
1546        //and because Bill Hamilton is not happy with the NeXus naming, we write it to 3 other places
1547        //groupName = "/instrument/parameters"
1548        //err = hdfWrite(fname, groupName, varName, wWavelength)
1549        //velocity_selector group causes Igor crash in some cases
1550        //groupName = "/instrument/velocity_selector"
1551        //err = hdfWrite(fname, groupName, varName, wWavelength)
1552        //
1553        //groupName = "/data"
1554        //varName = "lambda"
1555        //err = hdfWrite(fname, groupName, varName, wWavelength)
1556       
1557        KillWaves wWavelength
1558        //err not handled here
1559
1560        return(0)
1561End
1562
1563
1564
1565//wavelength spread (FWHM)
1566Function getWavelengthSpread(fname)
1567        String fname
1568       
1569        Variable value
1570       
1571        // your code returning value
1572        variable err
1573        string dfName = ""
1574        err = hdfRead(fname, dfName)
1575        //err not handled here
1576       
1577        //velocity_selector group causes Igor crash
1578        if(WaveExists($(dfName+":instrument:velocity_selector:LambdaResFWHM_percent")))  //canonical location
1579                Wave wWavelengthSpread = $(dfName+":instrument:velocity_selector:LambdaResFWHM_percent")
1580        elseif(WaveExists($(dfName+":instrument:parameters:LambdaResFWHM_percent")))
1581                Wave wWavelengthSpread = $(dfName+":instrument:parameters:LambdaResFWHM_percent")
1582        else
1583                print "Can't find Wavelength Spread in " + fname
1584        endif
1585        value = wWavelengthSpread[0]   
1586        KillWaves wWavelengthSpread
1587       
1588        return(value)
1589end
1590
1591//wavelength spread (FWHM)
1592Function WriteWavelengthDistrToHeader(fname,wavelengthSpread)
1593        String fname
1594        Variable wavelengthSpread
1595       
1596        // your writer here
1597        Wave wWavelengthSpread
1598        //nha ??? Should make this wave in our own DataFolder to avoid clashing names.
1599        Make /N=(1,1) wWavelengthSpread
1600        //velocity_selector group causes Igor crash
1601        String groupName = "/instrument/velocity_selector"
1602        String varName = "LambdaResFWHM_percent"
1603
1604        wWavelengthSpread[0] = wavelengthSpread
1605        // your code returning value
1606        variable err
1607        err = hdfWrite(fname, groupName, varName, wWavelengthSpread)
1608        KillWaves wWavelengthSpread
1609        //err not handled here
1610               
1611        return(0)
1612End
1613
1614// physical x-dimension of a detector pixel, in mm
1615Function getDetectorPixelXSize(fname)
1616        String fname
1617        Variable value
1618       
1619        // your code here returning value
1620        variable err
1621        string dfName = ""
1622        err = hdfRead(fname, dfName)
1623        //err not handled here
1624
1625        Wave wActiveArea = $(dfName+":instrument:detector:active_height")
1626        Wave w_x_bin = $(dfName+":instrument:detector:x_bin")
1627        Variable numPixels = dimsize(w_x_bin, 0)
1628        value = wActiveArea[0]/numPixels
1629        KillWaves wActiveArea
1630        KillWaves w_x_bin
1631       
1632        return(value)
1633end
1634
1635// physical y-dimension of a detector pixel, in mm
1636Function getDetectorPixelYSize(fname)
1637        String fname
1638        Variable value
1639       
1640        // your code here returning value
1641        variable err
1642        string dfName = ""
1643        err = hdfRead(fname, dfName)
1644        //err not handled here
1645
1646        Wave wActiveArea = $(dfName+":instrument:detector:active_width")
1647        Wave w_y_bin = $(dfName+":instrument:detector:y_bin")
1648        Variable numPixels = dimsize(w_y_bin, 0)
1649        value = wActiveArea[0]/numPixels
1650        KillWaves wActiveArea
1651        KillWaves w_y_bin
1652                       
1653        return(value)
1654end
1655
1656//transmission detector count (unused, return 0)
1657//
1658Function getTransDetectorCounts(fname)
1659        String fname
1660       
1661        Variable value
1662       
1663        // your code returning value
1664       
1665        return(0)
1666end
1667
1668
1669//total count time (seconds)
1670Function getCountTime(fname)
1671        String fname
1672        Variable value
1673        // your code returning value
1674        variable err
1675        string dfName = ""
1676        err = hdfRead(fname, dfName)
1677        //err not handled here
1678
1679        Wave wTime1 = $(dfName+":monitor:bm1_time")
1680        value = wTime1[0]       
1681        KillWaves wTime1
1682       
1683        return(value)
1684end
1685
1686
1687Function getPhysicalWidth(fname)
1688        String fname
1689        Variable value
1690        // your code returning value
1691        variable err
1692        string dfName = ""
1693        err = hdfRead(fname, dfName)
1694        //err not handled here
1695
1696        Wave wPhysicalWidth = $(dfName+":instrument:detector:active_width")
1697        value = wPhysicalWidth[0]/10
1698        KillWaves wPhysicalWidth
1699               
1700        return(value)
1701end
1702
1703Function/S getSICSVersion(fname)
1704        String fname
1705        String value
1706        // your code returning value
1707        variable err
1708        string dfName = ""
1709        err = hdfRead(fname, dfName)
1710        //err not handled here
1711
1712        Wave/T wSICSVersion = $(dfName+":sics_release")
1713        value = wSICSVersion[0]
1714        KillWaves wSICSVersion
1715       
1716        return(value)
1717end
1718
1719Function/S getHDFversion(fname)
1720        String fname
1721        String value
1722        // your code returning value
1723        variable err
1724        string dfName = ""
1725        string attribute = "HDF5_Version"
1726        err = hdfReadAttribute(fname, dfName, "/", 1, attribute)
1727//      string attribute ="signal"
1728//      err = hdfReadAttribute(fname, dfName, "/entry/data/hmm_xy", 2, attribute)
1729        //err not handled here
1730
1731        Wave/T wHDF5_Version = $(dfName+":"+attribute)
1732        value = wHDF5_Version[0]       
1733//      KillWaves wHDF5_Version
1734
1735        return(value)
1736end
1737
1738
1739//reads the wavelength from a reduced data file (not very reliable)
1740// - does not work with NSORTed files
1741// - only used in FIT/RPA (which itself is almost NEVER used...)
1742//
1743// DOES NOT NEED TO BE CHANGED IF USING NCNR DATA WRITER
1744Function GetLambdaFromReducedData(tempName)
1745        String tempName
1746       
1747        String junkString
1748        Variable lambdaFromFile, fileVar
1749        lambdaFromFile = 6.0
1750        Open/R/P=catPathName fileVar as tempName
1751        FReadLine fileVar, junkString
1752        FReadLine fileVar, junkString
1753        FReadLine fileVar, junkString
1754        if(strsearch(LowerStr(junkString),"lambda",0) != -1)
1755                FReadLine/N=11 fileVar, junkString
1756                FReadLine/N=10 fileVar, junkString
1757                lambdaFromFile = str2num(junkString)
1758        endif
1759        Close fileVar
1760       
1761        return(lambdaFromFile)
1762End
1763
1764/////   TRANSMISSION RELATED FUNCTIONS    ////////
1765//box coordinate are returned by reference
1766//
1767// box to sum over is bounded (x1,y1) to (x2,y2)
1768//
1769// this function returns the bounding coordinates as stored in
1770// the header
1771//
1772// if not using the NCNR Transmission module, this function default to
1773// returning 0000, and no changes needed
1774Function getXYBoxFromFile(fname,x1,x2,y1,y2)
1775        String fname
1776        Variable &x1,&x2,&y1,&y2
1777       
1778        // return your bounding box coordinates or default values of 0
1779        x1=0
1780        y1=0
1781        x2=0
1782        y2=0
1783
1784        // your code returning value
1785        variable err
1786        string dfName = ""
1787        err = hdfRead(fname, dfName)
1788        //err not handled here
1789
1790       
1791        Wave wX1 = $(dfName+":reduce:x1")
1792        if (waveexists(wX1) == 0)
1793                //Waves don't exists which means an XY box has not been set for this file.
1794                //Hence return 0 bounding boxes (default)
1795        else
1796                x1 = wX1[0]
1797                Wave wX2 = $(dfName+":reduce:x2")
1798                x2 = wX2[0]
1799                Wave wY1 = $(dfName+":reduce:y1")
1800                y1 = wY1[0]
1801                Wave wY2 = $(dfName+":reduce:y2")
1802                y2 = wY2[0]
1803        endif
1804       
1805        KillWaves/Z wX1, wX2, wY1, wY2
1806        return(0)
1807
1808End
1809
1810// Writes the coordinates of the box to the header after user input
1811//
1812// box to sum over bounded (x1,y1) to (x2,y2)
1813//
1814// if not using the NCNR Transmission module, this function is null
1815Function WriteXYBoxToHeader(fname,x1,x2,y1,y2)
1816        String fname
1817        Variable x1,x2,y1,y2
1818
1819        String groupName = "/reduce"
1820        variable err
1821               
1822        Wave wX1
1823        Make/O/N=(1,1) wX1
1824        Wave wX2
1825        Make/O/N=(1,1) wX2
1826        Wave wY1
1827        Make/O/N=(1,1) wY1
1828        Wave wY2
1829        Make/O/N=(1,1) wY2
1830               
1831        wX1[0] = x1
1832        wX2[0] = x2
1833        wY1[0] = y1
1834        wY2[0] = y2     
1835       
1836        String varName = "x1"   
1837        err = hdfWrite(fname, groupName, varName,wX1)
1838        varName = "x2"
1839        err = hdfWrite(fname, groupName, varName,wX2)
1840        varName = "y1"
1841        err = hdfWrite(fname, groupName, varName,wY1)
1842        varName = "y2"
1843        err = hdfWrite(fname, groupName, varName,wY2)   
1844       
1845        KillWaves wX1,wX2,wY1,wY2
1846       
1847        //err not handled here
1848        return(0)       
1849
1850End
1851
1852// for transmission calculation, writes an NCNR-specific alphanumeric identifier
1853// (suffix of the data file)
1854//
1855//AJJ June 3 2010 - Note!! For ANSTO data the "suffix" is just the filename.
1856Function WriteAssocFileSuffixToHeader(fname,assoc_fname)
1857        String fname,assoc_fname
1858               
1859// your writer here
1860        Wave/T wAssoc_fname
1861        //nha ??? Should make this wave in our own DataFolder to avoid clashing names.
1862        Make/T /N=(1,1) wAssoc_fname
1863       
1864        String varName =""
1865        String groupName = "/reduce"
1866        if(isTransFile(fname))
1867                varName = "empty_beam_file_name"
1868        elseif(isScatFile(fname))
1869                varName = "transmission_file_name"
1870        endif
1871
1872        wAssoc_fname[0] = assoc_fname
1873        // your code returning value
1874        variable err
1875        err = hdfWrite(fname, groupName, varName, wAssoc_fname)
1876        KillWaves wAssoc_fname
1877        //err not handled here
1878        return(0)
1879end
1880
1881Function/S GetFileAssoc(fname)
1882        String fname
1883       
1884        String assoc_fname
1885        String groupName = ":reduce:"
1886       
1887        String varName = ""
1888        if(isTransFile(fname))
1889                varName = "empty_beam_file_name"
1890        elseif(isScatFile(fname))
1891                varName = "transmission_file_name"
1892        endif
1893       
1894        variable err
1895        string dfName = ""
1896        err = hdfRead(fname, dfName)
1897        //err not handled here
1898
1899        Wave/T wAssoc_fname = $(dfName+groupName+varName)
1900        if (waveexists(wAssoc_fname) == 1)
1901                assoc_fname =wAssoc_fname[0]   
1902        else
1903                assoc_fname = ""
1904        endif
1905        KillWaves/Z wAssoc_fname
1906       
1907        return(assoc_fname)
1908end
1909
1910Function hdfReadAttribute(fname, dfName, nodeName, nodeType, attributeStr)
1911// this is a copy of hdfRead, and could be incorporated back into hdfRead.
1912       
1913        String fname, &dfName, nodeName, attributeStr
1914        variable nodeType
1915        String nxentryName
1916        variable err=0,fileID   
1917        String cDF = getDataFolder(1), temp
1918        String fname_temp = ParseFilePath(0, fname, ":", 1, 0)
1919
1920       
1921        String fileSuffix
1922       
1923        if(strsearch(fname_temp,".nx.hdf",0,2)>=0)
1924                fileSuffix=".nx.hdf"
1925        else
1926                err = 1
1927                abort "unrecognised file suffix. Not .nx.hdf"
1928        endif
1929       
1930        dfName = "root:packages:quokka:"+removeending(fname_temp,fileSuffix)
1931       
1932        //folder must already exist i.e. hdfRead must have already been called
1933        if(!DataFolderExists(dfName))
1934                // possibly call an hdfRead from here
1935                return err
1936        endif
1937       
1938        //test for the name of nxentry
1939        if(DataFolderExists(dfName+":"+removeending(fname_temp,fileSuffix)))
1940                nxentryName = removeending(fname_temp,fileSuffix)
1941        elseif(DataFolderExists(dfName+":"+"entry1"))
1942                nxentryName = "entry1"
1943        else
1944                print "NXentry not found"
1945                return err
1946        endif
1947       
1948        //this is the stupid bit.
1949        // If you're looking for attributes of the root node, then nodename = "/"
1950        // If you're looking for attributes     of the nxentry node, then e.g. nodename ="/entry/instrument"
1951        // /entry is replaced with nxentryName
1952        nodeName = ReplaceString("entry", nodeName, nxentryName)       
1953       
1954        //convert nodeName to data folder string
1955        String dfNodeName = nodeName
1956        dfNodeName = ReplaceString( "/", nodeName, ":")
1957        dfName = dfName + dfNodeName
1958        if(nodeType == 2) //data item (dataset)
1959                //remove the end of dfName so that it points to a folder and not a dataset
1960                variable length = strlen(dfName)
1961                variable position = strsearch(dfName, ":", length, 1) // search backwards to find last :
1962                // to do - truncate string to remove dataset
1963                string truncate = "\"%0." + num2str(position) + "s\""
1964                sprintf dfName, truncate, dfName
1965        endif
1966       
1967        setDataFolder dfName
1968       
1969        try     
1970                HDF5OpenFile /R /Z fileID  as fname
1971                if(!fileID)
1972                        err = 1
1973                        abort "couldn't load HDF5 file"
1974                endif
1975
1976                HDF5LoadData /O /Q /Z /TYPE=(nodeType) /A=attributeStr, fileID, nodeName
1977
1978                if (V_flag!=0)
1979                        print "couldn't load attribute " + attributeStr
1980                endif
1981        catch
1982
1983        endtry
1984        if(fileID)
1985                HDF5CloseFile /Z fileID
1986        endif
1987
1988// add the name of the root node to dfName
1989// in the case of sensitivity files aka DIV files, don't append a root node to dfName
1990        if(DataFolderExists(dfName+":"+removeending(fname_temp,fileSuffix)))
1991                dfName = dfName+":"+removeending(fname_temp,fileSuffix)  //for files written Dec 2009 and after
1992        elseif(DataFolderExists(dfName+":"+"entry1"))
1993                dfName = dfName+":entry1" //for files written before Dec 2009
1994        endif
1995
1996        setDataFolder $cDF
1997        return err
1998end
1999
2000Function hdfRead(fname, dfName)
2001        //Reads hdf5 file into root:packages:quokka:fname
2002        //N. Hauser. 16/12/08
2003        // 29/1/09 - hdf file must have .nx.hdf or .div suffix
2004       
2005        String fname, &dfName
2006        variable err=0,fileID
2007        String cDF = getDataFolder(1), temp
2008        String fname_temp = ParseFilePath(0, fname, ":", 1, 0)
2009               
2010        String fileSuffix
2011        if(strsearch(fname_temp,".nx.hdf",0,2)>=0)
2012                fileSuffix=".nx.hdf"
2013        //elseif(strsearch(fname_temp,".div",0,2)>=0)
2014        //      fileSuffix=".div"
2015        else
2016                err = 1
2017                abort "unrecognised file suffix. Not .nx.hdf"
2018        endif
2019       
2020        dfName = "root:packages:quokka:"+removeending(fname_temp,fileSuffix)
2021       
2022        //if(!DataFolderExists(dfName))
2023        //      return err
2024        //else         
2025                newDataFolder /O root:packages
2026                newDataFolder /O /S root:packages:quokka
2027                newDataFolder /O /S $dfName
2028        //endif
2029       
2030        try     
2031                HDF5OpenFile /R /Z fileID  as fname
2032                if(!fileID)
2033                        err = 1
2034                        abort "couldn't load HDF5 file"
2035                endif
2036                HDF5LoadGroup /O /R /Z  :, fileID, "."
2037        catch
2038
2039        endtry
2040        if(fileID)
2041                HDF5CloseFile /Z fileID
2042        endif
2043
2044        // add the name of the root node to dfName
2045        // in the case of sensitivity files aka DIV files, don't append a root node to dfName
2046        if(DataFolderExists(dfName+":"+removeending(fname_temp,fileSuffix)))
2047                dfName = dfName+":"+removeending(fname_temp,fileSuffix)  //for files written Dec 2009 and after
2048        elseif(DataFolderExists(dfName+":"+"entry1"))
2049                dfName = dfName+":entry1" //for files written before Dec 2009
2050        endif
2051
2052        setDataFolder $cDF
2053        return err
2054end
2055
2056Function hdfWrite(fname, groupName, varName, wav)
2057        //Write Wave 'wav' to hdf5 file 'fname'
2058        //N. Hauser. nha 8/1/09
2059               
2060        String fname, groupName, varName
2061        Wave wav
2062       
2063        variable err=0, fileID,groupID
2064        String cDF = getDataFolder(1), temp
2065        String fname_temp = ParseFilePath(0, fname, ":", 1, 0)
2066        String NXentry_name
2067                       
2068        try     
2069                HDF5OpenFile/Z fileID  as fname  //open file read-write
2070                if(!fileID)
2071                        err = 1
2072                        abort "HDF5 file does not exist"
2073                endif
2074               
2075                //get the NXentry node name
2076                HDF5ListGroup /TYPE=1 fileID, "/"
2077                //remove trailing ; from S_HDF5ListGroup
2078                NXentry_name = S_HDF5ListGroup
2079                NXentry_name = ReplaceString(";",NXentry_name,"")
2080                if(strsearch(NXentry_name,":",0)!=-1) //more than one entry under the root node
2081                        err = 1
2082                        abort "More than one entry under the root node. Ambiguous"
2083                endif
2084                //concatenate NXentry node name and groupName   
2085                groupName = "/" + NXentry_name + groupName
2086                HDF5OpenGroup /Z fileID , groupName, groupID
2087
2088        //      !! At the moment, there is no entry for sample thickness in our data file
2089        //      therefore create new HDF5 group to enable write / patch command
2090        //      comment out the following group creation once thickness appears in revised file
2091       
2092                if(!groupID)
2093                        HDF5CreateGroup /Z fileID, groupName, groupID
2094                        //err = 1
2095                        //abort "HDF5 group does not exist"
2096                else
2097                        // get attributes and save them
2098                        //HDF5ListAttributes /Z fileID, groupName    this is returning null. expect it to return semicolon delimited list of attributes
2099                        //Wave attributes = S_HDF5ListAttributes
2100                endif
2101       
2102                HDF5SaveData /O /Z /IGOR=0  wav, groupID, varName
2103                if (V_flag != 0)
2104                        err = 1
2105                        abort "Cannot save wave to HDF5 dataset" + varName
2106                endif   
2107               
2108                //HDF5SaveData /O /Z /IGOR=0 /A=attributes groupID, varName
2109                //if (V_flag != 0)
2110                ////    err = 1
2111                //      abort "Cannot save attributes to HDF5 dataset"
2112                //endif
2113        catch
2114
2115        endtry
2116       
2117        if(groupID)
2118                HDF5CloseGroup /Z groupID
2119        endif
2120       
2121        if(fileID)
2122                HDF5CloseFile /Z fileID
2123        endif
2124
2125        setDataFolder $cDF
2126        return err
2127end
2128
2129Function DoEdgeCorrection(type)
2130        String type
2131        variable err = 0
2132       
2133        WAVE typeData=$("root:Packages:NIST:"+type+":data")
2134       
2135        //nha workaround for high counts on edges
2136        //copy second column to first column
2137        typeData[][0] = typeData[p][1]
2138        //copy second last column to last column
2139        typeData[][191] = typeData[p][190]
2140        //copy second row to first row
2141        typeData[0][] = typeData[1][q]
2142        //copy second last row to last row
2143        typeData[191][] = typeData[190][q]
2144       
2145        return err     
2146end
2147
2148////// OCT 2009, facility specific bits from ProDiv()
2149//"type" is the data folder that has the corrected, patched, and normalized DIV data array
2150//
2151// the header of this file is rather unimportant. Filling in a title at least would be helpful/
2152//
2153Function Write_DIV_File(type)
2154        String type
2155       
2156        // Your file writing function here. Don't try to duplicate the VAX binary format...
2157        WriteHeaderAndWork(type)
2158       
2159        return(0)
2160End
2161
2162////// OCT 2009, facility specific bits from MonteCarlo functions()
2163//"type" is the data folder that has the data array that is to be (re)written as a full
2164// data file, as if it was a raw data file
2165//
2166// not really necessary
2167//
2168Function/S Write_RawData_File(type,fullpath,dialog)
2169        String type,fullpath
2170        Variable dialog         //=1 will present dialog for name
2171       
2172        // Your file writing function here. Don't try to duplicate the VAX binary format...
2173        Print "Write_RawData_File stub"
2174       
2175        return(fullpath)
2176End
Note: See TracBrowser for help on using the repository browser.