function make_fast_histograms, scaler, id_first, ch_first, id_second, ch_second, $
                               ttsp, apthresh, fwin
;
; Create the fast neutron (Cat4) histograms from GRaND neutron event data.
;
; Notes:
;  The science scaler CAT4 counts give the number of double pulse events 
;  processed in a collection interval (validated by comparing the pattern
;  of scaler counts and total events across the BLP scintillators).  
;  Consequently, event tallies for each BLP scintillator are multiplied by 
;  the ratio of the scaler counts to the number of recorded events.
;
; Output(stored in result):
;
;  factor        -   factor to convert the cat1 histogram to total
;                    counts, given by the ratio between the number of
;                    CAT1 events counted by the scaler to the number of
;                    CAT1 events stored in the event buffer. 
;  fast_spectrum -   fast neutron spectrum
;  fast_var      -   fast neutron spectrum variance (assuming the scaler counts
;                    and channelized event counts are Poisson random variants)
;  total_count   -   number of events tallied
;  first         -   first interaction event spectrum (early-late)
;  second        -   second interaction event spectrum (early-late)
;  ttsp          -   time to second pulse event spectrum
;  first_early   -   first interaction event spectrum, early time window
;  second_early  -   second interaction event spectrum, early time window
;  first_late    -   first interaction event spectrum, late time window
;  second_late   -   second interaction event spectrum, late time window
;
; Detector indices: [0: pz, 1:my, 2:py, 3:mz]
;
; Example: 
; Plot the PZ fast spectrum with propagated error bars:
; plot, total(bfast.fast_spectrum[*,0],2), /ylog, min=1
; oploterr,total(bfast.fast_spectrum[*,0],2), sqrt(total(bfast.fast_var[*,0],2)), 3
;
; Plot the PZ time series sum of channels 2 to 40 of the fast spectrum
; plot, total(bfast.fast_spectrum[2:40,0],1)
; oploterr,total(bfast.fast_spectrum[2:40,0],1), sqrt(total(bfast.fast_var[2:40,0],1)), 3
;
; Version 1.0, 18-May-2009
; T. H. Prettyman
; Planetary Science Institute
;
; Instrument:
;   Dawn GRaND (GRD)
;


result={ factor           : fltarr(4)       ,  $  
         fast_spectrum    : fltarr(64,4)    ,  $
         fast_var         : fltarr(64,4)    ,  $
         total_count      : lonarr(4)       ,  $
         first            : lonarr(64,4)    ,  $ 
         second           : lonarr(64,4)    ,  $
         ttsp             : lonarr(256,4)   ,  $
         first_early      : lonarr(64,4)    ,  $
         second_early     : lonarr(64,4)    ,  $
         first_late       : lonarr(64,4)    ,  $
         second_late      : lonarr(64,4)       }

first_early=lonarr(64)
second_early=lonarr(64)
first_late=lonarr(64)
second_late=lonarr(64)   
ttsp1=lonarr(256)      

for index_blp=0,3 do begin
 index_count=where(id_first eq index_blp and ttsp ne 0, fcount)
 index=where(ttsp ge 4 and id_first eq index_blp and id_second eq index_blp and ch_second gt apthresh[index_blp], count)
 if (count ne 0) then begin
  ttsp1=histogram(ttsp[index],min=0,max=255,nbins=256)
  index_early=where(ttsp ge fwin[0] and ttsp le fwin[1]    $
               and id_first eq index_blp and id_second eq index_blp and ch_second gt apthresh[index_blp], count)
  if (count ne 0) then begin
   first_early=histogram(ch_first[index_early],min=0,max=63,nbins=64)
   second_early=histogram(ch_second[index_early],min=0,max=63,nbins=64)
  endif
  index_late=where(ttsp ge fwin[2] and ttsp le fwin[3]    $
                 and id_first eq index_blp and id_second eq index_blp and ch_second gt apthresh[index_blp], count)
  if (count ne 0) then begin
   first_late=histogram(ch_first[index_late],min=0,max=63,nbins=64)
   second_late=histogram(ch_second[index_late],min=0,max=63,nbins=64)
  endif
 endif
 result.total_count[index_blp]=fcount
 result.first[*,index_blp]=first_early-first_late
 result.first_early[*,index_blp]=first_early
 result.first_late[*,index_blp]=first_late
 result.second[*,index_blp]=second_early-second_late
 result.second_early[*,index_blp]=second_early
 result.second_late[*,index_blp]=second_late
 result.ttsp[*,index_blp]=ttsp1
 c_scaler=float(scaler[7+index_blp])
 c_event=float(result.total_count[index_blp])
 result.factor[index_blp]=c_scaler/c_event
 result.fast_spectrum[*,index_blp]=float(result.first[*,index_blp])*result.factor[index_blp]
 result.fast_var[*,index_blp]=c_scaler/c_event*((1./c_scaler+1./c_event)*(first_early-first_late)^2+   $
                                             first_early+first_late)
 
endfor

return, result

end