;+
;NAME:
;	hxt_mcaccum
;PURPOSE:
; 	plots hxt data, prompts for start and end times for data 
;	accumulation and background, and returns an array with the
;       accumulated counts, for each interval and channel. Then
;       this will generate N sets of observations and sigmas
;       using poisson statistics. Restricted to 1 time interval.
;       Totals intervals before choosing samples, the old hxt_mcaccum
;       samples before totaling. 
;CALLING SEQUENCE:
;       Hxt_mcaccum, index, data, index1, hxi_index, counts, sigma, bcounts, $
;                     bsigma, trials=trials, img_trange=img_trange, $
;                     bck_trange=bck_trange, channel=channel, $
;                     accum_cnts=accum_cnts, ch_accum=ch_accum, $
;                     serr_only=serr_only, norm_error=norm_error, $
;                     no_serr=no_serr, st_times = st_times, $
;                     int_times = int_times
;INPUT:	
;	index=data info from hda file
;	data=counts in each subcollimator for each channel
;KEYWORDS:
;       trials = no. of MC trials for this interval. Default is 100
;       img_trange = a time range for the images, any yohkoh format
;       bck_range = a time_range for the background accumulation
;       channel= a channel number for the plots, default is to prompt the user
;       serr_only = if set, only use the systematic error term to
;                   obtain samples, the default is to add poisson
;                   variation and the systematic error term. The
;                   default systematic error is 0.03, and can be
;                   set using the sys_err keyword in HXT_MCIMG.
;       norm_error = if set, vary the SC counts using a normal distribution,
;                    with a standard deviation equal to norm_error*counts
;       no_serr = if set, do not add systematic error to the variation
;       So: counts = poidev(counts0)+sys_err*counts0*(2*(randomu(64,4)-0.5))
;           is the default variation. If norm_error is set,
;           then counts = counts0(1.0+norm_error(0)*randomn(64,4)+
;                         sys_err*counts0*(2*(randomu(64,4)-0.5)
;           If /serr_only is set, then
;           counts = sys_err*counts0*(2*(randomu(64, 4)-0.5)
;       And if the /no_serr keyword is set, the sys_err term is not included.
;       st_times= an array of accumulation start times, in any yohkoh format
;       int_times= interval times corresponding to st_times, may be a scalar,
;                  if this is not set, and st_times is set, the routine will
;                  use the differences between st_times. This must be an integer
;                  number of the appropriate resolution in seconds, if the
;                  value of int_times is less than a frame value.
;OUTPUT:
;       index1= index structure for the first frame of each image interval
;	hxi_index.chan= channel used
;		=0 L0 (15-24.4keV)
;		=1 M1 (24.4-35.2 keV)
;		=2 M2 (35.2 - 56.8 keV)
;		=3 HI (56.8 - 100.0 keV)
;	hxi_index.actim= image accumulation time in 10ths of secs
;	hxi_index.interval= actim
;	hxi_index.bkg_time= start time for background accumulation
;                           (millisec OF day)
;	hxi_index.bkg_day= day of start of background accumulation
;                          (days since 1-Jan-79)
;	hxi_index.bkg_actim= background accumulation time in 10ths of secs
;	hxi_index.bkg_interval= actim 
;       counts= the counts accumulated over the chosen intervals
;       sigma= sqrt(counts)
;       bcounts= the counts accumulated over the background interval
;       bsigma= sqrt(bcounts)
;DESCRIPTION:
;	Uses the function plot_lcur to select the data to be imaged,
;	and the background.
;MODIFICATION HISTORY:
;	Written Oct-93 by JMM,
;       Changed to give full resolution, 19-jul-94, JMM
;       9/21/94, account for 4.0 second time lag in hda files, jmm
;       1/9/95, changed calculation of sigma to be the same as
;               in Tom Metcalf's hxtaccumulate.pro, to be consistent
;               in pixon calculations.
;       1/13/95, added MC code here, as old hxt_mcimg is not valid
;       4/20/95, made consistent with new versions of hxtimg_accum, hxtsbt_mult,
;                And added channel and trange keywords, to enable batch runs
;       6/9/95, subtract background from counts before choosing intervals
;       4/15/96, Added serr_only, norm_error,
;       4/23/96, Added no_serr
;       1/6/97, Added st_times, int_times, noback
;-
PRO  Hxt_mcaccum, index, data, index1, hxi_index, counts, sigma, bcounts, $
                  bsigma, trials=trials, img_trange=img_trange, $
                  bck_trange=bck_trange, channel=channel, $
                  accum_cnts=accum_cnts, ch_accum=ch_accum, $
                  Serr_only=serr_only, Norm_error=norm_error, no_serr=no_serr, $
                  St_times=st_times, Int_times=int_times, noback=noback

   serr = hxi_index(0).syserr
   nc0 = 4                      ;no. of channels
   nc01 = nc0-1
   ndset = N_ELEMENTS(index)	;no. of data sets
   nds1 = ndset-1

;For full resolution, 4 minor frames for each data set
   nds4 = ndset*4 & nds41 = nds4-1

;No. of trials
   IF(KEYWORD_SET(trials)) THEN mct = fix(trials(0)) ELSE mct = 100

;i hate to have to do this, but here we decompress the whole data array
   data_x = float(hx_decomp(data))

;sum SC's to plot and you'll need a time array for each data point
   zz = hxt_chsums(index, data_x, /half_sec, tim_arr = tim_arr1, $
                   tim_sec = tim_arr, dt = dx4, /overwrite)
   time0 = tim_arr1(0)
   
;reform data to get into the proper form
   data_x = reform(data_x, nc0, 64, nds4, /overwrite)

;If i need to choose intervals by plotting
   no_time_selected = NOT (KEYWORD_SET(img_trange) OR KEYWORD_SET(st_times))
   no_bck_selected = NOT (KEYWORD_SET(bck_trange) OR KEYWORD_SET(noback))
   IF(no_time_selected OR no_bck_selected) THEN BEGIN
      IF(N_ELEMENTS(channel) GT 0) THEN ich = fix(channel(0)) ELSE BEGIN
         print, ' INPUT CHANNEL NUMBER FOR THE INTERVAL>>'
         print, ' 0 FOR LO CHANNEL'
         print, ' 1 FOR M1 CHANNEL'
         print, ' 2 FOR M2 CHANNEL'
         print, ' 3 FOR HI CHANNEL'
         read, ich
         ich = fix(ich)
      ENDELSE
   ENDIF ELSE ich = 0
   y0 = reform(zz(ich, *))/(dx4*64.0) ;y0 is the curve used to find the intervals
   IF(KEYWORD_SET(noback)) THEN BEGIN
      print, 'BACKGROUND SET TO ZERO'
      bcounts = fltarr(64, 4)
      bsigma = bcounts
      ib1 = 0 & ib2 = 3         ; i may need this stuff, for the structures
      ib14 = ib1/4 & ib24 = ib2/4
      bactim = total(dx4(ib1:ib2))
   ENDIF ELSE BEGIN      
;choose the background interval
      print, 'BACKGROUND SELECTION::'
      choose_interval, tim_arr1, ib1, ib2, data = y0, plot_title = 'HXT COUNTS/SECOND/SC', $
        time_range = bck_trange, dt_resolution = -1, alt_tsel = alt_tsel
      ib1 = ib1(0) & ib2 = ib2(0)

; Check time interval for non-flare mode data
      ib14 = ib1/4 & ib24 = ib2/4
      not_flare = where(gt_dp_mode(index(ib14:ib24)) NE 9)
      IF(not_flare(0) NE -1) THEN BEGIN
         message, /info, 'WARNING, NON-FLARE MODE BACKGROUND CHOSEN, IS ZERO FOR HI, M1, AND M2 CHANNELS'
      ENDIF
   
; time interval,
      bactim = total(dx4(ib1:ib2))
      tcounts = data_x(*, *, ib1:ib2)

;Speed things up, using total commands on 1 dimension
      nibs = ib2-ib1+1
      IF(nibs GT 1) THEN BEGIN
         bcounts = transpose(total(tcounts, 3))
         bsigma = transpose(sqrt(total(tcounts, 3)+ $
                                 total((serr*tcounts)^2, 3)))
      ENDIF ELSE BEGIN          ;no need to total
         bcounts = transpose(tcounts)
         bsigma = transpose(sqrt(tcounts+(serr*tcounts)^2)) ;oops, jmm,1-aug-95
      ENDELSE
   ENDELSE

;Ok, now you choose the image interval, subtract the background first
; you need a sum for the background
   IF(KEYWORD_SET(noback)) THEN bcps = fltarr(4) ELSE BEGIN
      IF(ib2 NE ib1) THEN bcps = total(zz(*, ib1:ib2), 2) ELSE bcps = zz(*, ib1)
   ENDELSE
   bcps = bcps/(64.0*bactim)    ;average background rate per SC

   print, 'IMAGE INTERVAL SELECTION::'
   IF(KEYWORD_SET(accum_cnts)) THEN BEGIN
      IF(accum_cnts EQ 1) THEN acc_cnts = 200.0 ELSE acc_cnts = accum_cnts(0)
      IF(N_ELEMENTS(ch_accum) GT 0) THEN ich0 = ch_accum(0) ELSE ich0 = 0
      y0 = reform(zz(ich0, *)/64.0)
      y0 = y0-bcps(ich0)*dx4
      choose_interval, tim_arr1, ss_st, ss_en, data = y0, accum_cnts = acc_cnts, $
        time_range = img_trange, plot_title = 'BACKGROUND SUBTRACTED HXT COUNTS/SC'
   ENDIF ELSE BEGIN
      choose_interval, tim_arr1, ss_st, ss_en, data = y0-bcps(ich), $
        plot_title = 'BACKGROUND SUBTRACTED HXT COUNTS/SEC/SC', $
        time_range = img_trange, dt_res = -1, st_times = st_times, int_times = int_times
   ENDELSE
   IF(N_ELEMENTS(ss_st) GT 1) THEN $
     print, 'ONLY ONE INTERVAL IS ALLOWED, THE FIRST INTERVAL CHOSEN IS USED'
   ss_st = ss_st(0)
   ss_en = ss_en(0)

; Check time interval for non-flare mode data
   ix14 = ss_st/4 & ix24 = ss_en/4
   not_flare = where(gt_dp_mode(index(ix14:ix24)) NE 9)
   IF(not_flare(0) NE -1) THEN BEGIN
      message, /info, 'WARNING, NON-FLARE MODE DATA IN CHOSEN INTERVAL IS ZERO FOR HI, M1, AND M2 CHANNELS'
   ENDIF

;You'll not have an index structure for each minor frame, so for
;each minor frame, find the appropriate index, and reset .time
;and .day
   ss_st_index = ss_st/4 & index1 = index(ss_st_index)
   index1.gen.time = tim_arr1(ss_st).time
   index1.gen.day = tim_arr1(ss_en).day
;And replicate index1
   index1 = replicate(index1(0), mct+1)

; now accumulate data, dti will be the time interval length, in seconds
   counts = fltarr(64, nc0, mct+1)
   sigma = counts

   dti = total(dx4(ss_st:ss_en))
   dt_intv = tim_arr(ss_en)-tim_arr(ss_st)+dx4(ss_en)
;go ahead and set the seed variable for the random generation,
;for some reason, poidev does not work in this routine without
;a seed being set.
   tcounts = data_x(*, *, ss_st:ss_en)
   dnn = ss_en-ss_st+1
   j0 = randomu(seed)
   FOR j = 0, mct DO BEGIN
      IF(dnn GT 1) THEN counts_x0 = transpose(total(tcounts, 3)) $
        ELSE counts_x0 = transpose(tcounts)
      IF(j GT 0) THEN BEGIN
         IF(KEYWORD_SET(no_serr)) THEN serr_term = 0.0 $
           ELSE serr_term = serr*counts_x0*(2.0*(randomu(seed, 64, nc0)-0.5))
         IF(KEYWORD_SET(norm_error)) THEN BEGIN
            counts(*, *, j)  = counts_x0*(1.0+norm_error*randomn(seed, 64, nc0))+serr_term
         ENDIF ELSE IF(KEYWORD_SET(serr_only)) THEN BEGIN
            counts(*, *, j) = counts_x0+serr_term
         ENDIF ELSE BEGIN
            counts(*, *, j) = poidev(counts_x0, seed = seed)+serr_term
         ENDELSE
      ENDIF ELSE counts(*, *, j) = counts_x0
   ENDFOR

   counts = counts > 0.0        ;gotta have it
   sigma = sqrt(counts+(serr*counts)^2)

; Replicate and fill hxi_index
   hxi_index.bkg_time = tim_arr1(ib1).time
   hxi_index.bkg_day = gt_day(tim_arr1(ib1))
   hxi_index.bkg_actim = fix(10.0*bactim)
   hxi_index.bkg_interval = fix(10.0*(tim_arr(ib2)-tim_arr(ib1)+dx4(ib2)))      
   hxi_index = replicate(hxi_index(0), nc0, mct+1)
   FOR i = 0, nc01 DO hxi_index(i, *).chan = i
   FOR j = 0, mct DO BEGIN
      hxi_index(*, j).actim = fix(10.0*dti)
      hxi_index(*, j).interval = fix(10.0*dt_intv)
   ENDFOR

   data_x = 0.0

; setgra gives the subscripts for the subcollimator switch
   setgra, sn2gra, gra2sn
   bcounts = bcounts(gra2sn, *)
   counts = counts(gra2sn, *, *)
   bsigma = bsigma(gra2sn, *)
   sigma = sigma(gra2sn, *, *)

   RETURN
END