function grd_read_l1a_science, pds_dir, after=after, fastwin=fastwin
;
; GRD_READ_L1A_SCI reads the GRaND L1A PDS science data within a 
; directory into an array of structures, representing a time series,
; for visualization and analysis.  This function was written to 
; illustrate how PDS L1A data are read using IDL and is not used
; for production of the L1B data.
;
; Data reduction:
;   Neutron event data are binned into histograms and stored in the
;   array of structures.
;
; Input:
;   pds_dir  - target directory containing the PDS data
;              (string) - refer to the example.
;   after    - (keyword) sets the after pulsing threshold for the CAT4 data analysis.  
;              If after is not set, then the default value is: [10,10,10,10].  If after 
;              has one element, the value is assigned to all four BLP scintillators.  
;              For example, after=7, results in [7,7,7,7].  Alternatively, different 
;              values can be set for each prism, for example, after=[10,15,15,10].
;   fastwin  - An array of values specifying the upper and lower time bins 
;              for subtracting accidentals from the prompt coincidences in the CAT4 
;              data analysis. The default is:
;              fastwin=transpose([[40L,90,200,250], [4,54,200,250], 
;                                 [4,54,200,250],   [40,90,200,250]])
;              The indexing scheme is fastwin[index_blp,index_parameter]
;
; Output:
;   result   - an array of structures providing a time series
;              of reduced data (see structure definition 
;              comments for sci_data below)
;
; Nomenclature:
;
;   .................
;   . ooooooooooooo .
;   . o           o .
;   . o           o .
;   . o    +Z     o .
;   . o   (PZ)    o .
;   . o           o .---> +Y (PY) 
;   . ooo       ooo .
;   .   o       o   .
;   .   o       o   .
;   .   ooooooooo   .
;   .               .
;   .................
;           |
;           v
;        +X (PX)
;
;  The coordinate system for GRaND is the same as that of the S/C.  
;  For the diagram above, the observer is looking in the -Z (MZ) direction 
;  and can see the outline of the phoswich assembly (o) on the +Z side of 
;  GRaND.  The phototubes are on the +X side and the scintillators 
;  are on the –X side.  During mapping at Vesta and Ceres, the asteroid’s 
;  surface is in the +Z direction.  The following provides a mapping 
;  between the direction nomenclature used during development and that 
;  used in flight:
;
;       Development       Flight     SPSI*    ID**
;          UP              PZ        1        0
;          DN              MZ        4        3
;          LF              MY        2        1
;          RT              PY        3        2
;  *Power supply index.  **ID is the index assigned to 
;  each of the plastic/phoswich scintillators in the event buffer
;  and the indices used in the arrays in the structure returned
;  by this function.
;
; Routines:
;   make_fast_histograms (function) 
;
; Example usage (PC):
;   path='F:\PDS\GRD_L0\ICO\07290_grd\'
;   directory=path+'GRD-L1A-071017-071018_090702\'
;   sci_data=GRD_READ_L1A_SCIENCE(directory, after=10, fastwin=[4,54,200,250])
;   plot, total(sci_data.bgo_hist,2), /ylog   ; plot of the BGO histogram
;   plot, total(sci_data.bgo_hist,1)          ; plot of the BGO total counts
;   plot, total(sci_data.fast[*,1],2)         ; plot of the CAT4 -Y fast neutron spectrum
;   plot, total(sci_data.second[*,1],2)       ; plot of the CAT4 -Y second interaction event spectrum
;   plot, total(sci_data.ttsp[*,1],2)         ; plot of the CAT4 -Y time-to-second-pulse event spectrum
;   plot, total(sci_data.bgo2[*,1],2)         ; plot of the CAT2 -Y BGO histogram
;   plot, total(sci_data.blp2[*,1],2)         ; plot of the CAT2 -Y BLP histogram
;   plot, total(sci_data.czt10[*,15],2)       ; plot of the CAT10 CZT histogram for CZT ID 15.
;
; Version 1.0, 17-May-2009:
;   Developed for the Dawn Data Workshop in Tennesee to
;   demonstrate the L1A data products from EMC and MGA. 
; Version 1.1, 24-Jun-2009:
;   Updated to construct (CAT4) fast neutron histograms and (CAT10) CZT
;   pulse height spectra.
;
; Written by
;   T. H. Prettyman
;   Planetary Science institute  
;
; Revision history:
;   *Routine does not execute on UNIX machines due to Windows path 
;    specification (DJL). Included OS-independent 
;    path separators; THP 25-Mar-2010 
;   *Modified routine to check for the existence of the target
;    directory; THP 30-Mar-2010
;   *Included packet sequence counter separately from scaler vector.
;    Scaler vector now has 23 entries; THP 2-Feb-2012
;   *Print out the error state message on open, THP 12-Jul-2013
;
; Instrument:
;   Dawn GRaND (GRD)
;

