; $Id: hfbs_readhfbs.pro,v 1.7 2005/02/09 16:58:01 dimeo Exp $
;###############################################################################
;
; NAME:
;  HFBS_READHFBS
;
; PURPOSE:
;  Reads in a raw HFBS data file.
;
; CATEGORY:
;  DAVE, HFBS, utility
;
; AUTHOR:
;   Robert M. Dimeo, Ph.D.
;   NIST Center for Neutron Research
;   100 Bureau Drive
;   Gaithersburg, MD 20899
;   Phone: (301) 975-8135
;   E-mail: robert.dimeo@nist.gov
;   http://www.ncnr.nist.gov/staff/dimeo
;
; LICENSE:
;  The software in this file is written by an employee of
;  National Institute of Standards and Technology
;  as part of the DAVE software project.
;
;  The DAVE software package is not subject to copyright protection
;  and is in the public domain. It should be considered as an
;  experimental neutron scattering data reduction, visualization, and
;  analysis system. As such, the authors assume no responsibility
;  whatsoever for its use, and make no guarantees, expressed or
;  implied, about its quality, reliability, or any other
;  characteristic. The use of certain trade names or commercial
;  products does not imply any endorsement of a particular product,
;  nor does it imply that the named product is necessarily the best
;  product for the stated purpose. We would appreciate acknowledgment
;  if the DAVE software is used of if the code in this file is
;  included in another product.
;
;###############################################################################
;
;
; Changed the normalization to the monitor.  Instead of just subtracting
; off a flat background and pulling the spectrum down to zero, we assume that
; there is a non-zero background in both the monitor and the detectors.  A
; background is subtracted off of both measurements until the level in the
; extremes of the normalized data agree with each other.
; This bit of code added (08/26/02) by R.M.Dimeo.  This now uses the set of
; minimization routines from Craig Markwardt called TNMIN.PRO.
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function norm_func,parms,_Extra = extra
; This function provides a metric for how similar the data is at the extremes
; of the dynamic range.
x = extra.x
y = extra.y
monitor = extra.monitor
bo = parms[0] & b1 = parms[1]
ynew = (y-bo)/(monitor-b1)
num = 40
nx = n_elements(x)
lowerMean = (moment(ynew[0:num-1]))[0]
upperMean = (moment(ynew[nx-1-num:nx-1]))[0]
f = abs(lowerMean-upperMean)
return,f
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function norm_func_one,parms,_Extra = extra
; This function provides a metric for how similar the data is at the extremes
; of the dynamic range.
x = extra.x
y = extra.y
monitor = extra.monitor
bmon = extra.bmon
bo = parms[0]
ynew = (y-abs(bo))/(monitor-abs(bmon))
num = 30
nx = n_elements(x)
lowerMean = (moment(ynew[0:num-1]))[0]
upperMean = (moment(ynew[nx-1-num:nx-1]))[0]
f = abs(lowerMean-upperMean)
return,f
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro normalizeData,x,y,yerr,monitor,monerr,$
                  pguess = pguess, $
                  presult = presult, $
                  ynormout=ynormout,$
                  ynormerrout=ynormerrout,$
                  maxiter = maxiter, $
                  errmsg = errmsg

;if n_elements(pguess) eq 0 then pguess = [0.01*max(y),0.0]
;if n_elements(pscale) eq 0 then pscale = [0.1*max(y),0.01*max(monitor)]
parinfo = replicate({limited:[0,0], $
                     limits:[0.D,0]},2)

parinfo[0].limited[0] = 1
parinfo[0].limits[0]  = 0.D
parinfo[1].limited[0] = 1
parinfo[1].limits[0]  = 0.D
functargs = {x:x,y:y,monitor:monitor}

if n_elements(maxiter) eq 0 then maxiter = 500
R = TNMIN('norm_func', autoderivative = 1,pguess,$
    functargs = functargs,quiet = 1,parinfo = parinfo,$
    status = status,maxiter = maxiter)
