source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/ANSTO_DataReadWrite.ipf @ 890

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

Changed the 2D resolution calculation to include gravity. This seems to work, but will still require some testing. It requires no modifications to the smearing calculation, still using parallel and perpendicular directions.

Added a Fake Update feature to the RealTime? update. There are specific, separate instructions for how to set this up. Usef (possibly) for summer schools or demos.

Adjusted the 2D MonteCarlo? simulation to a corrected beam center. The Gauss Peak was not symmetric around the old beam center.

Corrected the AAO resolution smearing functions to work with cursors.

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        Make/O/D/N=1 wTransmission
1027        String groupName = "/reduce"
1028        String varName = "Transmission"
1029        wTransmission[0] = trans //
1030        // your code returning value
1031        variable err
1032        err = hdfWrite(fname, groupName, varName, wTransmission)
1033        KillWaves wTransmission
1034        //err not handled here
1035               
1036        return(0)
1037End
1038
1039//box counts from stored transmission calculation
1040// return 0 if not using NCNR transmission module
1041Function getBoxCounts(fname)
1042        String fname
1043        Variable value
1044       
1045        // do nothing if not using NCNR Transmission module
1046        variable err
1047        string dfName = ""
1048        err = hdfRead(fname, dfName)
1049        //err not handled here
1050
1051        Wave wBoxCounts = $(dfName+":reduce:boxCounts")
1052        if (waveexists(wBoxCounts) == 0)
1053                //boxcounts not yet set in  reduce group
1054                //return 0
1055                value = 0
1056        else
1057                value = wBoxCounts[0]
1058        endif
1059
1060        KillWaves/Z wBoxCounts
1061       
1062        return(value)
1063end
1064
1065//whole detector transmission
1066// return 0 if not using NCNR transmission module
1067Function getSampleTransWholeDetector(fname)
1068        String fname
1069        Variable value
1070       
1071        // your code returning value
1072        // ??? don't know what to put here. nha
1073        value=0
1074        return(value)
1075end
1076
1077//SampleThickness in centimeters
1078Function getSampleThickness(fname)
1079        String fname
1080        Variable value
1081
1082        // your code returning value
1083       
1084        variable err
1085        string dfName = ""
1086        err = hdfRead(fname, dfName)
1087        //err not handled here
1088
1089        if(WaveExists($(dfName+":sample:SampleThickness"))) //canonical location - a bit ugly and verbose, but that's just my opinion
1090                Wave wThickness = $(dfName+":sample:SampleThickness")
1091        elseif(WaveExists($(dfName+":sample:thickness")))
1092                Wave wThickness = $(dfName+":sample:thickness")
1093        else
1094                print "Can't find Sample Thickness in " + fname
1095        endif
1096                       
1097        value = wThickness[0]/10
1098        //value = 1 //??? temporary fix. nha
1099        KillWaves wThickness
1100       
1101        return(value)
1102end
1103
1104//sample thickness in cm
1105Function WriteThicknessToHeader(fname,thickness)
1106        String fname
1107        Variable thickness
1108       
1109        // your writer here
1110        Wave wThickness
1111        //nha ??? Should make this wave in our own DataFolder to avoid clashing names.
1112        Make /N=(1,1) wThickness
1113        String groupName = "/sample"
1114        String varName = "SampleThickness"
1115        wThickness[0] = thickness*10 //
1116        // your code returning value
1117        variable err
1118        err = hdfWrite(fname, groupName, varName, wThickness) //does not exist ???
1119        KillWaves wThickness
1120        //err not handled here
1121               
1122        return(0)
1123End
1124
1125//Sample Rotation Angle (degrees)
1126Function getSampleRotationAngle(fname)
1127        String fname
1128        Variable value
1129       
1130        // your code returning value
1131        variable err
1132        string dfName = ""
1133        err = hdfRead(fname, dfName)
1134        //err not handled here
1135
1136        Wave wSample_rotation_angle = $(dfName+":sample:sample_theta") //is this correct
1137        value = wSample_rotation_angle[0]
1138        KillWaves wSample_rotation_angle
1139               
1140        return(value)
1141end
1142
1143//temperature (C)
1144Function getTemperature(fname)
1145        String fname
1146       
1147        Variable value
1148       
1149        // your code returning value
1150       
1151        return(value)
1152end
1153
1154//field strength (Oe)
1155Function getFieldStrength(fname)
1156        String fname
1157       
1158        Variable value
1159       
1160        // your code returning value
1161       
1162        return(value)
1163end
1164
1165//center of beam xPos in pixel coordinates
1166Function getBeamXPos(fname)
1167        String fname
1168        Variable value
1169        // your code returning value
1170        variable err
1171        string dfName = ""
1172        err = hdfRead(fname, dfName)
1173        //err not handled here
1174
1175        if(WaveExists($(dfName+":instrument:reduce:BeamCenterX"))) //canonical location
1176                        Wave wBeamXPos = $(dfName+":instrument:reduce:BeamCenterX")
1177        elseif(WaveExists($(dfName+":instrument:parameters:BeamCenterX")))
1178                        Wave wBeamXPos = $(dfName+":instrument:parameters:BeamCenterX")
1179        else
1180                        print "Can't find BeamCenterX in" $fname
1181        endif
1182        value = wBeamXPos[0]   
1183        KillWaves wBeamXPos
1184       
1185        return(value)   
1186
1187end
1188
1189//beam center X pixel location (detector coordinates)
1190Function WriteBeamCenterXToHeader(fname,beamCenterX)
1191        String fname
1192        Variable beamCenterX
1193       
1194        // your writer here
1195        Wave wBeamCenterX
1196        //nha ??? Should make this wave in our own DataFolder to avoid clashing names.
1197        Make /N=(1,1) wBeamCenterX
1198        String groupName = "/instrument/reduce"
1199        String varName = "BeamCenterX"
1200        wBeamCenterX[0] = beamCenterX //
1201        // your code returning value
1202        variable err
1203        err = hdfWrite(fname, groupName, varName, wBeamCenterX)
1204        KillWaves wBeamCenterX
1205        //err not handled here
1206       
1207        return(0)
1208End
1209
1210//center of beam Y pos in pixel coordinates
1211Function getBeamYPos(fname)
1212        String fname
1213        Variable value
1214        // your code returning value
1215        variable err
1216        string dfName = ""
1217        err = hdfRead(fname, dfName)
1218        //err not handled here
1219
1220        if(WaveExists($(dfName+":instrument:reduce:BeamCenterZ"))) //canonical location
1221                        Wave wBeamYPos = $(dfName+":instrument:reduce:BeamCenterZ")
1222        elseif(WaveExists($(dfName+":instrument:parameters:BeamCenterZ")))
1223                        Wave wBeamYPos = $(dfName+":instrument:parameters:BeamCenterZ")
1224        else
1225                        print "Can't find BeamCenterZ in" $fname
1226        endif
1227        value = wBeamYPos[0]   
1228        KillWaves wBeamYPos
1229               
1230        return(value)
1231end
1232
1233//beam center Y pixel location (detector coordinates)
1234Function WriteBeamCenterYToHeader(fname,beamCenterY)
1235        String fname
1236        Variable beamCenterY
1237       
1238        // your writer here
1239        Wave wBeamCenterY
1240        //nha ??? Should make this wave in our own DataFolder to avoid clashing names.
1241        Make /N=(1,1) wBeamCenterY
1242        String groupName = "/instrument/reduce"
1243        String varName = "BeamCenterZ"
1244        wBeamCenterY[0] = beamCenterY //
1245        // your code returning value
1246        variable err
1247        err = hdfWrite(fname, groupName, varName, wBeamCenterY)
1248        KillWaves wBeamCenterY
1249        //err not handled here
1250
1251        return(0)
1252End
1253
1254
1255//sample to detector distance (meters)
1256Function getSDD(fname)
1257        String fname
1258        Variable value
1259        // your code returning value
1260        variable err
1261        string dfName = ""
1262        err = hdfRead(fname, dfName)
1263        //err not handled here
1264
1265        //workaround for bad HDF5 dataset
1266        if(WaveExists($(dfName+":instrument:parameters:L2"))) //canonical location
1267                Wave wSourceToDetectorDist = $(dfName+":instrument:parameters:L2")
1268        elseif(WaveExists($(dfName+":instrument:parameters:L2mm")))
1269                Wave wSourceToDetectorDist = $(dfName+":instrument:parameters:L2mm")
1270        elseif(WaveExists($(dfName+":instrument:parameters:derived_parameters:L2mm"))) 
1271                Wave wSourceToDetectorDist = $(dfName+":instrument:parameters:derived_parameters:L2mm")
1272        else
1273                print "Can't find L2 in " + fname
1274        endif
1275       
1276        value = wSourceToDetectorDist[0]/1000   
1277        KillWaves wSourceToDetectorDist
1278               
1279        return(value)
1280end
1281
1282//sample to detector distance (meters)
1283Function WriteSDDToHeader(fname,sdd)
1284        String fname
1285        Variable sdd
1286       
1287// your writer here
1288        Wave wSDD
1289        //nha ??? Should make this wave in our own DataFolder to avoid clashing names.
1290        Make /N=(1,1) wSDD
1291        String groupName = "/instrument/parameters"
1292        String varName = "L2"
1293        wSDD[0] = sdd * 1000 //
1294        // your code returning value
1295        variable err
1296        err = hdfWrite(fname, groupName, varName, wSDD)
1297        KillWaves wSDD
1298        //err not handled here
1299       
1300        return(0)
1301End
1302
1303//lateral detector offset (centimeters)
1304Function getDetectorOffset(fname)
1305        String fname
1306        Variable value
1307        // your code returning value
1308        variable err
1309        string dfName = ""
1310        err = hdfRead(fname, dfName)
1311        //err not handled here
1312
1313        Wave wDetectorOffset = $(dfName+":instrument:detector:detector_x") //is this correct
1314        value = wDetectorOffset[0]/10
1315        KillWaves wDetectorOffset
1316       
1317        return(value)
1318end
1319
1320//Beamstop diameter (millimeters)
1321Function getBSDiameter(fname)
1322        String fname
1323        Variable value
1324        // your code returning value
1325        variable err
1326        string dfName = ""
1327        err = hdfRead(fname, dfName)
1328        //err not handled here
1329
1330        if(WaveExists($(dfName+":instrument:parameters:BSdiam"))) //canonical location
1331                Wave wBSdiameter = $(dfName+":instrument:parameters:BSdiam")
1332        elseif(WaveExists($(dfName+":instrument:parameters:BSXmm")))
1333                Wave wBSdiameter = $(dfName+":instrument:parameters:BSXmm")
1334        else
1335                print "Can't find Beamstop Diameter in " + fname
1336        endif
1337        value = wBSdiameter[0]
1338        KillWaves wBSdiameter
1339       
1340        return(value)   
1341end
1342
1343//beam stop diameter (millimeters)
1344Function WriteBeamStopDiamToHeader(fname,bs)
1345        String fname
1346        Variable bs     
1347       
1348        // your writer here
1349        Wave wBS
1350        //nha ??? Should make this wave in our own DataFolder to avoid clashing names.
1351        Make /N=(1,1) wBS
1352        String groupName = "/instrument/parameters"
1353        String varName = "BSdiam"
1354        wBS[0] = bs //
1355        // your code returning value
1356        variable err
1357        err = hdfWrite(fname, groupName, varName, wBS)
1358        KillWaves wBS
1359        //err not handled here
1360        return(0)
1361End
1362
1363//source aperture diameter (millimeters)
1364Function getSourceApertureDiam(fname)
1365        String fname
1366        Variable value
1367        // your code returning value
1368        variable err
1369        string dfName = ""
1370        err = hdfRead(fname, dfName)
1371        //err not handled here
1372
1373        if(WaveExists($(dfName+":instrument:parameters:EApX")))
1374                Wave wSourceApertureDiam = $(dfName+":instrument:parameters:EApX") // canonical location
1375        elseif(WaveExists($(dfName+":instrument:parameters:derived_parameters:EApXmm")))
1376                Wave wSourceApertureDiam = $(dfName+":instrument:parameters:derived_parameters:EApXmm")
1377        else
1378                print "Can't find Source Aperture Diameter in " + fname
1379        endif   
1380        value = wSourceApertureDiam[0]
1381        KillWaves wSourceApertureDiam
1382       
1383        return(value)
1384end
1385
1386//Source Aperture diameter (millimeters)
1387Function WriteSourceApDiamToHeader(fname,source)
1388        String fname
1389        Variable source
1390       
1391        // your writer here
1392        Wave wsource
1393        //nha ??? Should make this wave in our own DataFolder to avoid clashing names.
1394        Make /N=(1,1) wsource
1395        String groupName = "/instrument/parameters"
1396        String varName = "EApX"
1397        wsource[0] = source //
1398        // your code returning value
1399        variable err
1400        err = hdfWrite(fname, groupName, varName, wsource)
1401        KillWaves wsource
1402        //err not handled here
1403        return(0)
1404End
1405
1406//sample aperture diameter (millimeters)
1407Function getSampleApertureDiam(fname)
1408        String fname
1409        Variable value
1410        // your code returning value
1411        variable err
1412        string dfName = ""
1413        err = hdfRead(fname, dfName)
1414        //err not handled here
1415
1416        if(WaveExists($(dfName+":sample:diameter"))) //canonical location
1417                Wave wSampleApertureDiam = $(dfName+":sample:diameter")
1418        elseif(WaveExists($(dfName+":instrument:parameters:autoSampleAp:diameter"))) //canonical location
1419                Wave wSampleApertureDiam = $(dfName+":instrument:parameters:autoSampleAp:diameter")
1420        elseif (WaveExists($(dfName+":instrument:sample_aperture:geometry:shape:SApXmm")))
1421                Wave wSampleApertureDiam = $(dfName+":instrument:sample_aperture:geometry:shape:SApXmm")
1422        else
1423                print "Can't find Sample Aperture Diameter in " + fname
1424        endif   
1425        value = wSampleApertureDiam[0]
1426        KillWaves wSampleApertureDiam
1427
1428        return(value)   
1429end
1430
1431//Sample Aperture diameter (millimeters)
1432Function WriteSampleApDiamToHeader(fname,source)
1433        String fname
1434        Variable source
1435       
1436        // your writer here
1437        Wave wsource
1438        //nha ??? Should make this wave in our own DataFolder to avoid clashing names.
1439        Make /N=(1,1) wsource
1440        String groupName = "/sample"
1441        String varName = "diameter"
1442        wsource[0] = source //
1443        // your code returning value
1444        variable err
1445        err = hdfWrite(fname, groupName, varName, wsource)
1446        KillWaves wsource
1447        //err not handled here
1448        return(0)
1449End
1450
1451//source AP to Sample AP distance (meters)
1452Function getSourceToSampleDist(fname)
1453        String fname
1454       
1455        Variable value
1456       
1457        // your code returning value
1458        variable err
1459        string dfName = ""
1460               
1461        err = hdfRead(fname, dfName)
1462        //err not handled here
1463       
1464        if(WaveExists($(dfName+":instrument:parameters:L1"))) //canonical location
1465                Wave wSourceToSampleDist = $(dfName+":instrument:parameters:L1")       
1466        elseif(WaveExists($(dfName+":instrument:parameters:L1mm")))
1467                Wave wSourceToSampleDist = $(dfName+":instrument:parameters:L1mm")
1468        elseif(WaveExists($(dfName+":instrument:parameters:derived_parameters:L1mm")))
1469                Wave wSourceToSampleDist = $(dfName+":instrument:parameters:derived_parameters:L1mm")
1470        else
1471                print "Can't find L1 in " + fname
1472        endif
1473       
1474        value = wSourceToSampleDist[0]/1000
1475        KillWaves wSourceToSampleDist   
1476               
1477        return(value)
1478end
1479
1480//Source aperture to sample aperture distance (meters)
1481Function WriteSrcToSamDistToHeader(fname,SSD)
1482        String fname
1483        Variable SSD
1484       
1485// your writer here
1486        Wave wSSD
1487        //nha ??? Should make this wave in our own DataFolder to avoid clashing names.
1488        Make /N=(1,1) wSSD
1489        String groupName = "/instrument/parameters"
1490        String varName = "L1"
1491        wSSD[0] = SSD * 1000 //input in metres, converted to mm for saving to file.
1492        // your code returning value
1493        variable err
1494        err = hdfWrite(fname, groupName, varName, wSSD)
1495        KillWaves wSSD
1496        //err not handled here
1497       
1498        return(0)
1499End
1500
1501//wavelength (Angstroms)
1502Function getWavelength(fname)
1503        String fname
1504        Variable value
1505        // your code returning value
1506        variable err
1507        string dfName = ""
1508        err = hdfRead(fname, dfName)
1509        //err not handled here
1510
1511        //      Wave wWavelength = $(dfName+":data:LambdaA")
1512        //change 22/12/09 nha
1513        // all these locations to be deprecated
1514        if(WaveExists($(dfName+":instrument:velocity_selector:Lambda")))  // canonical location
1515                Wave wWavelength = $(dfName+":instrument:velocity_selector:Lambda")
1516        elseif(WaveExists($(dfName+":data:Lambda")))
1517                Wave wWavelength = $(dfName+":data:Lambda")
1518        elseif(WaveExists($(dfName+":data:LambdaA")))
1519                Wave wWavelength = $(dfName+":data:LambdaA")
1520        elseif(WaveExists($(dfName+":instrument:velocity_selector:LambdaA")))
1521                Wave wWavelength = $(dfName+":instrument:velocity_selector:LambdaA")
1522        else
1523                print "Can't find Lambda in " + fname
1524        endif
1525        value = wWavelength[0] 
1526        KillWaves wWavelength
1527       
1528        return(value)
1529end
1530
1531//wavelength (Angstroms)
1532Function WriteWavelengthToHeader(fname,wavelength)
1533        String fname
1534        Variable wavelength
1535       
1536        // your writer here
1537        Wave wWavelength
1538        //nha ??? Should make this wave in our own DataFolder to avoid clashing names.
1539        Make /N=(1,1) wWavelength
1540        String groupName = "/instrument/velocity_selector"
1541        String varName = "Lambda"
1542        wWavelength[0] = wavelength //
1543        // your code returning value
1544        variable err
1545        err = hdfWrite(fname, groupName, varName, wWavelength)
1546       
1547        //and because Bill Hamilton is not happy with the NeXus naming, we write it to 3 other places
1548        //groupName = "/instrument/parameters"
1549        //err = hdfWrite(fname, groupName, varName, wWavelength)
1550        //velocity_selector group causes Igor crash in some cases
1551        //groupName = "/instrument/velocity_selector"
1552        //err = hdfWrite(fname, groupName, varName, wWavelength)
1553        //
1554        //groupName = "/data"
1555        //varName = "lambda"
1556        //err = hdfWrite(fname, groupName, varName, wWavelength)
1557       
1558        KillWaves wWavelength
1559        //err not handled here
1560
1561        return(0)
1562End
1563
1564
1565
1566//wavelength spread (FWHM)
1567Function getWavelengthSpread(fname)
1568        String fname
1569       
1570        Variable value
1571       
1572        // your code returning value
1573        variable err
1574        string dfName = ""
1575        err = hdfRead(fname, dfName)
1576        //err not handled here
1577       
1578        //velocity_selector group causes Igor crash
1579        if(WaveExists($(dfName+":instrument:velocity_selector:LambdaResFWHM_percent")))  //canonical location
1580                Wave wWavelengthSpread = $(dfName+":instrument:velocity_selector:LambdaResFWHM_percent")
1581        elseif(WaveExists($(dfName+":instrument:parameters:LambdaResFWHM_percent")))
1582                Wave wWavelengthSpread = $(dfName+":instrument:parameters:LambdaResFWHM_percent")
1583        else
1584                print "Can't find Wavelength Spread in " + fname
1585        endif
1586        value = wWavelengthSpread[0]   
1587        KillWaves wWavelengthSpread
1588       
1589        return(value)
1590end
1591
1592//wavelength spread (FWHM)
1593Function WriteWavelengthDistrToHeader(fname,wavelengthSpread)
1594        String fname
1595        Variable wavelengthSpread
1596       
1597        // your writer here
1598        Wave wWavelengthSpread
1599        //nha ??? Should make this wave in our own DataFolder to avoid clashing names.
1600        Make /N=(1,1) wWavelengthSpread
1601        //velocity_selector group causes Igor crash
1602        String groupName = "/instrument/velocity_selector"
1603        String varName = "LambdaResFWHM_percent"
1604
1605        wWavelengthSpread[0] = wavelengthSpread
1606        // your code returning value
1607        variable err
1608        err = hdfWrite(fname, groupName, varName, wWavelengthSpread)
1609        KillWaves wWavelengthSpread
1610        //err not handled here
1611               
1612        return(0)
1613End
1614
1615// physical x-dimension of a detector pixel, in mm
1616Function getDetectorPixelXSize(fname)
1617        String fname
1618        Variable value
1619       
1620        // your code here returning value
1621        variable err
1622        string dfName = ""
1623        err = hdfRead(fname, dfName)
1624        //err not handled here
1625
1626        Wave wActiveArea = $(dfName+":instrument:detector:active_height")
1627        Wave w_x_bin = $(dfName+":instrument:detector:x_bin")
1628        Variable numPixels = dimsize(w_x_bin, 0)
1629        value = wActiveArea[0]/numPixels
1630        KillWaves wActiveArea
1631        KillWaves w_x_bin
1632       
1633        return(value)
1634end
1635
1636// physical y-dimension of a detector pixel, in mm
1637Function getDetectorPixelYSize(fname)
1638        String fname
1639        Variable value
1640       
1641        // your code here returning value
1642        variable err
1643        string dfName = ""
1644        err = hdfRead(fname, dfName)
1645        //err not handled here
1646
1647        Wave wActiveArea = $(dfName+":instrument:detector:active_width")
1648        Wave w_y_bin = $(dfName+":instrument:detector:y_bin")
1649        Variable numPixels = dimsize(w_y_bin, 0)
1650        value = wActiveArea[0]/numPixels
1651        KillWaves wActiveArea
1652        KillWaves w_y_bin
1653                       
1654        return(value)
1655end
1656
1657//transmission detector count (unused, return 0)
1658//
1659Function getTransDetectorCounts(fname)
1660        String fname
1661       
1662        Variable value
1663       
1664        // your code returning value
1665       
1666        return(0)
1667end
1668
1669
1670//total count time (seconds)
1671Function getCountTime(fname)
1672        String fname
1673        Variable value
1674        // your code returning value
1675        variable err
1676        string dfName = ""
1677        err = hdfRead(fname, dfName)
1678        //err not handled here
1679
1680        Wave wTime1 = $(dfName+":monitor:bm1_time")
1681        value = wTime1[0]       
1682        KillWaves wTime1
1683       
1684        return(value)
1685end
1686
1687
1688Function getPhysicalWidth(fname)
1689        String fname
1690        Variable value
1691        // your code returning value
1692        variable err
1693        string dfName = ""
1694        err = hdfRead(fname, dfName)
1695        //err not handled here
1696
1697        Wave wPhysicalWidth = $(dfName+":instrument:detector:active_width")
1698        value = wPhysicalWidth[0]/10
1699        KillWaves wPhysicalWidth
1700               
1701        return(value)
1702end
1703
1704Function/S getSICSVersion(fname)
1705        String fname
1706        String value
1707        // your code returning value
1708        variable err
1709        string dfName = ""
1710        err = hdfRead(fname, dfName)
1711        //err not handled here
1712
1713        Wave/T wSICSVersion = $(dfName+":sics_release")
1714        value = wSICSVersion[0]
1715        KillWaves wSICSVersion
1716       
1717        return(value)
1718end
1719
1720Function/S getHDFversion(fname)
1721        String fname
1722        String value
1723        // your code returning value
1724        variable err
1725        string dfName = ""
1726        string attribute = "HDF5_Version"
1727        err = hdfReadAttribute(fname, dfName, "/", 1, attribute)
1728//      string attribute ="signal"
1729//      err = hdfReadAttribute(fname, dfName, "/entry/data/hmm_xy", 2, attribute)
1730        //err not handled here
1731
1732        Wave/T wHDF5_Version = $(dfName+":"+attribute)
1733        value = wHDF5_Version[0]       
1734//      KillWaves wHDF5_Version
1735
1736        return(value)
1737end
1738
1739
1740//reads the wavelength from a reduced data file (not very reliable)
1741// - does not work with NSORTed files
1742// - only used in FIT/RPA (which itself is almost NEVER used...)
1743//
1744// DOES NOT NEED TO BE CHANGED IF USING NCNR DATA WRITER
1745Function GetLambdaFromReducedData(tempName)
1746        String tempName
1747       
1748        String junkString
1749        Variable lambdaFromFile, fileVar
1750        lambdaFromFile = 6.0
1751        Open/R/P=catPathName fileVar as tempName
1752        FReadLine fileVar, junkString
1753        FReadLine fileVar, junkString
1754        FReadLine fileVar, junkString
1755        if(strsearch(LowerStr(junkString),"lambda",0) != -1)
1756                FReadLine/N=11 fileVar, junkString
1757                FReadLine/N=10 fileVar, junkString
1758                lambdaFromFile = str2num(junkString)
1759        endif
1760        Close fileVar
1761       
1762        return(lambdaFromFile)
1763End
1764
1765/////   TRANSMISSION RELATED FUNCTIONS    ////////
1766//box coordinate are returned by reference
1767//
1768// box to sum over is bounded (x1,y1) to (x2,y2)
1769//
1770// this function returns the bounding coordinates as stored in
1771// the header
1772//
1773// if not using the NCNR Transmission module, this function default to
1774// returning 0000, and no changes needed
1775Function getXYBoxFromFile(fname,x1,x2,y1,y2)
1776        String fname
1777        Variable &x1,&x2,&y1,&y2
1778       
1779        // return your bounding box coordinates or default values of 0
1780        x1=0
1781        y1=0
1782        x2=0
1783        y2=0
1784
1785        // your code returning value
1786        variable err
1787        string dfName = ""
1788        err = hdfRead(fname, dfName)
1789        //err not handled here
1790
1791       
1792        Wave wX1 = $(dfName+":reduce:x1")
1793        if (waveexists(wX1) == 0)
1794                //Waves don't exists which means an XY box has not been set for this file.
1795                //Hence return 0 bounding boxes (default)
1796        else
1797                x1 = wX1[0]
1798                Wave wX2 = $(dfName+":reduce:x2")
1799                x2 = wX2[0]
1800                Wave wY1 = $(dfName+":reduce:y1")
1801                y1 = wY1[0]
1802                Wave wY2 = $(dfName+":reduce:y2")
1803                y2 = wY2[0]
1804        endif
1805       
1806        KillWaves/Z wX1, wX2, wY1, wY2
1807        return(0)
1808
1809End
1810
1811// Writes the coordinates of the box to the header after user input
1812//
1813// box to sum over bounded (x1,y1) to (x2,y2)
1814//
1815// if not using the NCNR Transmission module, this function is null
1816Function WriteXYBoxToHeader(fname,x1,x2,y1,y2)
1817        String fname
1818        Variable x1,x2,y1,y2
1819
1820        String groupName = "/reduce"
1821        variable err
1822               
1823        Wave wX1
1824        Make/O/N=(1,1) wX1
1825        Wave wX2
1826        Make/O/N=(1,1) wX2
1827        Wave wY1
1828        Make/O/N=(1,1) wY1
1829        Wave wY2
1830        Make/O/N=(1,1) wY2
1831               
1832        wX1[0] = x1
1833        wX2[0] = x2
1834        wY1[0] = y1
1835        wY2[0] = y2     
1836       
1837        String varName = "x1"   
1838        err = hdfWrite(fname, groupName, varName,wX1)
1839        varName = "x2"
1840        err = hdfWrite(fname, groupName, varName,wX2)
1841        varName = "y1"
1842        err = hdfWrite(fname, groupName, varName,wY1)
1843        varName = "y2"
1844        err = hdfWrite(fname, groupName, varName,wY2)   
1845       
1846        KillWaves wX1,wX2,wY1,wY2
1847       
1848        //err not handled here
1849        return(0)       
1850
1851End
1852
1853// for transmission calculation, writes an NCNR-specific alphanumeric identifier
1854// (suffix of the data file)
1855//
1856//AJJ June 3 2010 - Note!! For ANSTO data the "suffix" is just the filename.
1857Function WriteAssocFileSuffixToHeader(fname,assoc_fname)
1858        String fname,assoc_fname
1859               
1860// your writer here
1861        Wave/T wAssoc_fname
1862        //nha ??? Should make this wave in our own DataFolder to avoid clashing names.
1863        Make/T /N=(1,1) wAssoc_fname
1864       
1865        String varName =""
1866        String groupName = "/reduce"
1867        if(isTransFile(fname))
1868                varName = "empty_beam_file_name"
1869        elseif(isScatFile(fname))
1870                varName = "transmission_file_name"
1871        endif
1872
1873        wAssoc_fname[0] = assoc_fname
1874        // your code returning value
1875        variable err
1876        err = hdfWrite(fname, groupName, varName, wAssoc_fname)
1877        KillWaves wAssoc_fname
1878        //err not handled here
1879        return(0)
1880end
1881
1882Function/S GetFileAssoc(fname)
1883        String fname
1884       
1885        String assoc_fname
1886        String groupName = ":reduce:"
1887       
1888        String varName = ""
1889        if(isTransFile(fname))
1890                varName = "empty_beam_file_name"
1891        elseif(isScatFile(fname))
1892                varName = "transmission_file_name"
1893        endif
1894       
1895        variable err
1896        string dfName = ""
1897        err = hdfRead(fname, dfName)
1898        //err not handled here
1899
1900        Wave/T wAssoc_fname = $(dfName+groupName+varName)
1901        if (waveexists(wAssoc_fname) == 1)
1902                assoc_fname =wAssoc_fname[0]   
1903        else
1904                assoc_fname = ""
1905        endif
1906        KillWaves/Z wAssoc_fname
1907       
1908        return(assoc_fname)
1909end
1910
1911Function hdfReadAttribute(fname, dfName, nodeName, nodeType, attributeStr)
1912// this is a copy of hdfRead, and could be incorporated back into hdfRead.
1913       
1914        String fname, &dfName, nodeName, attributeStr
1915        variable nodeType
1916        String nxentryName
1917        variable err=0,fileID   
1918        String cDF = getDataFolder(1), temp
1919        String fname_temp = ParseFilePath(0, fname, ":", 1, 0)
1920
1921       
1922        String fileSuffix
1923       
1924        if(strsearch(fname_temp,".nx.hdf",0,2)>=0)
1925                fileSuffix=".nx.hdf"
1926        else
1927                err = 1
1928                abort "unrecognised file suffix. Not .nx.hdf"
1929        endif
1930       
1931        dfName = "root:packages:quokka:"+removeending(fname_temp,fileSuffix)
1932       
1933        //folder must already exist i.e. hdfRead must have already been called
1934        if(!DataFolderExists(dfName))
1935                // possibly call an hdfRead from here
1936                return err
1937        endif
1938       
1939        //test for the name of nxentry
1940        if(DataFolderExists(dfName+":"+removeending(fname_temp,fileSuffix)))
1941                nxentryName = removeending(fname_temp,fileSuffix)
1942        elseif(DataFolderExists(dfName+":"+"entry1"))
1943                nxentryName = "entry1"
1944        else
1945                print "NXentry not found"
1946                return err
1947        endif
1948       
1949        //this is the stupid bit.
1950        // If you're looking for attributes of the root node, then nodename = "/"
1951        // If you're looking for attributes     of the nxentry node, then e.g. nodename ="/entry/instrument"
1952        // /entry is replaced with nxentryName
1953        nodeName = ReplaceString("entry", nodeName, nxentryName)       
1954       
1955        //convert nodeName to data folder string
1956        String dfNodeName = nodeName
1957        dfNodeName = ReplaceString( "/", nodeName, ":")
1958        dfName = dfName + dfNodeName
1959        if(nodeType == 2) //data item (dataset)
1960                //remove the end of dfName so that it points to a folder and not a dataset
1961                variable length = strlen(dfName)
1962                variable position = strsearch(dfName, ":", length, 1) // search backwards to find last :
1963                // to do - truncate string to remove dataset
1964                string truncate = "\"%0." + num2str(position) + "s\""
1965                sprintf dfName, truncate, dfName
1966        endif
1967       
1968        setDataFolder dfName
1969       
1970        try     
1971                HDF5OpenFile /R /Z fileID  as fname
1972                if(!fileID)
1973                        err = 1
1974                        abort "couldn't load HDF5 file"
1975                endif
1976
1977                HDF5LoadData /O /Q /Z /TYPE=(nodeType) /A=attributeStr, fileID, nodeName
1978
1979                if (V_flag!=0)
1980                        print "couldn't load attribute " + attributeStr
1981                endif
1982        catch
1983
1984        endtry
1985        if(fileID)
1986                HDF5CloseFile /Z fileID
1987        endif
1988
1989// add the name of the root node to dfName
1990// in the case of sensitivity files aka DIV files, don't append a root node to dfName
1991        if(DataFolderExists(dfName+":"+removeending(fname_temp,fileSuffix)))
1992                dfName = dfName+":"+removeending(fname_temp,fileSuffix)  //for files written Dec 2009 and after
1993        elseif(DataFolderExists(dfName+":"+"entry1"))
1994                dfName = dfName+":entry1" //for files written before Dec 2009
1995        endif
1996
1997        setDataFolder $cDF
1998        return err
1999end
2000
2001Function hdfRead(fname, dfName)
2002        //Reads hdf5 file into root:packages:quokka:fname
2003        //N. Hauser. 16/12/08
2004        // 29/1/09 - hdf file must have .nx.hdf or .div suffix
2005       
2006        String fname, &dfName
2007        variable err=0,fileID
2008        String cDF = getDataFolder(1), temp
2009        String fname_temp = ParseFilePath(0, fname, ":", 1, 0)
2010               
2011        String fileSuffix
2012        if(strsearch(fname_temp,".nx.hdf",0,2)>=0)
2013                fileSuffix=".nx.hdf"
2014        //elseif(strsearch(fname_temp,".div",0,2)>=0)
2015        //      fileSuffix=".div"
2016        else
2017                err = 1
2018                abort "unrecognised file suffix. Not .nx.hdf"
2019        endif
2020       
2021        dfName = "root:packages:quokka:"+removeending(fname_temp,fileSuffix)
2022       
2023        //if(!DataFolderExists(dfName))
2024        //      return err
2025        //else         
2026                newDataFolder /O root:packages
2027                newDataFolder /O /S root:packages:quokka
2028                newDataFolder /O /S $dfName
2029        //endif
2030       
2031        try     
2032                HDF5OpenFile /R /Z fileID  as fname
2033                if(!fileID)
2034                        err = 1
2035                        abort "couldn't load HDF5 file"
2036                endif
2037                HDF5LoadGroup /O /R /Z  :, fileID, "."
2038        catch
2039
2040        endtry
2041        if(fileID)
2042                HDF5CloseFile /Z fileID
2043        endif
2044
2045        // add the name of the root node to dfName
2046        // in the case of sensitivity files aka DIV files, don't append a root node to dfName
2047        if(DataFolderExists(dfName+":"+removeending(fname_temp,fileSuffix)))
2048                dfName = dfName+":"+removeending(fname_temp,fileSuffix)  //for files written Dec 2009 and after
2049        elseif(DataFolderExists(dfName+":"+"entry1"))
2050                dfName = dfName+":entry1" //for files written before Dec 2009
2051        endif
2052
2053        setDataFolder $cDF
2054        return err
2055end
2056
2057Function hdfWrite(fname, groupName, varName, wav)
2058        //Write Wave 'wav' to hdf5 file 'fname'
2059        //N. Hauser. nha 8/1/09
2060               
2061        String fname, groupName, varName
2062        Wave wav
2063       
2064        variable err=0, fileID,groupID
2065        String cDF = getDataFolder(1), temp
2066        String fname_temp = ParseFilePath(0, fname, ":", 1, 0)
2067        String NXentry_name
2068                       
2069        try     
2070                HDF5OpenFile/Z fileID  as fname  //open file read-write
2071                if(!fileID)
2072                        err = 1
2073                        abort "HDF5 file does not exist"
2074                endif
2075               
2076                //get the NXentry node name
2077                HDF5ListGroup /TYPE=1 fileID, "/"
2078                //remove trailing ; from S_HDF5ListGroup
2079                NXentry_name = S_HDF5ListGroup
2080                NXentry_name = ReplaceString(";",NXentry_name,"")
2081                if(strsearch(NXentry_name,":",0)!=-1) //more than one entry under the root node
2082                        err = 1
2083                        abort "More than one entry under the root node. Ambiguous"
2084                endif
2085                //concatenate NXentry node name and groupName   
2086                groupName = "/" + NXentry_name + groupName
2087                HDF5OpenGroup /Z fileID , groupName, groupID
2088
2089        //      !! At the moment, there is no entry for sample thickness in our data file
2090        //      therefore create new HDF5 group to enable write / patch command
2091        //      comment out the following group creation once thickness appears in revised file
2092       
2093                if(!groupID)
2094                        HDF5CreateGroup /Z fileID, groupName, groupID
2095                        //err = 1
2096                        //abort "HDF5 group does not exist"
2097                else
2098                        // get attributes and save them
2099                        //HDF5ListAttributes /Z fileID, groupName    this is returning null. expect it to return semicolon delimited list of attributes
2100                        //Wave attributes = S_HDF5ListAttributes
2101                endif
2102       
2103                HDF5SaveData /O /Z /IGOR=0  wav, groupID, varName
2104                if (V_flag != 0)
2105                        err = 1
2106                        abort "Cannot save wave to HDF5 dataset" + varName
2107                endif   
2108               
2109                //HDF5SaveData /O /Z /IGOR=0 /A=attributes groupID, varName
2110                //if (V_flag != 0)
2111                ////    err = 1
2112                //      abort "Cannot save attributes to HDF5 dataset"
2113                //endif
2114        catch
2115
2116        endtry
2117       
2118        if(groupID)
2119                HDF5CloseGroup /Z groupID
2120        endif
2121       
2122        if(fileID)
2123                HDF5CloseFile /Z fileID
2124        endif
2125
2126        setDataFolder $cDF
2127        return err
2128end
2129
2130Function DoEdgeCorrection(type)
2131        String type
2132        variable err = 0
2133       
2134        WAVE typeData=$("root:Packages:NIST:"+type+":data")
2135       
2136        //nha workaround for high counts on edges
2137        //copy second column to first column
2138        typeData[][0] = typeData[p][1]
2139        //copy second last column to last column
2140        typeData[][191] = typeData[p][190]
2141        //copy second row to first row
2142        typeData[0][] = typeData[1][q]
2143        //copy second last row to last row
2144        typeData[191][] = typeData[190][q]
2145       
2146        return err     
2147end
2148
2149////// OCT 2009, facility specific bits from ProDiv()
2150//"type" is the data folder that has the corrected, patched, and normalized DIV data array
2151//
2152// the header of this file is rather unimportant. Filling in a title at least would be helpful/
2153//
2154Function Write_DIV_File(type)
2155        String type
2156       
2157        // Your file writing function here. Don't try to duplicate the VAX binary format...
2158        WriteHeaderAndWork(type)
2159       
2160        return(0)
2161End
2162
2163////// OCT 2009, facility specific bits from MonteCarlo functions()
2164//"type" is the data folder that has the data array that is to be (re)written as a full
2165// data file, as if it was a raw data file
2166//
2167// not really necessary
2168//
2169Function/S Write_RawData_File(type,fullpath,dialog)
2170        String type,fullpath
2171        Variable dialog         //=1 will present dialog for name
2172       
2173        // Your file writing function here. Don't try to duplicate the VAX binary format...
2174        Print "Write_RawData_File stub"
2175       
2176        return(fullpath)
2177End
Note: See TracBrowser for help on using the repository browser.