; Get the path separator for the target platform
psep=path_sep()

; Hardwired parameters
n_bits_czt=9      ; number of bits for the CZT histograms (max is 11)      
mge=3876L         ; maximum number of gamma ray events per science data record
mne=2800L         ; maximum number of neutron events per science data record

nshift_czt=11-n_bits_czt
n_channel_czt=2L^n_bits_czt

; 
; Structure definition
sci_data = { scet_utc    : 'yyyy-mm-ddThr:mn:sc'      ,  $ ; UTC spacecraft event time (SCET)
             sclk        : 0L                         ,  $ ; SCLK
             psc         : 0L                         ,  $ ; Packet sequence counter
             scaler      : ulonarr(23)                ,  $ ; Scaler data
             bgo_hist    : ulonarr(1024)              ,  $ ; BGO histogram
             czt10       : fltarr(n_channel_czt,16)   ,  $ ; CZT CAT10 single interaction histogram
             phos_pz     : ulonarr(256)               ,  $ ; +Z phoswich Cat1
             phos_mz     : ulonarr(256)               ,  $ ; -Z phoswich Cat1
             blp2        : lonarr(64,4)               ,  $ ; BLP CAT2 spectrum
             bgo2        : lonarr(64,4)               ,  $ ; BGO CAT2 spectrum
             fast        : lonarr(64,4)               ,  $ ; Fast - First interaction spectrum (scaler-corrected)
             second      : lonarr(64,4)               ,  $ ; Fast - Second interaction spectrum
             ttsp        : lonarr(256,4)                 } ; Fast - Time to second pulse

;
; Fast neutron parameters (fragile, not much error checking here)
apthresh=lonarr(4)
apthresh=[13,15,15,14]
if keyword_set(after) then begin
 if (n_elements(after) eq 1) then begin
  apthresh[*]=after ; threshold for after-pulsing
 endif else begin
  apthresh=after
 endelse
endif
apup=[64,63,63,63]  ; hardwired to max value
fwin=transpose([[40L,90,200,250], [4,54,200,250], [4,54,200,250], [40,90,200,250]])  ; fwin[index_blp,index_parameter]
if keyword_set(fastwin) then fwin=fastwin ; fast neutron time windows

; Default result
result=-1

; Find the data set prefix and do some checking
test=file_test(pds_dir)
if not test then begin
 print, ' ERROR: The directory does not exist. PDS_DIR = ', pds_dir
 return, result
endif
pos=strpos(pds_dir,'GRD-')
if (pos eq -1) then begin
 print, ' ERROR: The directory does not have the prefix GRD-. PDS_DIR = ', pds_dir
 return, result
endif
check=strlen(pds_dir)-1-pos
if (check ne 28) then begin
 print, ' ERROR: The PDS file name does not have the right length. '
 return, result
endif
prefix=strmid(pds_dir,pos,28)
search_pattern=prefix+'*'
files=file_search(pds_dir,search_pattern)
if (n_elements(files) ne 34) then begin   ; 34 files (TBL, DAT, LBL) are expected.
 print, ' ERROR: The archive is not complete. '
 return, result
endif
; From here on out there will be less error checking (hope for the best!)

; Find out how many records are in the time series...
file=pds_dir+'LEVEL1A_AUX'+psep+prefix+'-SCI-SCL.LBL'
val=get_pds_value(file,'FILE_RECORDS')
if (val eq -1) then begin
 print, ' ERROR: FILE_RECORDS not found by get_pds_value.'
 return, result
endif
reads, val, n_records, format='(i11)'
n_records=long(n_records)

; Create an array of structures
result=replicate(sci_data,n_records)

; Read the scaler data
FORMAT_STRING='(a19,i11,24(i6))'
scet_utc=' '
sclk=0L & psc=0L
scaler=ulonarr(23)
file=pds_dir+'LEVEL1A_AUX'+psep+prefix+'-SCI-SCL.TAB'
openr, lun, file, /get_lun, error=err
if (err ne 0) then begin
 print, !error_state.msg
 print, ' ERROR: Problem opening -SCI-SCL.TBL '
 return, result
endif
line=' '
for i=0,n_records-1 do begin
 readf, lun, scet_utc, sclk, psc, scaler, format=FORMAT_STRING
 result[i].scet_utc=scet_utc
 result[i].sclk=sclk
 result[i].psc=psc
 result[i].scaler=scaler
endfor
free_lun, lun


; Read the BGO histogram
FORMAT_STRING='(a19,i11,1024(i6))'
scet_utc=' '
sclk=0L
bgo_hist=ulonarr(1024)
file=pds_dir+'LEVEL1A_GAMMA'+psep+prefix+'-BGO.TAB'
openr, lun, file, /get_lun, error=err
if (err ne 0) then begin
 print, !error_state.msg
 print, ' ERROR: Problem opening -BGO.TBL '
 return, result