presult = R
;R = DAVE_AMOEBA(1.0e-5, SCALE=pscale, P0 = pguess, $
;    function_name = 'norm_func',FUNCTION_VALUE=fval,$
;    x = x,y = y,monitor = monitor,nmax = 10000)

IF N_ELEMENTS(R) EQ 1 THEN begin
  errmsg = 1
  ynormerrout = 0.0
  ynormout = 0.0
  return
endif
errmsg = 0
ynormout = (y-r[0])/(monitor-r[1])
; Now propagate the error (Gaussian statistics) resulting
; from the normalization.
u = y
v = monitor
uerr = yerr
verr = sqrt(v);monerr
ynormerrout = (u/v)*sqrt((uerr/u)^2+(verr/v)^2)
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro normalizeData_one,x,y,yerr,monitor,monerr,$
                      bmon = bmon, $
                      pguess = pguess, $
                      pscale = pscale, $
                      presult = presult, $
                      ynormout=ynormout,$
                      ynormerrout=ynormerrout,$
                      errmsg = errmsg

if n_elements(pguess) eq 0 then pguess = 0.1*max(y)
if n_elements(pscale) eq 0 then pscale = 0.025*max(y)
functargs = {x:x,y:y,monitor:monitor}
R = TNMIN('cm_norm_func', autoderivative = 1,[0D, 0D],$
    functargs = functargs,quiet = 1)
;R = DAVE_AMOEBA(1.0e-6, SCALE=pscale, P0 = pguess, $
;    function_name = 'norm_func_one',FUNCTION_VALUE=fval,$
;    x = x,y = y,monitor = monitor,nmax = 10000,bmon = bmon)
presult = abs(R)
IF N_ELEMENTS(R) EQ 1 and R[0] eq -1 THEN begin
  errmsg = 1
  ynormerrout = 0.0
  ynormout = 0.0
  return
endif
errmsg = 0
ynormout = (y-abs(r[0]))/(monitor-bmon)
; Now propagate the error (Gaussian statistics) resulting
; from the normalization.
u = y
v = monitor
uerr = yerr
verr = sqrt(v);monerr
ynormerrout = (u/v)*sqrt((uerr/u)^2+(verr/v)^2)
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;==========================================
PRO hfbs_readhfbs,file,camType,energy,data,error,ndet,nchan,comments,integ_int,$
                  SUMWBM = thisSumWBM,$
                  SUMTBM = thisSumTBM,$
                  SUMFC = thisSumFC,$
                  LIVETIME = thisLiveTime,$
                  Norm2One = thisNorm2One,$
                  NoConvert2Energy = thisNoConvert2Energy,$
                  DetectorRates = thisDetectorRates,$
                  Doppler = thisDopplerRate,$
                  NORM2MON = thisNorm2Mon,$
                  ncycles = countDop,$
                  temperature = thisTemperature,$
                  DETRates = detrates
; reads in an HFBS file
; These next couple of lines were added so that the user does not have to determine which cam
; to use when analyzing data.  The triangle cam was put in the instrument in June so that all
; previous data automatically uses the sine cam and all subsequent data uses the triangle cam.
;
; Modified 9/28/2000 by RMD
; Modified 10/11/2000 by RMD to include the keyword parameters
;
; If NoNorm2One is present then program will not normalize data to one
; if NoConvert2Energy is present then program will not convert to energy
;
; Modified 12/12/01 (RMD):	Error propagation corrected for monitor normalization
; Modified 01/23/02 (RMD):  Replaced some "FOR" loops with array operations, thus increasing the speed.
;
if n_elements(thisNorm2One) eq 0 then thisNorm2One = 0
detRates = fltarr(16)

res = strpos(file,'.hfbs')
compnum = double(strmid(file,res-11,8))
if compnum ge 20000601d or compnum le 19990824d then begin
  camType = 1	; triangle cam
endif else begin
  camType = 0	; sine cam
endelse