endif
line=' '
for i=0,n_records-1 do begin
 readf, lun, scet_utc, sclk, bgo_hist, format=FORMAT_STRING
 result[i].scet_utc=scet_utc
 result[i].sclk=sclk
 result[i].bgo_hist=bgo_hist
endfor
free_lun, lun

; Read the +Z phoswich histogram
FORMAT_STRING='(a19,i11,256(i6))'
scet_utc=' '
sclk=0L
phos_pz=ulonarr(256)
file=pds_dir+'LEVEL1A_NEUTRON'+psep+prefix+'-PHOS_PZ.TAB'
openr, lun, file, /get_lun, error=err
if (err ne 0) then begin
 print, !error_state.msg
 print, ' ERROR: Problem opening -PHOS_PZ.TBL '
 return, result
endif
line=' '
for i=0,n_records-1 do begin
 readf, lun, scet_utc, sclk, phos_pz, format=FORMAT_STRING
 result[i].scet_utc=scet_utc
 result[i].sclk=sclk
 result[i].phos_pz=phos_pz
endfor
free_lun, lun

; Read the -Z phoswich histogram
FORMAT_STRING='(a19,i11,256(i6))'
phos_mz=ulonarr(256)
file=pds_dir+'LEVEL1A_NEUTRON'+psep+prefix+'-PHOS_MZ.TAB'
openr, lun, file, /get_lun, error=err
if (err ne 0) then begin
 print, !error_state.msg
 print, ' ERROR: Problem opening -PHOS_MZ.TBL '
 return, result
endif
line=' '
for i=0,n_records-1 do begin
 readf, lun, scet_utc, sclk, phos_mz, format=FORMAT_STRING
 result[i].phos_mz=phos_mz
endfor
free_lun, lun


; Read the event mode neutron data
utc='yyyy-mm-ddThr:mn:sc '
sclk=0L
scaler=make_array(23,/long,value=0)
id_first=make_array(mne,/byte,value=0)
ch_first=make_array(mne,/byte,value=0)
id_second=make_array(mne,/byte,value=0)
ch_second=make_array(mne,/byte,value=0)
ttsp=make_array(mne,/byte,value=0)
file=pds_dir+'LEVEL1A_NEUTRON'+psep+prefix+'-EMN.DAT'
openr, lun, file, /get_lun, /swap_if_little_endian, err=err
if (err ne 0) then begin
 print, !error_state.msg
 print, ' ERROR: Problem opening -EMN.DAT '
 return, result
endif
for i=0, n_records-1 do begin
 readu, lun, utc, sclk, scaler, id_first, ch_first, id_second, ch_second, ttsp
 
 ; Build the CAT4 histograms for all of the BLP sensors
 fast=make_fast_histograms(scaler,id_first,ch_first,id_second,ch_second,ttsp,apthresh,apup,fwin)
 result[i].fast=fast.fast_spectrum
 result[i].second=fast.second
 result[i].ttsp=fast.ttsp
 
endfor
free_lun, lun


; Read the BLP2_PZ histogram
FORMAT_STRING='(a19,i11,64(i6))'
temp=ulonarr(64)
file=pds_dir+'LEVEL1A_NEUTRON'+psep+prefix+'-BLP2_PZ.TAB'
openr, lun, file, /get_lun, error=err
if (err ne 0) then begin
 print, !error_state.msg
 print, ' ERROR: Problem opening -BLP2_PZ.TBL '
 return, result
endif
line=' '
for i=0,n_records-1 do begin
 readf, lun, scet_utc, sclk, temp, format=FORMAT_STRING
 result[i].blp2[*,0]=temp
endfor
free_lun, lun

; Read the BLP2_MY histogram
FORMAT_STRING='(a19,i11,64(i6))'
temp=ulonarr(64)
file=pds_dir+'LEVEL1A_NEUTRON'+psep+prefix+'-BLP2_MY.TAB'
openr, lun, file, /get_lun, error=err
if (err ne 0) then begin
 print, !error_state.msg
 print, ' ERROR: Problem opening -BLP2_MY.TBL '
 return, result
endif
for i=0,n_records-1 do begin
 readf, lun, scet_utc, sclk, temp, format=FORMAT_STRING
 result[i].blp2[*,1]=temp
endfor
free_lun, lun

; Read the BLP2_PY histogram
FORMAT_STRING='(a19,i11,64(i6))'
temp=ulonarr(64)
file=pds_dir+'LEVEL1A_NEUTRON'+psep+prefix+'-BLP2_PY.TAB'
openr, lun, file, /get_lun, error=err
if (err ne 0) then begin
 print, !error_state.msg
 print, ' ERROR: Problem opening -BLP2_PY.TBL '
 return, result
endif
for i=0,n_records-1 do begin
 readf, lun, scet_utc, sclk, temp, format=FORMAT_STRING
 result[i].blp2[*,2]=temp
endfor
free_lun, lun

; Read the BLP2_MZ histogram
FORMAT_STRING='(a19,i11,64(i6))'
temp=ulonarr(64)
file=pds_dir+'LEVEL1A_NEUTRON'+psep+prefix+'-BLP2_MZ.TAB'
openr, lun, file, /get_lun, error=err
if (err ne 0) then begin
 print, !error_state.msg
 print, ' ERROR: Problem opening -BLP2_MZ.TBL '
 return, result
endif
for i=0,n_records-1 do begin
 readf, lun, scet_utc, sclk, temp, format=FORMAT_STRING
 result[i].blp2[*,3]=temp
endfor
free_lun, lun



; Read the BGO2_PZ histogram
FORMAT_STRING='(a19,i11,64(i6))'
temp=ulonarr(64)
file=pds_dir+'LEVEL1A_NEUTRON'+psep+prefix+'-BGO2_PZ.TAB'
openr, lun, file, /get_lun, error=err
if (err ne 0) then begin
 print, !error_state.msg
 print, ' ERROR: Problem opening -BGO2_PZ.TBL '
 return, result
endif
line=' '
for i=0,n_records-1 do begin
 readf, lun, scet_utc, sclk, temp, format=FORMAT_STRING
 result[i].bgo2[*,0]=temp
endfor
free_lun, lun

; Read the BGO2_MY histogram
FORMAT_STRING='(a19,i11,64(i6))'
temp=ulonarr(64)
file=pds_dir+'LEVEL1A_NEUTRON'+psep+prefix+'-BGO2_MY.TAB'
openr, lun, file, /get_lun, error=err
if (err ne 0) then begin
 print, !error_state.msg
 print, ' ERROR: Problem opening -BGO2_MY.TBL '
 return, result
endif
for i=0,n_records-1 do begin
 readf, lun, scet_utc, sclk, temp, format=FORMAT_STRING
 result[i].bgo2[*,1]=temp
endfor
free_lun, lun

; Read the BGO2_PY histogram
FORMAT_STRING='(a19,i11,64(i6))'
temp=ulonarr(64)
file=pds_dir+'LEVEL1A_NEUTRON'+psep+prefix+'-BGO2_PY.TAB'
openr, lun, file, /get_lun, error=err
if (err ne 0) then begin
 print, !error_state.msg
 print, ' ERROR: Problem opening -BGO2_PY.TBL '
 return, result
endif
for i=0,n_records-1 do begin
 readf, lun, scet_utc, sclk, temp, format=FORMAT_STRING
 result[i].bgo2[*,2]=temp
endfor
free_lun, lun

; Read the BGO2_MZ histogram
FORMAT_STRING='(a19,i11,64(i6))'
temp=ulonarr(64)
file=pds_dir+'LEVEL1A_NEUTRON'+psep+prefix+'-BGO2_MZ.TAB'
openr, lun, file, /get_lun, error=err
if (err ne 0) then begin
 print, !error_state.msg
 print, ' ERROR: Problem opening -BGO2_MZ.TBL '
 return, result
endif
for i=0,n_records-1 do begin
 readf, lun, scet_utc, sclk, temp, format=FORMAT_STRING
 result[i].bgo2[*,3]=temp
endfor
free_lun, lun


; Read the event mode gamma data and form CAT10 histograms 
utc='yyyy-mm-ddThr:mn:sc '
sclk=0L
scaler=make_array(23,/long,value=0)
id_czt=make_array(mge,/byte,value=0)
ch_czt=make_array(mge,/uint,value=0)
ch_bgo=make_array(mge,/uint,value=0)
czthist=fltarr(n_channel_czt,16)
file=pds_dir+'LEVEL1A_GAMMA'+psep+prefix+'-EMG.DAT'
openr, lun, file, /get_lun, /swap_if_little_endian, err=err
if (err ne 0) then begin
 print, !error_state.msg
 print, ' ERROR: Problem opening -EMG.DAT '
 return, result
endif
for i=0, n_records-1 do begin
 readu, lun, utc, sclk, scaler, id_czt, ch_czt, ch_bgo
 for j=0,16-1 do begin
  index=where(id_czt eq j and ch_czt ne 0 and ch_bgo eq 0, count)
  if (count gt 1) then begin
   ch=ishft(ch_czt[index],-nshift_czt)
   czthist[*,j]=histogram(ch,min=0,max=n_channel_czt-1,nbins=n_channel_czt)
  endif
 endfor
 result[i].czt10=czthist
endfor


return, result

end