OPENR, lun, file, /GET_LUN
nlines = 0
nheader = 9
dumstr = ''
header=strarr(nheader)
for i=0,nheader-1 do begin
  readf,lun,dumstr
  header[i] = dumstr
  nlines = nlines+1
endfor
comments = header

thisLiveTime = (strmid(comments[7],21))
timebin = float(strmid(header[8],21,2))
;timebin = 50.0

; Read in the Log Doppler Frequency vector
READF, lun, dumstr & nlines = nlines+1
dopFreq = 0.0
READF, lun, dumstr & nlines = nlines+1
sumDop = float(dumstr)
countDop = 1
while strmid(dumstr,0,1) ne '#' do begin
  readf, lun, dumstr & nlines = nlines+1
  countDop = countDop + 1
  if (strmid(dumstr,0,1) ne '#') then  begin
  dopFreq = float(dumstr)
  sumDop = temporary(sumDop) + dopFreq
  endif
endwhile
countDop = temporary(countDop) - 1
DopFreq = sumDop/countDop
thisDopplerRate = DopFreq
realchannels = 0.0

; Read in the Log white beam vector
sumWBM = 0.0
WBM = 0.0
for i = 0,countDop - 1 do begin
  READF,lun, dumstr & nlines = nlines+1
  WBM = WBM+float(dumstr)
  sumWBM = WBM
endfor
thisSumWBM = sumWBM

changetoenergy = 1
; Read in the Log trans beam vector
sumTBM = 0.0
TBM = 0.0
READF,lun, dumstr  & nlines = nlines+1  ; read the text header for the TBM
for i = 0,countDop - 1 do begin
  READF,lun, dumstr & nlines = nlines+1
  TBM = TBM+float(dumstr)
  sumTBM = TBM
endfor
thisSumTBM = sumTBM
readf,lun,dumstr & nlines = nlines + 1
free_lun,lun

if strmid(dumstr,0,3) eq '# D' then begin
  openr,lun,file
  for i = 0,nlines-2 do begin
    readf,lun,dumstr
  endfor
  skipthis = 1
endif else begin
  openr,lun,file
  for i = 0,nlines-2 do begin
    readf,lun,dumstr
  endfor
  skipthis = 0
endelse

; Read in the log sample temperature
if skipthis eq 0 then begin
  tempPtr = ptr_new(/allocate_heap)
  READF,lun, dumstr     ; '#Log sample temp'
  READF,lun, dumstr
  if float(dumstr) ne 0.0 then begin
    temperature = fltarr(countdop)
    temperature[0] = float(dumstr)
    for i = 1,countdop-1 do begin
      readf,lun,dumstr
      temperature[i] = float(dumstr)
    endfor
    *tempPtr = temperature
  endif
endif
timeMinutes = (double(thisLiveTime))/60d
if (n_elements(*tempPtr) ne 0) then begin
  thisTemperature = temperature
endif else begin
  thisTemperature = 0.0
endelse
ptr_free,tempPtr

CASE strmid(header[2],21,2) OF
've':  BEGIN
      timeMinutes = (double(thisLiveTime))/60d
      ndet = 20
      nchan = 4096
      energy = dblarr(nchan)
      de = 100.0/(nchan-1.0)
      energy = -50.0+de*dindgen(nchan)
      energy = 4.0*temporary(energy)
      ; These next five lines replace a nested for loop.  These lines
      ; use array operations rather than looping to increase speed.
      ; A 3% decrease in time to create data was found with a particular
      ; data set.
      datStr = {instring:'',data:fltarr(2,nchan)}
      str_arr = replicate(datStr,ndet)
      readf,lun,str_arr
      data = transpose(str_arr.data[1,*])
      error = sqrt(data)
      free_lun, lun

      wherezero = where(error eq 0,thisnum)
      if thisnum gt 0 then begin
        error[wherezero] = 1.0
      endif
      whereok = where(energy le 50. and energy ge -50.,nchan)
      newenergy = dblarr(nchan)
      newenergy[0:nchan-1] = energy[whereok]
      newdata = dblarr(ndet,nchan)
      newerror = dblarr(ndet,nchan)
      newdata[*,0:nchan-1] = data[*,whereok]
      newerror[*,0:nchan-1] = error[*,whereok]
      energy = dblarr(nchan)
      energy = newenergy
      data = dblarr(ndet,nchan)
      error = dblarr(ndet,nchan)
      data = newdata
      error = newerror
	  ; calculate the count rates in each detector (cps)
      detRates = (total(data[1:16,*],2))/float(thisLivetime)

        lamo = 6.27
        eo = 81.81/lamo^2
        vo = 3956.0/lamo

        totdetcounts = dblarr(ndet)
        totdetcounts = total(double(data[0:ndet-1,0:nchan-1]),1)
        thisDetectorRates = dblarr(ndet)
        thisDetectorRates = totdetcounts/timeMinutes
        moncount = 0l
        moncount = total(data[0,0:nchan-1])

; normalize the detectors to p(E) and multiply by integral of p(E)
       finitecounts = where(data[19,*] ne 0.0)
; replace the loop with an appropriate array operation

       uvec = 1+bytarr(16)
       pe = data[19,finitecounts]/(total(data[19,finitecounts]*deriv(energy[finitecounts])))
       data[1:16,finitecounts] = data[1:16,finitecounts]/(uvec#pe)
       error[1:16,finitecounts] = error[1:16,finitecounts]/(uvec#pe)

; normalize the monitor to its P(v)
       finitecounts = where(data[18,*] ne 0.0)
       data[0,finitecounts] = (data[0,finitecounts]/data[18,finitecounts])*$
                              total(data[18,finitecounts]*deriv(energy[finitecounts]))
       error[0,finitecounts] = (error[0,finitecounts]/data[18,finitecounts])*$
                              total(data[18,finitecounts]*deriv(energy[finitecounts]))
        PvMon = FLTARR(1,nchan)
        tempmon = FLTARR(nchan)
        tempmon = REPLICATE(1,nchan)
        PvMon[0,0:nchan-1] = tempmon[0:nchan-1]
      END
'ps':  BEGIN
      ndet = 20
      nchan = 4096
      data = FLTARR(ndet,nchan)
      error = FLTARR(ndet,nchan)

      FOR i = 0,ndet-1 DO BEGIN
         READF, lun, dumstr1
         FOR j = 0,nchan-1 DO BEGIN
           READF, lun, ch, intensity
           data[i,j] = intensity
           error[i,j] = SQRT(intensity)
         ENDFOR
      ENDFOR

        FREE_LUN, lun
        totdetcounts = INTARR(ndet)

        FOR i=0,ndet-1 DO BEGIN
          totdetcounts[i] = TOTAL(data[i,0:nchan-1])
        ENDFOR

        energy = FLTARR(nchan)
        lamo = 6.27
        eo = 81.81/lamo^2
        vo = 3956.0/lamo
        chan2ueV = (2.0*eo*1.0e3/vo)/(2048.0/20.0)
        energy = 1.0+FINDGEN(nchan)


        moncount = 0l
        moncount = TOTAL(data[0,0:nchan-1])
      END

'ti': BEGIN
      ndet = 20
      nchan = 4096
      data = FLTARR(ndet,nchan)
      error = FLTARR(ndet,nchan)
      FOR i = 0,ndet-1 DO BEGIN
         READF, lun, dumstr
         FOR j = 0,nchan-1 DO BEGIN
           READF, lun, ch, intensity
           data[i,j] = intensity
           error[i,j] = SQRT(intensity)
         ENDFOR
      ENDFOR

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; close the HFBS data file
      free_lun,lun
data[0:ndet-1,nchan-1] = 0.0*findgen(ndet)
; Do we want to convert to energy or leave as time-of-arrival?
if n_elements(thisNoConvert2Energy) eq 0 then begin

      channels = fltarr(nchan)
      channels = findgen(nchan)
      guesschannels = fix(1.0/(1.0e-6*DopFreq*timebin))
      realchannels = guesschannels
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; begin the time to energy conversion
        lamo = 6.27
        eo = 81.81/lamo^2
        vo = 3956.0/lamo
        dms = 2.25
        dsad = 4.15
        dsafc = 4.24    ; distance to fission chamber
        nchannels = fix(dsafc/(vo*TimeBin*1.0e-6))
        if realchannels lt nchan then begin
          data[0,0:realchannels-1]=shift(data[0,0:realchannels-1],nchannels)
          error[0,0:realchannels-1]=shift(error[0,0:realchannels-1],nchannels)
        endif
        camData = fltarr(2,nchan)
        dataNew = fltarr(16,nchan)
        energy = fltarr(nchan)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Perform a circular shift of the data up to realchannels
; in order to center the data in the effective time window
        ytest = fltarr(nchan)
        for i=5,13 do begin
          ytest[0:nchan-1]=ytest[0:nchan-1]+data[i,0:nchan-1]
        endfor
        numtest = 1300
        y1max = max(ytest[0:numtest-1]) & y2max = max(ytest[numtest:4095])
        y1ind = where(ytest[0:numtest-1] eq y1max)
        y2ind = where(ytest[numtest:4095] eq y2max)
        y2ind = y2ind + numtest

        y1ind = y1ind[0] & y2ind = y2ind[0]

        nhalf = fix((y2ind-y1ind)/2.0)
        realhalf = fix(0.5*realchannels)
        halfshift = (nhalf-realhalf)

if DopFreq gt 8.0 and compnum ge 19990824d then begin
        for i = 0,16 do begin
          data[i,0:realchannels-1]=shift(data[i,0:realchannels-1],halfshift)
          error[i,0:realchannels-1]=shift(error[i,0:realchannels-1],halfshift)
        endfor
endif
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
        ytest = fltarr(nchan)
        for i=5,13 do begin
          ytest[0:nchan-1]=ytest[0:nchan-1]+data[i,0:nchan-1]
        endfor
; Added the following line for numtest on 11/1/2000
;       numtest = 1300
        numtest = 1000
        y1max = max(ytest[0:numtest-1]) & y2max = max(ytest[numtest:4095])
        y1ind = where(ytest[0:numtest-1] eq y1max)
        y2ind = where(ytest[numtest:4095] eq y2max)
        y2ind = y2ind + numtest

        arbchannels = 4096
        nfold = fix((y2ind-y1ind)/2)+y1ind+arbchannels/2

        nfold = nfold[0]
        nchan = 4096

        channels = fltarr(nchan+arbchannels)
        channels = findgen(nchan+arbchannels)

        energy = fltarr(nchan)
        energy = indgen(nchan)
        newdata = fltarr(ndet,nchan+arbchannels)
        jfinal = abs(nfold-fix(nchan))
        longdata = fltarr(ndet,nchan+arbchannels)
        longdata[0:ndet-1,arbchannels/2:nchan+arbchannels/2-1]=data[0:ndet-1,0:nchan-1]

        for j=0,nfold do begin
          if nfold+j lt nchan+arbchannels then begin
            newdata[0:ndet-1,j+arbchannels/2]=longdata[0:ndet-1,nfold+j]+$
                                              longdata[0:ndet-1,nfold-j]
          endif else begin
            newdata[0:ndet-1,j+arbchannels/2]=0.0
          endelse
        endfor


     nfold2 = arbchannels/2+fix(realchannels/2)
     numchannels = nfold2-nfold
     ynew = fltarr(ndet,nchan+arbchannels)
     ynew = newdata
     y = fltarr(ndet,nchan+arbchannels)
     y = newdata
     for j=1,numchannels-1 do begin
       ynew[0:ndet-1,nfold+j]=y[0:ndet-1,nfold+j]+y[0:ndet-1,nfold2+numchannels-j]
     endfor

     ynew[0:ndet-1,nfold2:nchan+arbchannels-1]=0.0*findgen(ndet,nchan+arbchannels-nfold2)

     indlo = nfold-y1ind-(y2ind-y1ind)/2
     indhi = nfold2
     indlo = indlo[0]
     indhi = indhi[0]
     numind = indhi-indlo+1

     ydatnew = fltarr(ndet,numind)
     ydatnew[0:ndet-1,0:numind-1]=ynew[0:ndet-1,indlo:indhi]

       realchannels = numind
       artchannels = round(1.0/(2.0e-6*TimeBin*DopFreq))
       ndiff = realchannels-artchannels+1

       realchannels = realchannels+ndiff
       camdata = fltarr(2,2*(realchannels))
       camdata = sine_cam_funQ(2*(realchannels),DopFreq)

      ;;;Determine energy from cam profile
      ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
       y1new = fltarr(ndet,realchannels)
       y1new[0:ndet-1, 0:realchannels-1-ndiff] = ydatnew[0:ndet-1, 0:realchannels-1-ndiff]

       velocity = fltarr(realchannels)

       if CamType eq 0 then begin               ; sine cam
         scalefactor = 1.0
         camdata = sine_cam_funQ(2*(realchannels),DopFreq)
         velocity[0:realchannels-1] = scalefactor*camData[1,0:realchannels-1]
       endif else begin                         ; triangle cam
         camData = triangle_cam_funQ(2*realchannels,DopFreq)
         velocity[0:realchannels-1] = camData[1,0:realchannels-1]
       endelse
       energy = fltarr(realchannels)

       k1 = 5.22707e-3 & k2 = 1261.64

       energy[0:realchannels-1] = k1*(k2*velocity[0:realchannels-1]+$
                                  (velocity[0:realchannels-1])^2)

       if realchannels mod 2 eq 0 then begin
         finalchannels = realchannels/2
         newchannels = realchannels-finalchannels+1
       endif else begin
         finalchannels = (realchannels+1)/2
         newchannels = realchannels-finalchannels
       endelse

       newchannels = realchannels-ndiff
       newenergy = fltarr(newchannels)
       newenergy[0:newchannels-1] = energy[0:newchannels-1]
       newdata = fltarr(ndet,newchannels)
       newdata[0:ndet-1,0:newchannels-1] = y1new[0:ndet-1,0:newchannels-1]

        nchan = newchannels
        energy = fltarr(nchan)
        energy[0:nchan-1] = newenergy[0:nchan-1]
        data = fltarr(ndet,nchan)
        error = fltarr(ndet,nchan)


        ;  shift the low angle detectors the appropriate amount so that
        ;  their elastic peaks line up with the high angle detectors

        for i = 1,16 do begin
          ymax = max(newdata[i,0:nchan-1]) & ymax = ymax[0]
          yind = where(newdata[i,0:nchan-1] eq ymax) & yind = yind[0]
          emin = min(abs(energy[0:nchan-1])) & emin = emin[0]
          eind = where(abs(energy[0:nchan-1]) eq abs(emin)) & eind = eind[0]
          newdata[i,0:nchan-1] = shift(newdata[i,0:nchan-1],(eind-yind))
        endfor
        data  = newdata
        error = sqrt(data)
        totdetcounts = dblarr(ndet)
        for i = 0,ndet-1 do begin
          totdetcounts[i] = total(double(data[i,0:nchan-1]))
        endfor
        thisDetectorRates = dblarr(ndet)
        thisDetectorRates = totdetcounts/timeMinutes
  endif else begin
    energy = fltarr(nchan)
    energy = findgen(nchan)
  endelse

       END
'hy':  BEGIN
        void = dialog_message(['Hybrid mode not supported yet.'])
        RETURN
      END
ELSE:  void = dialog_message(['Not a valid HFBS file.'])
ENDCASE
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;if thisNorm2Mon eq 2 then begin	; normalize to total fission chamber counts
if thisNorm2Mon eq 999 then begin
  ndet = 17
  x = fltarr(nchan)
  x = energy
  wherefinite = where(data[0,*] gt 0.0,count)
  newdata = fltarr(ndet,count)
  newerror = fltarr(ndet,count)
  newdata[0:ndet-1,0:count-1] = data[0:ndet-1,wherefinite]
  newerror[0:ndet-1,0:count-1] = error[0:ndet-1,wherefinite]
  nchan = count
  y = fltarr(ndet,nchan)
  yerr = fltarr(ndet,nchan)
  for i = 1,16 do begin
    y[i,0:nchan-1] = newdata[i,0:nchan-1]/total(newdata[0,0:nchan-1])
    yerr[i,0:nchan-1] = newerror[i,0:nchan-1]/total(newdata[0,0:nchan-1])
  endfor
  data = fltarr(ndet,count)
  error = fltarr(ndet,count)
  data = y
  error = yerr
  data[0,*] = newdata[0,*]
  error[0,*] = newerror[0,*]

  energy = dblarr(nchan)
  energy = x[wherefinite]

endif
if thisNorm2Mon eq 1 then begin	; normalize to the beam monitor
  x = dblarr(nchan)
  x = double(energy)
  ndet = 17
  y = dblarr(ndet,nchan)
  yerr = dblarr(ndet,nchan)

  y = double(data[0:16,0:nchan-1])
  yerr = double(error[0:16,0:nchan-1])

  yMon = dblarr(nchan) & yMonerr = dblarr(nchan)
  yMon[0:nchan-1] = double(data[0,0:nchan-1])
  yMonerr = sqrt(yMon)

  wherefinite = where(ymon gt 0.0,count)
  yout = dblarr(ndet,count)
  xout = dblarr(count)
  yerrout = dblarr(ndet,count)
  yMonout = dblarr(count)
  yMonerrout = dblarr(count)

  xout[0:count-1] = x[wherefinite]
  yout[0:ndet-1,0:count-1] = y[0:ndet-1,wherefinite]
  yerrout[0:ndet-1,0:count-1] = yerr[0:ndet-1,wherefinite]
; Now perform the normalization to the monitor.  Instead of just subtracting
; off a flat background and pulling the spectrum down to zero, we assume that
; there is a non-zero background in both the monitor and the detectors.  A
; background is subtracted off of both measurements until the level in the
; extremes of the normalized data is within a few percent of each other.
; This bit of code added (08/26/02) by R.M.Dimeo.

; First cycle through detectors 5-14 and get an average value of
; the monitor background....note, only count those for which DAVE_AMOEBA
; converged!
;detlo = 5 & dethi = 14
;nrange = dethi-detlo+1
;counter = 0
;for j = detlo,dethi do begin
;  yin = data[j,wherefinite]
;  yinerr = error[j,wherefinite]
;  monitor = yMon[wherefinite]
;  xin = energy[wherefinite]
;  monerr = yMon[wherefinite]
;  normalizeData,xin,yin,yinerr,$
;                monitor,monerr,$
;                errmsg = errmsg,$
;                presult = presult,$
;                ynormout = ynormout,$
;                ynormerrout = ynormerrout
;  if presult[0] ne 1 then begin
;    if counter eq 0 then begin
;      barray = presult[1]
;    endif else begin
;      barray = [barray,presult[1]]
;    endelse
;    counter = counter+1
;  endif
;endfor
;bmon = (moment(barray))[0]
;print,'Average monitor background percentage: ',100.0*bmon/max(monitor)

;for j = 1,ndet-1 do begin
;  yin = data[j,wherefinite]
;  yinerr = error[j,wherefinite]
;  monitor = yMon[wherefinite]
;  xin = energy[wherefinite]
;  monerr = yMon[wherefinite]
;  if j lt 5 then maxiter = 1500 else maxiter = 500
;  if j eq 1 then pguess = [0.d,0.d]
;  normalizeData,xin,yin,yinerr,$
;                monitor,monerr,$
;                errmsg = errmsg,$
;                pguess = pguess, $
;                presult = presult,$
;                ynormout = ynormout,$
;                maxiter = maxiter, $
;                ynormerrout = ynormerrout
;  pguess = presult
;;  normalizeData_one,xin,yin,yinerr,$
;;                    monitor,monerr,$
;;                    bmon = bmon, $
;;                    errmsg = errmsg,$
;;                    presult = presult,$
;;                    ynormout = ynormout,$
;;                    ynormerrout = ynormerrout
;;  print,presult[0]/max(yin)
;  if errmsg eq 0 then begin
;    yout[j,*] = ynormout
;    yerrout[j,*] = ynormerrout
;  endif else begin
;    yout[j,*] = yin
;    yerrout[j,*] = yinerr
;  endelse
;endfor


;; Fit a straight line to the first few points and subtract this off from
;; the entire data set
  fewPts = 10
  for i = 1,ndet-1 do begin
    xline = fltarr(fewPts)
    yline = fltarr(fewPts)
    xline[0:fewPts-1] = xout[0:fewPts-1]
    yline[0:fewPts-1] = yout[i,0:fewPts-1]
    result = linfit(xline,yline)
    a = result[0] & b = result[1]
    ymin = b*xline[0]+a
    yout[i,*] = yout[i,*]-ymin
  endfor

;  yerrout[0:ndet-1,0:count-1] = yerr[0:ndet-1,wherefinite]
  yMonout[0:count-1] = yMon[wherefinite]
  yMonerrout[0:count-1] = yMonerr[wherefinite]

; Replace the loop with array operations for monitor normalization
  uvec = 1+bytarr(ndet-1)
  u = yout[1:ndet-1,0:count-1]
  v = uvec#yMonOut
  sigu = yerrout[1:ndet-1,0:count-1]
  sigv = uvec#yMonerrout

  z = u/v
  sigz = z*sqrt((sigu/u)^2+(sigv/v)^2)
  yout[1:ndet-1,0:count-1] = z
  yerrout[1:ndet-1,0:count-1] = sigz

  data = fltarr(ndet,count)
  data = yout
  error = fltarr(ndet,count)
  error = yerrout
  energy = fltarr(count)
  energy = xout
  nchan = count
endif
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
if strmid(header[2],21,2) eq 'ti' then begin
  ndet = 17
  data1 = fltarr(ndet,nchan)
  error1 = fltarr(ndet,nchan)
  data1[0:ndet-1,0:nchan-1] = data[0:ndet-1,0:nchan-1]
  error1[0:ndet-1,0:nchan-1] = error[0:ndet-1,0:nchan-1]
  data = fltarr(ndet,nchan) & data = data1
  error = fltarr(ndet,nchan) & error = error1

  integ_intensity = dblarr(ndet)
  ; scale peak maxima to 1.0
    for i = 0,ndet-1 do begin
  ;    print,'Det #',i
  ;   commented out the integration in the following line for a test (11/28/2000)
      integ_intensity[i] = int_tabulated(energy,data[i,*])

      if thisNorm2One eq 1 then begin
        scalefactor = max(data[i,0:nchan-1])+0.0*findgen(nchan)
        data[i,0:nchan-1] = data[i,0:nchan-1]/scalefactor
        error[i,0:nchan-1] = error[i,0:nchan-1]/scalefactor
      endif
    endfor

  integ_int = fltarr(ndet-1)
  integ_int[0:ndet-2] = integ_intensity[1:ndet-1]
  thisSumFC = total(data[0,*])
endif
for i = 1,ndet-1 do begin
  if thisNorm2One eq 1 then begin
    scalefactor = max(data[i,*])+0.0*findgen(nchan)
    data[i,0:*] = data[i,*]/scalefactor
    error[i,*] = error[i,*]/scalefactor
  endif
endfor
detRates = detRates
return
end

