;+
;PROCEDURE:	mvn_sta_bkg_load
;PURPOSE:	
;	Loads background data into APID c6 common block 
;INPUT:		
;
;KEYWORDS:
;	trange		dbl(2)	time range for background subtraction
;	AAA		flt	default=1., droop non-linear parameter - doesn't seem to matter that much

;  the following are primarily for testing purposes
;	no_update	0/1	when set, will not change background in any common block
;					this allows the routine to generate idlsave_d1_full_str without impacting bkg in common block
;	add		0/1	when set, this will not zero out background before adding the calculated background 
;					this allows backgrounds to be calculated on sequential runs as parameters are tuned
;	no_bkg3		0/1	when set, turns off bkg3 background removal - bkg3 not working yet
;	no_bkg4		0/1	when set, turns off bkg4 background removal
;	no_bkg5		0/1	when set, turns off bkg5 background removal
;	no_bkg6		0/1	when set, turns off bkg6 background removal
;	no_bkg7		0/1	when set, turns off bkg7 background removal
;	no_bkg8		0/1	when set, turns off bkg8 background removal
;	no_bkg9		0/1	when set, turns off bkg9 background removal
;	no_bkg10	0/1	when set, turns off bkg10 background removal
;	save_d1_full	0/1	when set, saves d1_full and bkg arrays in idlsave file - diagnostic for bkg3,bkg5 development - not recommended 
;
;
;CREATED BY:	J. McFadden
;VERSION:	1
;LAST MODIFICATION:  20/08/18		
;MOD HISTORY:
;
;NOTES:	  
; bkg2 is TBD - sputtered ions from the spacecraft caused by keV solar wind protons
;	this background is very tenuous, appears at 1-40eV, and at 20<mass<30 - most likely sputtered CO2++ and CO+
;	probably will not bother trying to develop this background estimation
; bkg3 is TBD - scattered ions off entrance posts and grids appearing at other energies and angles - proportional to rate
;	this is primarily observed for solar wind protons below the spectral peak - about 3% of the beam counts
;	the angle distribution appears to be complex so not clear whether this can be worked out
; bkg4 are internally backscattered cold ions near periapsis when the mechanical attenuator is closed
; 	removal assumes all heavy ions outside of anodes 5 to 9, and with energy < Esh, are background
;	Esh is the highest energy where the counts in c6 exceed 1/50. of the peak counts and exceed 20 counts.
;	the numbers 50 and 20 were empirically determined and should be tuned to make the algorithm work optimally 
;	tf_bkg4 is use to prevent double counting in background removal
; bkg5 is TBD (disabled) caused by event straggling (energy loss) in the carbon foils and from inelastic scattering from grids at the start foil
;	bkg5a are proton stragglers - proportional to event rate
; 	bkg5b are alpha  stragglers - proportional to event rate
;	these background are anode dependent - probably dependent on details of the grids
; bkg6 is TBD (disabled) background caused by start_no_stop followed closely by a valid event - proportional to rate squared
;	at high rate, the TOF circuit may not be completely reset before the next event producing a statistical high mass tail
;	bkg6 high mass tail is generally observed at masses just above the high flux mass peak (typically in solar wind)
;	however, valid events that contribute to bkg6 may also be due to other forms of background -- primarily coincidence events 
; bkg7 is background from coincident events where different ions (or radiation) produce the start and stop signals - proportional to rate squared
;	bkg7 has errors due to variations in flux on <4ms time scales which are corrected by mvn_sta_bkg_correct.pro
;	mvn_sta_bkg_correct.pro will also correct for the small variation in coincident background with mass
; bkg8 is TBD (disabled) background caused by start_no_stop followed closely by a TOF=0 event - proportional to rate squared
;	bkg8 is similar to bkg6
;	TOF=0 event are caused by scattered electrons that produce simultaneous start and stop signals
;	because bkg8 depends on rate^2, it requires the correction calculated by mvn_sta_bkg_correct.pro
; bkg9 is background created by delayed Start signal from scattered molecular ion fragment 
;	that failed to produce a start at the carbon foil but produced a Start electron when it hit another TOF analyzer surface
; bkg10 is TBD (disabled) bkg10 is internally produced background from sputtered ions generated by backscattered start-foil secondary electrons  
;	the primary sputtered ion peak is C++, however O+, O++ and C+ are also observed but are much smaller
;	the size of this peak depends on the incident ion mass (secondary e- production) and amount of disolved CO2 on surfaces
;	therefore the size of these peaks will vary with time, anode, and composition
;
; this algorithm currently generates bkg4, bkg7, bkg9
; 	bkg2 will likely never calculated since it is such a weak signal and depends on a complex dynamic of FOV, s/c orientation and solar wind flux
;	bkg3 may never be developed due depends on a complex dynamic of FOV, s/c orientation and solar wind flux
; 	bkg5,bkg6 will be treated in a separate program and use outputs from mvn_sta_bkg_correct.pro
;	bkg8 will be treated in a separate program using outputs from mvn_sta_bkg_correct.pro
;	bkg10 will be treated in a separate program using outputs from mvn_sta_bkg_correct.pro
;-

pro mvn_sta_bkg_load,trange=trange,AAA=AAA,add=add,no_update=no_update,no_bkg4=no_bkg4,no_bkg5=no_bkg5,no_bkg6=no_bkg6,no_bkg7=no_bkg7,no_bkg8=no_bkg8,no_bkg9=no_bkg9,save_d1_full=save_d1_full

starttime = systime(1)

common mvn_c0,mvn_c0_ind,mvn_c0_dat 
common mvn_c6,mvn_c6_ind,mvn_c6_dat 
common mvn_c8,mvn_c8_ind,mvn_c8_dat 
common mvn_ca,mvn_ca_ind,mvn_ca_dat 
common mvn_d0,mvn_d0_ind,mvn_d0_dat 
common mvn_d1,mvn_d1_ind,mvn_d1_dat 
common mvn_d8,mvn_d8_ind,mvn_d8_dat 
common mvn_d9,mvn_d9_ind,mvn_d9_dat 
common mvn_da,mvn_da_ind,mvn_da_dat 
common mvn_db,mvn_db_ind,mvn_db_dat 

; only require c6,db,d8,d9 data to be loaded

if size(mvn_c6_dat,/type) ne 8 then begin
	print,'Error - c6 data not loaded'
	return
endif
if size(mvn_db_dat,/type) ne 8 then begin
	print,'Error - db data not loaded'
	return
endif
if size(mvn_d9_dat,/type) ne 8 then begin
	print,'Error - d9 data not loaded'
	return
endif
if size(mvn_d8_dat,/type) ne 8 then begin
	print,'Error - d8 data not loaded'
	return
endif
if size(mvn_da_dat,/type) ne 8 then begin
	print,'Error - da data not loaded'
	return
endif

if size(mvn_d1_dat,/type) eq 8 then begin
	mtime_d1 = (mvn_d1_dat.time+mvn_d1_dat.end_time)/2.
	d1_loaded = 1 
endif else d1_loaded = 0

if size(mvn_d0_dat,/type) eq 8 and size(mvn_c8_dat,/type) eq 8 and size(mvn_ca_dat,/type) eq 8 then begin
	print,'Performing comprehensive background for all available data products'
	max_ind_c8=n_elements(mvn_c8_dat.time)-1
	max_ind_ca=n_elements(mvn_ca_dat.time)-1
	mtime_c6 = (mvn_c6_dat.time+mvn_c6_dat.end_time)/2.
	mtime_c0 = (mvn_c0_dat.time+mvn_c0_dat.end_time)/2.
	mtime_c8 = (mvn_c8_dat.time+mvn_c8_dat.end_time)/2.
	mtime_ca = (mvn_ca_dat.time+mvn_ca_dat.end_time)/2.
	mtime_d0 = (mvn_d0_dat.time+mvn_d0_dat.end_time)/2.
endif else begin
	print,'Aborting mvn_sta_bkg_load. All data products must be loaded: c0,c6,c8,ca,d0, - d1 is optional'
	print,'c0',size(mvn_c0_dat,/type)
	print,'c6',size(mvn_c6_dat,/type)
	print,'c8',size(mvn_c8_dat,/type)
	print,'ca',size(mvn_ca_dat,/type)
	print,'d0',size(mvn_d0_dat,/type)
	print,'d1',size(mvn_d1_dat,/type)
	return
endelse

; the following would have to be eliminated if another background program was to be run before this program

if not keyword_set(no_update) and not keyword_set(add) then begin
   if not keyword_set(trange) then begin
	mvn_c0_dat.bkg[*,*,*]=0
	mvn_c6_dat.bkg[*,*,*]=0
	mvn_c8_dat.bkg[*,*,*]=0
	mvn_ca_dat.bkg[*,*,*]=0
	mvn_d0_dat.bkg[*,*,*,*]=0
	if size(mvn_d1_dat,/type) eq 8 then mvn_d1_dat.bkg[*,*,*,*]=0
   endif else begin
	ind_jj=where(mvn_c6_dat.time gt trange[0] and mvn_c6_dat.time lt trange[1])
	mvn_c0_dat.bkg[ind_jj,*,*]=0
	mvn_c6_dat.bkg[ind_jj,*,*]=0
	mvn_c8_dat.bkg[ind_jj,*,*]=0
	mvn_ca_dat.bkg[ind_jj,*,*]=0
	mvn_d0_dat.bkg[ind_jj,*,*,*]=0
	if size(mvn_d1_dat,/type) eq 8 then mvn_d1_dat.bkg[ind_jj,*,*,*]=0
   endelse
endif

	nenergy = mvn_c6_dat.nenergy

;***************************************************************************************
; must run mvn_sta_dead_load,/test,/make_common for this to work

common mvn_sta_dead,dat_dead	

if size(dat_dead,/type) ne 8 then begin
	print,'Error - mvn_sta_dead common not created. RUN: mvn_sta_dead_load,/test,/make_common'
	return
endif

if max(dat_dead.rate) eq 0. then begin
	print,'Error - rate data indicates background is zero'
	return
endif

;	st_eff = replicate(1.,64)#replicate(.70,16) 						; start eff of protons, assume energy independent, initially assume anode/time independent but change it later
;	sp_eff = replicate(1.,64)#replicate(.47,16) 						; stop  eff of protons, assume energy independent, initially assume anode/time independent but change it later
;	st_no_stop_eff = replicate(1.,64)#replicate(.70,16) 					; start eff of protons, assume energy independent, initially assume anode/time independent but change it later

;	no_sp_eff = .45									; fraction of start_no_stop in rst - may need droop correction
;	sp_eff = .47									; fraction of start_no_stop in rst - may need droop correction
	integ_t = 0.0037988282								; accumulation time for a single sample, roughly 4s/1024.
	b_ns  = 5.844 									; tof bins per ns
	tof_offset = 16.

;   scale factor for bkg7 - the following accounts for anode variations in coincidence efficiency - empirically determined
;   this is likely due to signal degradation due to resistance in the 2ns delay lines in the anode array
	coin_an = replicate(1.,16)-0.25*(abs((findgen(16)-7.5))/7.5)
	scale = transpose(reform(replicate(1.,64l*16*64)#coin_an,64,16,64,16),[0,1,3,2])
;	scale = 1.0d


; Note on d8 = ['R1_Time_A','R1_Time_B','R1_Time_C','R1_Time_D','R1_Time_RST','R1_Time_NoStart',$
;		'R1_Time_Unqual','R1_Time_Qual','R1_Time_AnRej','R1_Time_MaRej','R1_Time_A&B','R1_Time_C&D']

; the following change allows mvn_sta_l2_tplot.pro to be run after this mvn_sta_c6_bkg_load.pro
;	previous code assumed mvn_sta_l2_tplot.pro had been run and used get_data to fill tmp1a to tmp7a
;	mvn_d8_dat has 4 sec time resolution

	tmp1a = {x:mvn_d8_dat.time,y:mvn_d8_dat.rates[*,7]/(mvn_d8_dat.rates[*,11]>1.e-10)}	; sp_eff, qual/C&D
	tmp2a = {x:mvn_d8_dat.time,y:mvn_d8_dat.rates[*,7]/(mvn_d8_dat.rates[*,10]>1.e-10)}	; st_eff, qual/A&B
	tmp3a = {x:mvn_d8_dat.time,y:mvn_d8_dat.rates[*,6]}					; unqual
	tmp4a = {x:mvn_d8_dat.time,y:mvn_d8_dat.rates[*,4]}					; rst - time_reset or total rate
	tmp5a = {x:mvn_d8_dat.time,y:mvn_d8_dat.rates[*,5]}					; no_start
	tmp6a = {x:mvn_d8_dat.time,y:mvn_d8_dat.rates[*,10]}					; A&B
	tmp7a = {x:mvn_d8_dat.time,y:mvn_d8_dat.rates[*,7]}					; qual

	get_data,'mvn_sta_dl_eff_qual',data=tmp8a						; this is made by mvn_sta_dead_load.pro



;***************************************************************************************
; the following is part of an inflight calibration

; the shift in the tof peak caused by imperfect delay lines
; 		 	 0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15
	mass_offset = [.05, .03, .05, .08, .00, .00, .00, .00, .15, .13, .18, .20, .20, .17, .20, .20]

;***************************************************************************************
; "stop anode efficiency"=sp_an_eff, determined from inflight data
; assumes this is constant for the mission and independent of count rate (droop only observed in start efficiency)
; this was measured in the sheath and/or solar wind

; 20160124 	anode 6,11
	sp_an_eff = [.415, .415, .440, .435, .420, .405, .405, .425, .460, .475, .420, .400, .435, .460, .440, .435]			; before anode reject table correction

;***************************************************************************************
; droop_rate is a measure of relative efficiency - end anodes have more droop due to losses in delay lines

;	anode number	0     1     2     3     4     5     6     7     8     9    10 	 11    12    13    14    15
; 20160401 	anode 0,1,2,7,12,13,15
	droop_rate = [1.00, 1.20, 1.00, 1.00, 1.00, 1.00, 1.00, 1.20, 1.00, 1.00, 1.00, 1.00, 1.40, 1.00, 1.00, 1.00]*1.		; anode-11 assume same as 10
	droop_rate = [1.22, 1.25, 1.25, 1.00, 1.00, 1.00, 1.00, 1.45, 1.00, 1.00, 1.00, 1.00, 1.40, 1.10, 1.00, 1.15]*1.		; anode-11 assume same as 10
; 20180106 	anode 1,(5,6),8,9,13,15
	droop_rate = [1.22, 1.25, 1.25, 1.00, 1.00, 1.40, 1.40, 1.45, 1.00, 1.30, 1.00, 1.00, 1.40, 1.10, 1.00, 1.15]*1.		; anode-11 assume same as 10
; 20171223 	anode 4,6,7,12,13,15
	droop_rate = [1.22, 1.25, 1.25, 1.00, 1.00, 1.40, 1.40, 1.45, 1.00, 1.30, 1.00, 1.00, 1.40, 1.10, 1.00, 1.05]*1.		; anode-11 assume same as 10
; 20180113 	anode 8,9,
	droop_rate = [1.22, 1.25, 1.25, 1.00, 1.00, 1.40, 1.40, 1.45, 1.20, 1.25, 1.00, 1.00, 1.40, 1.10, 1.00, 1.15]*1.		; anode-11 assume same as 10
; 20160303 	anode 4,15
	droop_rate = [1.22, 1.25, 1.25, 1.00, 1.38, 1.40, 1.40, 1.45, 1.20, 1.25, 1.00, 1.00, 1.40, 1.10, 1.00, 1.05]*1.		; anode-11 assume same as 10
; 20160313 	anode 4,15
	droop_rate = [1.22, 1.25, 1.25, 1.10, 1.38, 1.40, 1.40, 1.45, 1.20, 1.25, 1.00, 1.00, 1.40, 1.10, 1.00, 1.05]*1.		; anode-11 assume same as 10
; 20170623 	anode 0,10,11,15
	droop_rate = [1.22, 1.25, 1.25, 1.10, 1.38, 1.40, 1.40, 1.45, 1.20, 1.25, 1.15, 1.40, 1.40, 1.10, 1.00, 1.05]*1.		; anode-11 assume same as 10


;***************************************************************************************
; Empirically determined formula for use in calculating bkg9
; Molecular ions produce a low mass shoulder due to a delayed TOF Start 
; For some events the molecular ion fails to produce a Start event as it exits the foil 
; The molecule fragments at the foil with one fragment hitting an internal TOF analyzer surface producing a delayed Start electron
; The other molecular fragment produces a normal Stop signal
; The shoulder shifts its TOF profile with the shift of the main TOF peak (as the Acceleration potential is changed).
; Assume the low mass shoulder TOF profile is proportional to the TOF peak (ie an effective shorter TOF for these background events)
; tof2 is a mapping of valid tof events to shoulder tof events determined empirically
; bkg9 is only important at low energy where O2 and CO2 fluxes are large

tof1 = reform(mvn_c6_dat.tof_arr[1,31,*])		; low energy tof array
tof1 = (tof1+tof_offset)/b_ns				; convert to nanosec
tof2 = fltarr(64,64)
;fra2 = .49
fra2 = .48
for i=0,63 do begin
	for j=0,63 do begin
		tof2[i,j] = (tof1[j] gt (1.-fra2)*tof1[i])*((tof1[i]-tof1[j])>0.)/(fra2*tof1[i])
	endfor
endfor

;***************************************************************************************
; assume c6 exists for all times - interpolate all other events to this time base

npts = n_elements(mvn_c6_dat.time)

if keyword_set(trange) then begin
	ind_jj=where(mvn_c6_dat.time gt trange[0] and mvn_c6_dat.time lt trange[1])
	npts=n_elements(ind_jj)
	datetime = time_string(trange[0])
	yrmodatm = strmid(datetime,0,4)+strmid(datetime,5,2)+strmid(datetime,8,2)+'_'+strmid(datetime,11,2)+strmid(datetime,14,2)+strmid(datetime,17,2)
endif else begin
	ind_jj=lindgen(npts)
	datetime = time_string(mvn_c6_dat.time[0])
	yrmodatm = strmid(datetime,0,4)+strmid(datetime,5,2)+strmid(datetime,8,2)+'_'+strmid(datetime,11,2)+strmid(datetime,14,2)+strmid(datetime,17,2)
endelse

; these are for tracking errors in interpolating d1,ca into a complete 64Ex16Dx16A array for conincidence background - bkg7
	n_ca=0
	n_d1_4s = 0
	n_d1_16s = 0
	n_ca_4s = 0

	if not keyword_set(AAA) then AAA = 1.			; this is a controlling factor for non-linear droop in bkg7, probably not very sensitive to values

; these are for testing purposes
	st_to_rate_min = fltarr(npts)
	rate_max = fltarr(npts)
	st_to_rate_diff = fltarr(npts)
	rate0_arr = fltarr(npts)
	n_iter_arr = fltarr(npts)

;***************************************************************************************
;***************************************************************************************
; the following is for test purposes in developing bkg5,bkg3

if keyword_set(save_d1_full) then begin

; the following arrays were too large to work with
;	d1_full_arr = fltarr(npts,32,16,16,64)
;	d1_full_bkg = fltarr(npts,32,16,16,64)
;	d1_full_tim = dblarr(npts)

	d1_full_ct1 = fltarr(npts,32,16,64)
	d1_full_bg1 = fltarr(npts,32,16,64)
	d1_full_ct2 = fltarr(npts,32,16,16,8)
	d1_full_bg2 = fltarr(npts,32,16,16,8)
	d1_full_tim = dblarr(npts)
	d1_full_end = dblarr(npts)
	d1_full_nrg = fltarr(npts,32,64)
	d1_full_twt = fltarr(npts,32,64)
	d1_full_tof = fltarr(npts,32,64)
	d1_full_mas = fltarr(npts,32,64)
endif

; get dl background values from tplot variables - these are rates
;	dl_bkg_ab, dl_bkg_cd are determined from apid d9 rates in mvn_sta_dead_load.pro
;		ab_bk = total(mvn_d9_dat.rates[ind_d9,10,indbk9_ab[0:9]])/10. > 25.
;		cd_bk = total(mvn_d9_dat.rates[ind_d9,11,indbk9_cd[0:9]])/10. > 90.

	get_data,'mvn_sta_dl_bkg_ab',data=dl_bkg_ab_arr
	get_data,'mvn_sta_dl_bkg_cd',data=dl_bkg_cd_arr
	get_data,'mvn_sta_dl_bkg_fq',data=dl_bkg_fq_arr

	if size(dl_bkg_cd_arr,/type) ne 8 then begin
		print,'Error: dl_bkg_cd_arr values are zero'
		print,'    Must run:   mvn_sta_dead_load,/test,/make_common'
		return
	endif
	if max(dl_bkg_cd_arr.y) le 10 then begin
		print,'Error: dl_bkg_cd_arr values are zero'
		print,'    Must run:   mvn_sta_dead_load,/test,/make_common'
		return
	endif


; make da_bkg_sm avegeraged over 11 points (40s)
; 	these are rates - used by bkg6

	nn_da = n_elements(mvn_da_dat.time)
	da_bkg = fltarr(nn_da)
	for pp=0,nn_da-1 do begin
		da_ind = sort(mvn_da_dat.rates[pp,*])
		da_bkg[pp]=total(mvn_da_dat.rates[pp,da_ind[0:9]])/10.
	endfor	
	store_data,'da_bkg',data={x:(mvn_da_dat.time+mvn_da_dat.end_time)/2.,y:da_bkg}
	tsmooth2,'da_bkg',11
	get_data,'da_bkg_sm',data=da_bkg_sm


;***************************************************************************************
;***************************************************************************************
;***************************************************************************************
;***************************************************************************************
; main loop starts here

print,'Starting main loop'


for jj=0l,npts-1 do begin
     ii = ind_jj[jj]
	if jj mod 1000 eq 1 then print,jj,ii,npts,systime(1)-starttime		; tracking progress

	alt = (total(mvn_c6_dat.pos_sc_mso[ii,*]^2))^.5-3386.
	att = mvn_c6_dat.att_ind[ii]
	mode = mvn_c6_dat.mode[ii]

; get nearest deadtime and droop from common mvn_sta_dead, dat_dead  -- must have run: mvn_sta_dead_load,/test,/make_common

	nearest = min(abs(dat_dead.time - mvn_c6_dat.time[ii]),ind)

	rate = reform(dat_dead.rate[ind,*,*])			; dblarr(npts,64,16)  energy-deflector event rate		
	valid = reform(dat_dead.valid[ind,*,*])			; dblarr(npts,64,16)  energy-deflector valid counts				
	droop = reform(dat_dead.droop[ind,*,*])			; dblarr(npts,64,16)  energy-deflector droop correction			
	dead = reform(dat_dead.dead[ind,*,*])			; dblarr(npts,64,16)  energy-deflector dead correction		
	anode = reform(dat_dead.anode[ind,*,*])			; dblarr(npts,64,16)  energy-anode distribution of counts, normalized at each energy - assumes anode distribution does not depend on deflector or mass
	rate_valid_ratio = (rate/1024.+5.)/(((rate/1024.)<valid)+5.)
	rate_d = rate*dead

; get dl background values from tplot variables - used by bkg6

	nearest = min(abs(dl_bkg_ab_arr.x - mvn_c6_dat.time[ii]),ind_tmp)
	dl_bkg_ab = dl_bkg_ab_arr.y[ind_tmp]
	dl_bkg_cd = dl_bkg_cd_arr.y[ind_tmp]
	dl_bkg_fq = dl_bkg_fq_arr.y[ind_tmp]

; get da rst background

	nearest = min(abs(da_bkg_sm.x -2. - mvn_c6_dat.time[ii]),ind_tmp)
	da_bk = da_bkg_sm.y[ind_tmp]

; get nearest d8 values from common mvn_d8,mvn_d8_ind,mvn_d8_dat
; Rates: 0:TA, 1:TB, 2:TC, 3:TD, 4:Trst, 5:nostart, 6:unqual, 7:qual, 8:ano_rej, 9:mass_rej, 10:A&B, 11:C&D, 

	nearest = min(abs(mvn_d8_dat.time - mvn_c6_dat.time[ii]),ind_d8)

	d8_arr = reform(mvn_d8_dat.rates[ind_d8,*])			; dblarr(npts,64,16)  energy-deflector event rate
	d8_ta = d8_arr[0]		
	d8_tb = d8_arr[1]		
	d8_tc = d8_arr[2]		
	d8_td = d8_arr[3]		
	d8_trst = d8_arr[4]		
	d8_nostart = d8_arr[5]		
	d8_unqual = d8_arr[6]		
	d8_qual = d8_arr[7]		
	d8_anrej = d8_arr[8]		
	d8_marej = d8_arr[9]		
	d8_ab = d8_arr[10]		
	d8_cd = d8_arr[11]		

;	m_rate=0.5e5
	def_st=0.85

;	st_to_rst = 0.2<((d8_ta>d8_tb)/(d8_trst))<def_st					; 0.2 is a guess, probably should be higher
;	st_to_rst = 0.2>(((d8_ta>d8_tb)-dl_bkg_ab)/((d8_trst-da_bk)>1))<def_st			; 0.2 is a guess, probably should be higher
	st_to_rst = 0.2>(((d8_ta<d8_tb)-dl_bkg_ab)/((d8_trst-da_bk)>1))<def_st			; 0.2 is a guess, probably should be higher
		
; determine start_to_rate for use by bkg6

		width0 = 1.0
		rate0 = 200000.
;		min_rate=10000.
		min_rate=1.

; poisson
		tof_grid = 0.65				; tof grid transmission
		thick_foil_gain = 2.			; unknown - contributes to mcp_eff_stop
		thick_foil = .9				; thick foils transmission of 11 keV electrons
		foil_coverage_start = 0.95
		foil_coverage_stop  = 0.90*tof_grid	; stop foils have less coverage
		lamda=2.13  
		pois1_fraction=0.90	
		mcp_eff=.60
		mcp_eff_stop=0.81
		mcp_eff2 = mcp_eff_stop*thick_foil

		kk=0 & pois0 = lamda^kk*exp(-lamda)/factorial(kk)
		kk=1 & pois1 = lamda^kk*exp(-lamda)/factorial(kk)
		kk=2 & pois2 = lamda^kk*exp(-lamda)/factorial(kk)
		kk=3 & pois3 = lamda^kk*exp(-lamda)/factorial(kk)
		kk=4 & pois4 = lamda^kk*exp(-lamda)/factorial(kk)
		kk=5 & pois5 = lamda^kk*exp(-lamda)/factorial(kk)
		kk=6 & pois6 = lamda^kk*exp(-lamda)/factorial(kk)

	pois1=pois1_fraction*pois1
	start_eff = foil_coverage_start*(pois1*mcp_eff+pois2*(1-(1-mcp_eff)^2)+pois3*(1-(1-mcp_eff)^3)+pois4*(1-(1-mcp_eff)^4)+pois5*(1-(1-mcp_eff)^5)+pois6*(1-(1-mcp_eff)^6))
	stop_eff  = foil_coverage_stop *(pois1*mcp_eff2+pois2*(1-(1-mcp_eff2)^2)+pois3*(1-(1-mcp_eff2)^3)+pois4*(1-(1-mcp_eff2)^4)+pois5*(1-(1-mcp_eff2)^5)+pois6*(1-(1-mcp_eff2)^6))

; print,start_eff, stop_eff

		start_to_rate = st_to_rst*replicate(1.,64,16)
;		tmp71 = total(start_to_rate*rate)/total(rate)
		tmp71 = total(def_st*rate)/total(rate)
		n_iter = 0

; print,n_iter,(tmp71-st_to_rst),st_to_rst,rate0,max(rate),minmax(start_to_rate)

    if st_to_rst lt def_st then begin

	while (abs(tmp71-st_to_rst) gt 0.001) and (n_iter lt 15) do begin

		zz1 = width0*(rate_d/(1*rate0) - 1./(1.001-exp(-rate/min_rate)) )
		zz2 = width0*(rate_d/(2*rate0) - 1./(1.001-exp(-rate/min_rate)) )
		zz3 = width0*(rate_d/(3*rate0) - 1./(1.001-exp(-rate/min_rate)) )
		zz4 = width0*(rate_d/(4*rate0) - 1./(1.001-exp(-rate/min_rate)) )
		zz5 = width0*(rate_d/(5*rate0) - 1./(1.001-exp(-rate/min_rate)) )
		zz6 = width0*(rate_d/(6*rate0) - 1./(1.001-exp(-rate/min_rate)) )

		start_eff = foil_coverage_start*( pois1*mcp_eff *(1.-erf(zz1))/2.$
			+pois2*(1-(1-mcp_eff)^2)*(1.-erf(zz2))/2. $
			+pois3*(1-(1-mcp_eff)^3)*(1.-erf(zz3))/2. $
			+pois4*(1-(1-mcp_eff)^4)*(1.-erf(zz4))/2. $
			+pois5*(1-(1-mcp_eff)^5)*(1.-erf(zz5))/2. $
			+pois6*(1-(1-mcp_eff)^6)*(1.-erf(zz6))/2. $
			)

		stop_eff = .45

		start_to_rate = start_eff/(start_eff + (1-start_eff)*stop_eff)

		tmp71 = total(start_to_rate*rate)/total(rate)
		rate0 = rate0*(1. - (-.5>(5.0*(tmp71-st_to_rst))<.5))
		n_iter = n_iter+1

	endwhile

    endif else print,'st_to_rst>def_st, st_to_rst=',st_to_rst,time_string(mvn_c6_dat.time[ii])

	st_to_rate_min[jj] = min(start_to_rate)
	rate_max[jj] = max(rate)
	st_to_rate_diff[jj]=tmp71-st_to_rst
	rate0_arr[jj]=rate0
	n_iter_arr[jj] = n_iter
;	print,n_iter,(tmp71-st_to_rst),st_to_rst,rate0,max(rate),minmax(start_to_rate)

; anode_sq is needed for coincident events
	anode_sq = total(anode^2,2)#replicate(1.,16)	

; the following are used to estimate the background rate from penetrating radiation - needed for SEP events
; use the 100 lowest count rate, energy-def bin out of 1024, averaged over three 4sec periods
; brate is a diagnostic
					
	ratem = reform(dat_dead.rate[(ind-1)>0,*,*])						
	ratep = reform(dat_dead.rate[(ind+1)<npts-1,*,*])
	sort_rate0 = rate[sort(rate)]					
	brate0 = total(sort_rate0[0:99])			
	sort_ratem = ratem[sort(ratem)]
	bratem = total(sort_ratem[0:99])		
	sort_ratep = ratep[sort(ratep)]
	bratep = total(sort_ratep[0:99])	
	brate = (brate0+bratem+bratep)/300.		


	mlut_ind = mvn_c6_dat.mlut_ind[ii]
	swp_ind = mvn_c6_dat.swp_ind[ii]
	tof_arr = reform(mvn_c6_dat.tof_arr[mlut_ind,*,*])
	twt_arr = reform(mvn_c6_dat.twt_arr[mlut_ind,*,*])

; tmp1a-tmp7a are from product d8 with 4 sec cadence
; tmp8a is tplot str 'mvn_sta_dl_eff_qual' created by mvn_sta_dead_load.pro

	nearest = min(abs(tmp1a.x - (mvn_c6_dat.time[ii]+2.)),ind_d8)

	nearest = min(abs(tmp8a.x - (mvn_c6_dat.time[ii]+2.)),ind_dl)

	tmp1 = tmp1a.y[ind_d8]		; st_eff
	tmp2 = tmp2a.y[ind_d8]		; sp_eff
	tmp3 = tmp3a.y[ind_d8]		; unqual  (rate)
	tmp4 = tmp4a.y[ind_d8]		; rst or  (total rate)
	tmp5 = tmp5a.y[ind_d8]		; nostart (rate)
	tmp6 = tmp6a.y[ind_d8]		; A&B	  (rate)
	tmp7 = tmp7a.y[ind_d8]		; qual    (rate)
	tmp8 = tmp8a.y[ind_dl]		; qu_eff  (qualified eff from mvn_sta_dead_load.pro)
; 	tmp9 = see below		; qu_eff  (the qualified efficiency from d8 product using d9 for background subtractions)

;***************************************************************************************
; tmp9 is the average qualified efficiency from d8 product using d9 product for background - generally a better estimate than tmp8 
; 	valid is a 64x16 energy-deflector array determined by mvn_sta_dead_load.pro
; 	code uses the larger of tmp8,tmp9:  qu_eff[*,*] = tmp9>tmp8

; db is the mass histogram rate - some qualified TOF events are outside the range allowed into data products
	nearest = min(abs((mvn_db_dat.time+mvn_db_dat.end_time)/2. - (mvn_c6_dat.time[ii]+2.)),ind_db)
	db_arr = reform(mvn_db_dat.data[ind_db,*])
	qu_valid = total(db_arr[42:815])/total(db_arr)				; not all valid events are sent to data products

; d9 is the energy dependent version of d8 averaged over 128sec periods 
; this is used to get background as needed for SEP events, or after-emission background after detector saturation
	nearest = min(abs((mvn_d9_dat.time+mvn_d9_dat.end_time)/2.-(mvn_c6_dat.time[ii]+2.)),ind_d9)
	indbk9 = sort(mvn_d9_dat.rates[ind_d9,7,*])
	bk7 = total(mvn_d9_dat.rates[ind_d9,7,indbk9[0:9]])/10.			; avg background rate from 10 lowest rate channels in d9, needed for high bkg SEP events

; tmp9 estimate of qu_eff limited to 1/1.1 ~ 0.9 to prevent statistical fluctuation errors when bk7~tmp7
; 	4.* converts (tmp7-bk7) rate to counts
	tmp9 = total(valid)/((4.*qu_valid*(tmp7-bk7)) > 1.1*total(valid))	 

; these lines were for testing - expect tmp9 to be >tmp8, identify those times where this is not the case
;	if (tmp8 gt tmp9) then print,tmp8,tmp9,qu_valid,total(valid),tmp7,bk7,brate,total(rate*integ_t)
; 	if (tmp8 gt tmp9) then print,total(rate)/1024.,total(sort_rate0[0:1023])*integ_t,total(sort_rate0[0:500])*integ_t,total(sort_rate0[0:99])*integ_t
	
;***************************************************************************************
;***************************************************************************************

; use d1 data product if loaded


if d1_loaded then begin

	cnts_c6 = reform(mvn_c6_dat.data[ii,*,*])
	time_c6 = mtime_c6(ii)

	nearest = min(abs(mtime_c0-time_c6),ii_c0)
	cnts_c0 = reform(mvn_c0_dat.data[ii_c0,*,*])
	time_c0 = mtime_c0[ii_c0] & delta_t_c0 = mvn_c0_dat.delta_t[ii_c0]

	nearest = min(abs(mtime_d1-time_c6),ii_d1)
	cnts_d1 = reform(mvn_d1_dat.data[ii_d1,*,*,*])
	time_d1 = mtime_d1[ii_d1] & delta_t_d1 = mvn_d1_dat.delta_t[ii_d1]

	nearest = min(abs(mtime_d0-time_c6),ii_d0)
;	cnts_d0 = reform(mvn_d0_dat.data[ii_d0,*,*,*])
	time_d0 = mtime_d0[ii_d0] & delta_t_d0 = mvn_d0_dat.delta_t[ii_d0]

	nearest = min(abs(mtime_c8-time_c6),ii_c8)
	cnts_c8 = reform(mvn_c8_dat.data[ii_c8,*,*])
	time_c8 = mtime_c8[ii_c8] 
	if (time_c8-time_c6) gt 3. then cnts_c8 = (cnts_c8+reform(mvn_c8_dat.data[0>(ii_c8-1),*,*]))/2.
	if (time_c8-time_c6) lt -3 then cnts_c8 = (cnts_c8+reform(mvn_c8_dat.data[max_ind_c8<(ii_c8+1),*,*]))/2.

	nearest = min(abs(mtime_ca-time_c6),ii_ca)
	cnts_ca = reform(mvn_ca_dat.data[ii_ca,*,*])
	time_ca = mtime_ca[ii_ca] 
	if (time_ca-time_c6) gt 3. then cnts_ca = (cnts_ca+reform(mvn_ca_dat.data[0>(ii_ca-1),*,*]))/2.
	if (time_ca-time_c6) lt -3 then cnts_ca = (cnts_ca+reform(mvn_ca_dat.data[max_ind_ca<(ii_ca+1),*,*]))/2.

	if (delta_t_d1 le 5.) and (abs(time_d1-time_c6) lt 2. or (abs(time_d1-time_c6) lt 5. and abs(total(cnts_d1)-total(cnts_c6)) lt total(cnts_c6)^.5)) then begin			; use this when d1 is at c6 cadence

	     if abs(time_d1-time_c8) lt 2. and abs(total(cnts_d1)-total(cnts_c8)) lt .5*total(cnts_d1)^.5 then begin
		norm_c8  = transpose(reform(replicate(1.,4)#reform(total(reform(transpose(cnts_c8),4,4,32),1),4*32),16,32))
		def4to16 = reform(reform(cnts_c8/(norm_c8+.00000001),32*16)#replicate(1.,16*64),32,16,16,64)
	     endif else begin
		norm_c8  = transpose(reform(replicate(1.,4)#reform(total(reform(transpose(cnts_c8>1.),4,4,32),1),4*32),16,32))
		def4to16 = reform(reform((cnts_c8>1.)/(norm_c8),32*16)#replicate(1.,16*64),32,16,16,64)
	     endelse

	     if abs(time_d1-time_c6) lt 2. and abs(total(cnts_d1)-total(cnts_c6)) lt .5*total(cnts_d1)^.5 then begin
		norm_c6  = reform(transpose(reform(reform(total(reform(cnts_c6,32,8,8),2),32*8)#replicate(1.,8),32,8,8),[0,2,1]),32,64)
		mas8to64 = transpose(reform(reform(cnts_c6/(norm_c6+.00000001),32*64)#replicate(1.,16*16),32,64,16,16),[0,2,3,1])
	     endif else begin
		norm_c6  = reform(transpose(reform(reform(total(reform((cnts_c6>.1),32,8,8),2),32*8)#replicate(1.,8),32,8,8),[0,2,1]),32,64)
		mas8to64 = transpose(reform(reform((cnts_c6>.1)/(norm_c6),32*64)#replicate(1.,16*16),32,64,16,16),[0,2,3,1])
	     endelse

		d1_def = transpose(reform(replicate(1.,4)#reform(transpose(reform(cnts_d1,32,4,16,8),[1,0,2,3]),4*32*16*8),16,32,16,8),[1,0,2,3])
		d1_tmp = reform(transpose(reform(reform(d1_def,32l*16*16*8)#replicate(1.,8),32,16,16,8,8),[0,1,2,4,3]),32,16,16,64)

		d1_full = def4to16*d1_tmp*mas8to64				; 32Ex16Dx16Ax64M

		if abs(total(cnts_d1) - total(d1_full)) gt 1. then begin
			print,'Error in 4sec d1 normalization, times differ by:',time_d1-time_c6,' at ',time_string(time_c6),time_c8-time_c6
			print,'Normalization changes counts by:',abs(total(cnts_d1)-total(d1_full)),total(d1_full),total(cnts_d1),total(cnts_c6),total(cnts_c8)
			print,'d1_tmp   ',abs(total(cnts_d1)-total(d1_tmp/32.))
			print,'def4to16 ',abs(total(cnts_d1)-total(def4to16*d1_tmp/8.))
			print,'mas8to64 ',abs(total(cnts_d1)-total(mas8to64*d1_tmp/4.))
		endif

		d1_full = d1_full*(total(cnts_c6)/total(cnts_d1))
		n_d1_4s = n_d1_4s + 1

	endif else if (delta_t_d1 gt 5. and delta_t_d1 lt 17.) and (abs(time_d1-time_c6) lt 8.) and (abs(total(cnts_d1)/4.-total(cnts_c6)) lt 2.*total(cnts_c6)^.5) then begin

		norm_c8  = transpose(reform(replicate(1.,4)#reform(total(reform(transpose(cnts_c8>1.),4,4,32),1),4*32),16,32))
		def4to16 = reform(reform((cnts_c8>1.)/(norm_c8),32*16)#replicate(1.,16*64),32,16,16,64)

		norm_c6  = reform(transpose(reform(reform(total(reform((cnts_c6>.1),32,8,8),2),32*8)#replicate(1.,8),32,8,8),[0,2,1]),32,64)
		mas8to64 = transpose(reform(reform((cnts_c6>.1)/(norm_c6),32*64)#replicate(1.,16*16),32,64,16,16),[0,2,3,1])

		d1_def = transpose(reform(replicate(1.,4)#reform(transpose(reform(cnts_d1,32,4,16,8),[1,0,2,3]),4*32*16*8),16,32,16,8),[1,0,2,3])
		d1_tmp = reform(transpose(reform(reform(d1_def,32l*16*16*8)#replicate(1.,8),32,16,16,8,8),[0,1,2,4,3]),32,16,16,64)

		d1_full = def4to16*d1_tmp*mas8to64		; 32Ex16Dx16Ax64M

		if abs(total(cnts_d1) - total(d1_full)) gt 1. then begin
			print,'Error in 16sec d1 normalization, times differ by:',time_d1-time_c6,' at ',time_string(time_c6)
			print,'Normalization changes counts by:',abs(total(cnts_d1)-total(d1_full)),total(d1_full),total(cnts_d1),total(cnts_c8),total(cnts_c6)*4.,total(cnts_c6)^.5
			print,'d1_tmp   ',abs(total(cnts_d1)-total(d1_tmp/32.))
			print,'def4to16 ',abs(total(cnts_d1)-total(def4to16*d1_tmp/8.))
			print,'mas8to64 ',abs(total(cnts_d1)-total(mas8to64*d1_tmp/4.))
		endif

		d1_full = d1_full*(total(cnts_c6)/total(cnts_d1))		; 32Ex16Dx16Ax64M
		n_d1_16s = n_d1_16s + 1

	endif else if abs(time_ca - time_c6) lt 6. then begin		; use ca to construct d1_full 

;	     if abs(time_ca-time_c8) lt 2. and abs(total(cnts_ca)-total(cnts_c8)) lt .5*total(cnts_ca)^.5 then begin
;		norm_c8  = transpose(reform(replicate(1.,4)#reform(total(reform(transpose(cnts_c8),4,4,32),1),4*32),16,32))
;		def4to16 = reform(reform(cnts_c8/(norm_c8+.00000001),32*16)#replicate(1.,16*64),32,16,16,64)
tfdef=0
;	     endif else begin
		norm_c8  = transpose(reform(replicate(1.,4)#reform(total(reform(transpose(cnts_c8>1.),4,4,32),1),4*32),16,32))
		def4to16 = reform(reform((cnts_c8>1.)/(norm_c8),32*16)#replicate(1.,16*64),32,16,16,64)
tfdef=1
;	     endelse

	     if abs(time_ca-time_c6) lt 2. and abs(total(cnts_ca)-total(cnts_c6)) lt .5*total(cnts_ca)^.5 then begin
		norm_c6  = reform(total(cnts_c6,2)#replicate(1.,64),32,64)
		mas1to64 = transpose(reform(reform(cnts_c6/(norm_c6+.00000001),32*64)#replicate(1.,16*16),32,64,16,16),[0,2,3,1])
tfmas=0
	     endif else begin
		norm_c6  = reform(total((cnts_c6>.1),2)#replicate(1.,64),32,64)
		mas1to64 = transpose(reform(reform((cnts_c6>.1)/(norm_c6),32*64)#replicate(1.,16*16),32,64,16,16),[0,2,3,1])
tfmas=1
	     endelse

	     if abs(time_ca-time_c6) lt 2. and abs(total(cnts_ca)-total(cnts_c6)) lt .5*total(cnts_ca)^.5 then begin
		en16to32 = reform((total(cnts_c6,2)/(reform(replicate(1.,2)#total(reform(total(cnts_c6,2),2,16),1),32)+.0000001))#replicate(1.,64),32,64)
		ca_32E = en16to32*reform(replicate(1.,2)#reform(cnts_ca,16*64),32,64)
tfe32=0
	     endif else begin
		en16to32 = reform((total((cnts_c6>1.),2)/(reform(replicate(1.,2)#total(reform(total((cnts_c6>1.),2),2,16),1),32)))#replicate(1.,64),32,64)
		ca_32E = en16to32*reform(replicate(1.,2)#reform(cnts_ca,16*64),32,64)
tfe32=1
	     endelse

		ca_def = transpose(reform(replicate(1.,4)#reform(transpose(reform(ca_32E,32,4,16),[1,0,2]),4*32*16),16,32,16),[1,0,2])
		ca_tmp = reform(reform(ca_def,32l*16*16)#replicate(1.,64),32,16,16,64)

		ca_full = def4to16*ca_tmp*mas1to64

		if abs(total(cnts_ca) - total(ca_full)) gt 1. then begin
			print,'Error in ca normalization, times differ by:',time_ca-time_c6,' at ',time_string(time_c6)
			print,'Normalization changes counts by:',abs(total(cnts_ca)-total(ca_full)),total(ca_full),total(cnts_ca),total(cnts_c8),total(cnts_c6),total(cnts_c6)^.5
			print,'ca_32E   ',abs(total(cnts_ca)-total(ca_32E))
			print,'ca_tmp   ',abs(total(cnts_ca)-total(ca_tmp/256.))
			print,'def4to16 ',abs(total(cnts_ca)-total(def4to16*ca_tmp/64.))
			print,'mas1to64 ',abs(total(cnts_ca)-total(mas1to64*ca_tmp/4.))
			print,'tfdef,tfmas,tfe32=',tfdef,tfmas,tfe32
			n_ca=n_ca+1
			print,n_ca
		endif

		d1_full = ca_full*(total(cnts_c6)/total(cnts_ca))
		n_ca_4s = n_ca_4s + 1

	endif else begin

		d1_full = fltarr(32,16,16,64)							; null array
		print,'Error in d1 normalization - no valid data at: ',time_string(time_c6)

	endelse

endif else begin

	cnts_c6 = reform(mvn_c6_dat.data[ii,*,*])
	time_c6 = mtime_c6(ii)

	nearest = min(abs(mtime_c0-time_c6),ii_c0)
	cnts_c0 = reform(mvn_c0_dat.data[ii_c0,*,*])
	time_c0 = mtime_c0[ii_c0] & delta_t_c0 = mvn_c0_dat.delta_t[ii_c0]

	nearest = min(abs(mtime_d0-time_c6),ii_d0)
;	cnts_d0 = reform(mvn_d0_dat.data[ii_d0,*,*,*])
	time_d0 = mtime_d0[ii_d0] & delta_t_d0 = mvn_d0_dat.delta_t[ii_d0]

	nearest = min(abs(mtime_c8-time_c6),ii_c8)
	cnts_c8 = reform(mvn_c8_dat.data[ii_c8,*,*])
	time_c8 = mtime_c8[ii_c8] 
	if (time_c8-time_c6) gt 3. then cnts_c8 = (cnts_c8+reform(mvn_c8_dat.data[0>(ii_c8-1),*,*]))/2.
	if (time_c8-time_c6) lt -3 then cnts_c8 = (cnts_c8+reform(mvn_c8_dat.data[max_ind_c8<(ii_c8+1),*,*]))/2.

	nearest = min(abs(mtime_ca-time_c6),ii_ca)
	cnts_ca = reform(mvn_ca_dat.data[ii_ca,*,*])
	time_ca = mtime_ca[ii_ca] 
	if (time_ca-time_c6) gt 3. then cnts_ca = (cnts_ca+reform(mvn_ca_dat.data[0>(ii_ca-1),*,*]))/2.
	if (time_ca-time_c6) lt -3 then cnts_ca = (cnts_ca+reform(mvn_ca_dat.data[max_ind_ca<(ii_ca+1),*,*]))/2.

	if abs(time_ca - time_c6) lt 6. then begin		; use ca to construct d1_full 

;	     if abs(time_ca-time_c8) lt 2. and abs(total(cnts_ca)-total(cnts_c8)) lt .5*total(cnts_ca)^.5 then begin
;		norm_c8  = transpose(reform(replicate(1.,4)#reform(total(reform(transpose(cnts_c8),4,4,32),1),4*32),16,32))
;		def4to16 = reform(reform(cnts_c8/(norm_c8+.00000001),32*16)#replicate(1.,16*64),32,16,16,64)
tfdef=0
;	     endif else begin
		norm_c8  = transpose(reform(replicate(1.,4)#reform(total(reform(transpose(cnts_c8>1.),4,4,32),1),4*32),16,32))
		def4to16 = reform(reform((cnts_c8>1.)/(norm_c8),32*16)#replicate(1.,16*64),32,16,16,64)
tfdef=1
;	     endelse

	     if abs(time_ca-time_c6) lt 2. and abs(total(cnts_ca)-total(cnts_c6)) lt .5*total(cnts_ca)^.5 then begin
		norm_c6  = reform(total(cnts_c6,2)#replicate(1.,64),32,64)
		mas1to64 = transpose(reform(reform(cnts_c6/(norm_c6+.00000001),32*64)#replicate(1.,16*16),32,64,16,16),[0,2,3,1])
tfmas=0
	     endif else begin
		norm_c6  = reform(total((cnts_c6>.1),2)#replicate(1.,64),32,64)
		mas1to64 = transpose(reform(reform((cnts_c6>.1)/(norm_c6),32*64)#replicate(1.,16*16),32,64,16,16),[0,2,3,1])
tfmas=1
	     endelse

	     if abs(time_ca-time_c6) lt 2. and abs(total(cnts_ca)-total(cnts_c6)) lt .5*total(cnts_ca)^.5 then begin
		en16to32 = reform((total(cnts_c6,2)/(reform(replicate(1.,2)#total(reform(total(cnts_c6,2),2,16),1),32)+.0000001))#replicate(1.,64),32,64)
		ca_32E = en16to32*reform(replicate(1.,2)#reform(cnts_ca,16*64),32,64)
tfe32=0
	     endif else begin
		en16to32 = reform((total((cnts_c6>1.),2)/(reform(replicate(1.,2)#total(reform(total((cnts_c6>1.),2),2,16),1),32)))#replicate(1.,64),32,64)
		ca_32E = en16to32*reform(replicate(1.,2)#reform(cnts_ca,16*64),32,64)
tfe32=1
	     endelse

		ca_def = transpose(reform(replicate(1.,4)#reform(transpose(reform(ca_32E,32,4,16),[1,0,2]),4*32*16),16,32,16),[1,0,2])
		ca_tmp = reform(reform(ca_def,32l*16*16)#replicate(1.,64),32,16,16,64)

		ca_full = def4to16*ca_tmp*mas1to64

		if abs(total(cnts_ca) - total(ca_full)) gt 1. then begin
			print,'Error in ca normalization, times differ by:',time_ca-time_c6,' at ',time_string(time_c6)
			print,'Normalization changes counts by:',abs(total(cnts_ca)-total(ca_full)),total(ca_full),total(cnts_ca),total(cnts_c8),total(cnts_c6),total(cnts_c6)^.5
			print,'ca_32E   ',abs(total(cnts_ca)-total(ca_32E))
			print,'ca_tmp   ',abs(total(cnts_ca)-total(ca_tmp/256.))
			print,'def4to16 ',abs(total(cnts_ca)-total(def4to16*ca_tmp/64.))
			print,'mas1to64 ',abs(total(cnts_ca)-total(mas1to64*ca_tmp/4.))
			print,'tfdef,tfmas,tfe32=',tfdef,tfmas,tfe32
			n_ca=n_ca+1
			print,n_ca
		endif

		d1_full = ca_full*(total(cnts_c6)/total(cnts_ca))
		n_ca_4s = n_ca_4s + 1

	endif else begin

		d1_full = fltarr(32,16,16,64)							; null array
		print,'Error in d1 normalization - no valid data at: ',time_string(time_c6)

	endelse

endelse

	valid_full = reform((reform(valid,64*16)/(reform(replicate(1.,2)#reform(total(reform(valid,2,32,16),1),32*16),64*16)+.000001))#replicate(1.,16*64),64,16,16,64)
	d1_full_64E = valid_full*reform(replicate(1.,2)#reform(d1_full,32l*16*16*64),64,16,16,64)

;***************************************************************************************
;***************************************************************************************
; bkg2  - scattering off s/c - not working yet
;***************************************************************************************
;***************************************************************************************
; bkg3  - scattering off entrace grids/posts - not working yet
;***************************************************************************************
;***************************************************************************************
; bkg4 are internally backscattered cold ions near periapsis when the mechanical attenuator is closed
; 	assume all heavy ions outside of anodes 5 to 9, and with energy < Esh
;	Esh is the highest energy where the counts in c6 exceed 1/50. of the peak counts and exceed 20 counts.
;	the numbers 50 and 20 are arbitrary and should be tuned to make the algorithm work optimally 
;	tf_bkg4 is use to prevent double counting in background removal

   if not keyword_set(no_bkg4) then begin						; this just turns on/off this background

	tf_bkg4 = replicate(1,32,16,16,64) 
	if alt lt 500. and att ge 2 and mode ne 6 then begin
		c6_cnt = reform(total(mvn_c6_dat.data[ii,*,*],3))
		c6_nrg = reform(mvn_c6_dat.energy[swp_ind,*,0])
		pk_cnt = max(c6_cnt,pk_ind)
		pk_nrg = c6_nrg[pk_ind]
;		tf_ind = where(c6_cnt gt pk_cnt/50. and c6_cnt gt 10.,count)			; 50 and 10 are arbitrary, should be tuned
		tf_ind = where((c6_cnt gt pk_cnt/30.) or (c6_nrg lt 2.*pk_nrg),count)		; 30 and 2. are arbitrary, should be tuned
		if count gt 0 then begin
			tf_min = min(tf_ind)
			bkg4 = d1_full
			if tf_min ne 0 then bkg4[0:tf_min-1,*,*,*]=0.
			bkg4[tf_min:31,*,5:9,32:63]=0.
			bkg4[tf_min:31,*,*,0:31]=0.					; only remove heavies
			tf_bkg4[tf_min:31,*,0:4,32:63]=0
			tf_bkg4[tf_min:31,*,10:15,32:63]=0

; reorder for c8 -- most of the counts are in def_bin=11 for att=2, and def_bin=15 for att=3
		   if mode ne 1 and mode ne 7 then begin
			bkgtmp = total(bkg4[*,8:11,*,*],2)
			ratio1 = cnts_c8[*,11]/(total(total(bkgtmp,3),2) > cnts_c8[*,11] > .000001) 
			ratio2 = cnts_c8[*,10]/((1.-ratio1)*total(total(bkgtmp,3),2) > cnts_c8[*,10] > .000001) 
			trbkg4 = transpose(bkg4,[1,0,2,3])
			trbkg4[11,*,*,*] = bkgtmp * (ratio1#replicate(1.,16*64))
			trbkg4[10,*,*,*] = bkgtmp * ((ratio2*(1.-ratio1))#replicate(1.,16*64))
			trbkg4[9,*,*,*] = bkgtmp * (((1.-ratio2)*(1.-ratio1))#replicate(1.,16*64))
			trbkg4[8,*,*,*] = 0.
			bkg4 = transpose(trbkg4,[1,0,2,3])
		   endif else begin
			bkgtmp = total(bkg4[*,12:15,*,*],2)
			bkgtmp_tot = total(total(bkgtmp,3),2)
			ratio1 = cnts_c8[*,15]/(bkgtmp_tot > cnts_c8[*,15] > .000001) 
			ratio2 = cnts_c8[*,14]/((1.-ratio1)*bkgtmp_tot > cnts_c8[*,14] > .000001) 
			ratio3 = cnts_c8[*,13]/((1.-ratio2)*(1.-ratio1)*bkgtmp_tot > cnts_c8[*,13] > .000001) 
			trbkg4 = transpose(bkg4,[1,0,2,3])
			trbkg4[15,*,*,*] = bkgtmp * (ratio1#replicate(1.,16*64))
			trbkg4[14,*,*,*] = bkgtmp * ((ratio2*(1.-ratio1))#replicate(1.,16*64))
			trbkg4[13,*,*,*] = bkgtmp * ((ratio3*(1.-ratio2)*(1.-ratio1))#replicate(1.,16*64))
			trbkg4[12,*,*,*] = bkgtmp * (((1.-ratio3)*(1.-ratio2)*(1.-ratio1))#replicate(1.,16*64))
			bkgtmp2 = total(bkg4[*,8:11,*,*],2)
			trbkg4[11,*,*,*] = bkgtmp2 * .25			;assume it is constant
			trbkg4[10,*,*,*] = bkgtmp2 * .25
			trbkg4[9,*,*,*]  = bkgtmp2 * .25
			trbkg4[8,*,*,*]  = bkgtmp2 * .25
			bkg4 = transpose(trbkg4,[1,0,2,3])
		   endelse
		endif else bkg4=fltarr(32,16,16,64)
	endif else bkg4=fltarr(32,16,16,64)

   endif else begin
	bkg4=fltarr(32,16,16,64)
	tf_bkg4 = replicate(1,32,16,16,64) 
   endelse


;***************************************************************************************
;***************************************************************************************
; bkg5 are proton and alpha stragglers
; this section is turned off - the method was useful in understanding the source, but inadequate since anode variations dominated - will try another method

; bkg5a are proton stragglers - bkg5a ~rate^1 of valids - algorithm is still being developed 
; bkg5b are alpha  stragglers - bkg5b ~rate^1 of valids - algorithm is still being developed 
;	Might want to rewrite this to operate on tof_arr rather than mass_arr
;  	Note: bkg8 (~rate^2) looks like stragglers but are due to incomplete discharge of TOF capacitor from a Start-no-Stop event followed by a valid proton event

   if 0 and not keyword_set(no_bkg5) then begin

	offset_arr 	= transpose(reform(replicate(1.,32l*16*64)#mass_offset,32,16,64,16),[0,1,3,2])					; 64Ex16Dx16A
	mass = transpose(reform(reform(mvn_c6_dat.mass_arr[swp_ind,*,*],32l*64)#replicate(1.,16*16),32,64,16,16),[0,2,3,1]) - offset_arr

	twt_full = transpose(reform(reform(twt_arr,32*64)#replicate(1.,16*16),32,64,16,16),[0,2,3,1])

; TBD - the mass array should depend on rate due to TOF circuit -- the mass peaks shift to higher mass at high rates due to incomplete discharger of TOF capacitor 
;	mass = mass + rate 						

	max_anode = max(total(total(total(d1_full[*,*,*,0:7],4),2),1),ind_anode)

	p_cts = reform(reform(total(d1_full[*,*,*,0:7],4),32l*16*16)#replicate(1.,64),32,16,16,64)

;  Because different Start foil thicknesses affect straggling, the straggling function is assumed to have linear anode dependence
;	anode	     0     1     2     3     4     5     6     7     8     9    10    11    12    13    14    15  
	an_scale = [1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00]
	anode_scale = transpose(reform(replicate(1.,32l*16*64)#an_scale,32,16,64,16),[0,1,3,2])*p_cts/(transpose(reform(reform(total(p_cts,3),32l*16*64)#replicate(1.,16),32,16,64,16),[0,1,3,2])+.00000001)

;  this formula is imperically determined
;  aa0 is a discrete peak likely caused by grid inelastic scattering since it was not present in EM model without the field emission suppression grid
;  aa1,aa2 are carbon foil stragglers assumed to fit a pair of power law functions  
  
;	aa7 = 2.28e-5*.7
	aa7 = 1.60e-5
;	aa0 = (2.5 le mass and mass lt 3.5) *p_cts*0025*exp(-(mass-3.)^(2)/.3^2)*3.^(-2)		
	aa0 = (2.5 le mass and mass lt 3.5) *p_cts*2.78e-4*exp(-(mass-3.)^(2)/.3^2)		
;	aa1 = (1.5-mass_offset(ind_anode) le mass and mass lt 7.)  *p_cts*aa7*(mass/7.)^(-2.1)         	 
	aa1 = (1.5-offset_arr le mass and mass lt 7.)  *p_cts*aa7*(mass/7.)^(-2.1)         	 
	aa2 = (7. le mass and mass lt 100.) *p_cts*aa7*(mass/7.)^(-1.2)  

	bkg5a = anode_scale*(aa0+aa1+aa2)*twt_full	

;  bkg5b is due to alpha stragglers - TBD

; TBD - the mass array should depend on rate due to TOF circuit -- the mass peaks shift to higher mass at high rates due to incomplete discharger of TOF capacitor 
;	mass = mass + rate 						

	max_anode = max(total(total(total(d1_full[*,*,*,8:15],4),2),1),ind_anode)

	a_cts = reform(reform(total(d1_full[*,*,*,8:15],4),32l*16*16)#replicate(1.,64),32,16,16,64)

;  Because different Start foil thicknesses affect straggling, the straggling function is assumed to have linear anode dependence
;     these scale factors might be different for protons and alphas
;	anode	     0     1     2     3     4     5     6     7     8     9    10    11    12    13    14    15  
	an_scale = [1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00]
	anode_scale = transpose(reform(replicate(1.,32l*16*64)#an_scale,32,16,64,16),[0,1,3,2])*a_cts/(transpose(reform(reform(total(a_cts,3),32l*16*64)#replicate(1.,16),32,16,64,16),[0,1,3,2])+.00000001)

;  this formula is imperically determined
;  bb0 is a discrete peak likely caused by grid inelastic scattering since it was not present in EM model without the field emission suppression grid
;  bb1,bb2 are carbon foil stragglers assumed to fit a pair of power law functions  

; parameters are guessed - this needs tuning
	bb7 = 1.60e-5
	bb0 = (4.0 le mass and mass lt 5.6) *a_cts*2.78e-4*exp(-(mass-4.8)^(2)/.5^2)			
	bb1 = (2.5-offset_arr le mass and mass lt 8.)  *a_cts*bb7*(mass/8.)^(-2.1)         	 
	bb2 = (8. le mass and mass lt 100.) *a_cts*bb7*(mass/8.)^(-1.2)  

	bkg5b = anode_scale*(bb0+bb1+bb2)*twt_full	

	bkg5b = 0.

	bkg5_full = bkg5a + bkg5b	

    endif else bkg5_full=0.

;***************************************************************************************
; bkg5 are proton and alpha stragglers - bkg5a ~rate^1 of valids - algorithm is still being developed 
; this alternate bkg5 section is turned off - the method was useful in understanding the source, but inadequate since anode variations dominated - will try another method


	valid_full = reform((reform(valid,64*16)/(reform(replicate(1.,2)#reform(total(reform(valid,2,32,16),1),32*16),64*16)+.000001))#replicate(1.,16*64),64,16,16,64)
	d1_full_64E = valid_full*reform(replicate(1.,2)#reform(d1_full,32l*16*16*64),64,16,16,64)

   if 0 and not keyword_set(no_bkg5) then begin

	bkg5_full = fltarr(64,16,16,64)

		b_ns  = 5.844 
		tof_offset=21						; check
		int_time = 0.00380
		tau7 = 6.5   		; ns, this works if bkg6 is zero
		scale1 = 0.0030		; this works if bkg6 is zero

		tau7 = 4.5   		; ns, this works if for low rates if bkg6 is zero, but is 2-3 low at high rates
		scale1 = 0.0040		; this works if bkg6 is zero

		scale1 = 0.0000		; ignore this

		tau7 = 4.5   		; ns, this works if for low rates if bkg6 is zero, but is 2-3 low at high rates
		scale1 = 0.0025		; this works if bkg6 is zero

;		tau7 = 6.5   		; 
;		scale1 = 0.0020		; this works if bkg6 is zero
;		scale1 = 0.0015		; this works if bkg6 is zero

;		tau7 = 0.5   		; ns, is used if bkg6 handles low flux case
;		scale1 = 0.03
;		scale1 = 0.00

		tau7 = 4.5   		; 
		scale1 = 0.0030		; 

		tau8 = .65   		; ns
;		scale2 = 0.0029		; this works if bkg6 is zero
;		scale2 = 0.0024		;
		scale2 = 0.0020		; 20200619
		delta_tof = 8.6

	for m=0,63 do begin		; energy loop
		tof_tmp = reform(tof_arr[fix(m/2),*])			; tof_arr[32,64]
		twt_tmp = reform(twt_arr[fix(m/2),*])
		dtof1 = (1./b_ns)*replicate(1.,64)#(tof_offset+tof_tmp-twt_tmp/2.)
		dtof2 = (1./b_ns)*replicate(1.,64)#(tof_offset+tof_tmp+twt_tmp/2.)
		tof  = (1./b_ns)*(tof_offset+tof_tmp)#replicate(1.,64)
		tof22  = (1./b_ns)*(tof_offset+tof_tmp+twt_tmp/2.)#replicate(1.,64)
		tof33 = tof + delta_tof

		ttt1 = scale1*exp(((tof-dtof1)<0)/tau7) * (tof+1. lt dtof1)
		ttt2 = scale1*exp(((tof-dtof1)<0)/tau7) * (tof+1. lt dtof1)

		ttt3 = scale2*exp(-(tof33-dtof1)^2/tau8^2)
		ttt4 = scale2*exp(-(tof33-dtof2)^2/tau8^2)

		dttt1 = (dtof2-dtof1)*(ttt1+ttt2)/2.
		dttt2 = (dtof2-dtof1)*(ttt3+ttt4)/2.
; not needed
if 0 then begin
		for s=0,63 do begin				; this makes sure these events are not thrown away but shifted down in mass
			ind = where(dttt[s,*] gt 0.,count)
			if count ge 1 then dttt[s,ind[0]-1] = -1.*total(dttt[s,*])
;			ind = where(dttt[*,s] gt 0.,count)
;			if count ge 1 then dttt[ind[0]-1,s] = -1.*total(dttt[s,*])
		endfor	
endif
;  Because different Start foil thicknesses affect straggling, the straggling function may have a linear anode dependence
;	anode	     0     1     2     3     4     5     6     7     8     9    10    11    12    13    14    15  
		an_scale1 = [1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00]
		an_scale2 = [1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00]
		anode_scale1 = transpose(reform(replicate(1.,16*64)#an_scale1,16,64,16),[0,2,1])
		anode_scale2 = transpose(reform(replicate(1.,16*64)#an_scale2,16,64,16),[0,2,1])

		d1_tmp1 = reform(d1_full_64E[m,*,*,*],256,64)			; this calculates bkg6 for all masses
		d1_tmp1[*,16:63]=0.
;		d1_tmp1[*,8:63]=0.
		bkg5_tmp1 = reform(dttt1##d1_tmp1,16,16,64)*anode_scale1  
		bkg5_tmp2 = reform(dttt2##d1_tmp1,16,16,64)*anode_scale2 
		bkg5_full[m,*,*,*] = bkg5_tmp1 + bkg5_tmp2 
;		bkg5_full[m,*,*,*] = bkg5_tmp1 

	endfor					

    endif else bkg5_full[*]=0.

;********************************************************************************************
; bkg6 background is caused by start_no_stop followed closely by a valid event 
;	this section is turned off - the method was useful in understanding the source, but inadequate since anode variations dominated - will try another method
; 	valid events can be caused by other forms of background including coincidence events
;	bkg6 for unphysical high mass events (>60 amu) is about 10% of coincidence rates

	valid_full = reform((reform(valid,64*16)/(reform(replicate(1.,2)#reform(total(reform(valid,2,32,16),1),32*16),64*16)+.000001))#replicate(1.,16*64),64,16,16,64)
	d1_full_64E = valid_full*reform(replicate(1.,2)#reform(d1_full,32l*16*16*64),64,16,16,64)

;	d1_full_64E_tmp = reform(replicate(.5,2)#reform(d1_full,32l*16*16*64),64,16,16,64)
;	d1_full_64E_tmp = d1_full_64E
;	d1_full_64E_tmp[*,*,*,24:63] = 0.	; retrict bkg6 to only for low mass ions

   if 0 and not keyword_set(no_bkg6) then begin

	bkg6_full = fltarr(64,16,16,64)
		ratio_startnostop_valid = 1.17				; (1-.65*.75)/(.65*.75*.9)  assumes no stop droop
		b_ns  = 5.844 
		tof_offset=21						; check
		int_time = 0.00380

;		tofmax = 25 & tau = 320	& dtof=2.0 & exptof=1					; ns   ns   dead2=660ns or 620ns???
;		tofmax = 16 & tau = 200	& dtof=2.0 & exptof=1					; ns
;		tofmax = 12 & tau = 230	& dtof=2.0 & exptof=1					; right amplitude but tofmax needs to extend to higher mass
;		tofmax = 20 & tau = 150	& dtof=2.0 & exptof=1					; cutoff at 6.5amu, want 5.5
;		tofmax = 18 & tau = 150	& dtof=2.0 & exptof=1					; cutoff at 5.9 want 5.5
;		tofmax = 17 & tau = 200	& dtof=2.0 & exptof=1					; cutoff at 5.5 want 5.5
;		tofmax = 17 & tau = 250	& dtof=2.0 & exptof=1					; cutoff at 5.5 want 5.5
;		tofmax = 16 & tau = 300	& dtof=1.5 & exptof=1					; cutoff at 5.5 want 5.5
;		tofmax = 16 & tau = 300	& dtof=0.5 & exptof=1					; cutoff at 5.5 want 5.5
;		tofmax = 16 & tau = 300	& dtof=1.5 & exptof=2.					; cutoff at 5.5 want 5.5
;		tofmax = 20 & tau = 280	& dtof=1.5 & exptof=2.				; cutoff at 5.5 want 5.5
;		tofmax = 20 & tau = 250	& dtof=1.5 & exptof=2.				; cutoff at 5.5 want 5.5
;		tofmax = 20 & tau = 250	& dtof=1.5 & exptof=1.5				; cutoff at 5.5 want 5.5
;		tofmax = 18 & tau = 250	& dtof=1.5 & exptof=1.8				; cutoff at 5.5 want 5.5
;		tofmax = 20 & tau = 250	& dtof=1.5 & exptof=1.8				; cutoff at 5.5 want 5.5
		tofmax = 20 & tau = 250	& dtof=0.7 & exptof=1.8				; cutoff at 5.5 want 5.5

if 0 then begin
;  this ignores the accumulated charge at high rates for sequential stopnostart events
	for m=0,63 do begin		; energy loop
		tof_tmp = reform(tof_arr[fix(m/2),*])			; tof_arr[32,64]
		twt_tmp = reform(twt_arr[fix(m/2),*])
		dtof1 = (1./b_ns)*replicate(1.,64)#(tof_offset+tof_tmp-twt_tmp/2.)
		dtof2 = (1./b_ns)*replicate(1.,64)#(tof_offset+tof_tmp+twt_tmp/2.)
		tof  = (1./b_ns)*(tof_offset+tof_tmp)#replicate(1.,64)
		tof22  = (1./b_ns)*(tof_offset+tof_tmp+twt_tmp/2.)#replicate(1.,64)
;		ttt1 = tau*(0.>alog( tofmax/(((tof+tofmax)<dtof1>(tof+dtof))-tof ) ))^exptof 
;		ttt2 = tau*(0.>alog( tofmax/(((tof+tofmax)<dtof2>(tof+dtof))-tof ) ))^exptof 
		ttt1 = tau*(0.>alog( tofmax/(((tof+tofmax)<dtof1>(tof22))-tof ) ))^exptof 
		ttt2 = tau*(0.>alog( tofmax/(((tof+tofmax)<dtof2>(tof22))-tof ) ))^exptof 
		dttt = (ttt1-ttt2)>0.	
		for s=0,63 do begin				; this makes sure these events are not thrown away but shifted down in mass
			ind = where(dttt[s,*] gt 0.,count)
			if count ge 1 then dttt[s,ind[0]-1] = -1.*total(dttt[s,*])
;			ind = where(dttt[*,s] gt 0.,count)
;			if count ge 1 then dttt[ind[0]-1,s] = -1.*total(dttt[s,*])
		endfor	
		startnostop_rate = 1.e-9*ratio_startnostop_valid*reform(reform(valid[m,*])#replicate(1.,16*64),16,16,64)/int_time

   if 0 then begin
; 	this single iteration should work at low rates
;		d1_tmp = reform(startnostop_rate*d1_full_64E_tmp[m,*,*,*],16*16,64)		; this restricted bkg6 to only low mass
		d1_tmp = reform(startnostop_rate*d1_full_64E[m,*,*,*],16*16,64)			; this calculates bkg6 for all masses
		bkg6_full[m,*,*,*] = reform(dttt##d1_tmp,16,16,64)
   endif else begin
; 	this multi iteration is needed for high rates
;	mval = max(total(total(total(d1_full_64E,4),3),2),max_mind)
		d1_tmp1 = reform(startnostop_rate*d1_full_64E[m,*,*,*],16*16,64)			; this calculates bkg6 for all masses
		bkg6_tmp1 = reform(dttt##d1_tmp1,16,16,64)					
		d1_tmp2 = reform(startnostop_rate*(d1_full_64E[m,*,*,*]-bkg6_tmp1),16*16,64)			
		bkg6_tmp2 = reform(dttt##d1_tmp2,16,16,64)
		d1_tmp3 = reform(startnostop_rate*(d1_full_64E[m,*,*,*]-bkg6_tmp2),16*16,64)					
		bkg6_tmp3 = reform(dttt##d1_tmp3,16,16,64)
		d1_tmp4 = reform(startnostop_rate*(d1_full_64E[m,*,*,*]-bkg6_tmp3),16*16,64)					
		bkg6_tmp4 = reform(dttt##d1_tmp4,16,16,64)
		bkg6_full[m,*,*,*] = bkg6_tmp4

;	if m eq max_mind then print,total(abs(bkg6_tmp1)),total(abs(bkg6_tmp2)),total(abs(bkg6_tmp3)),total(abs(bkg6_tmp4)),$
;		total(abs(bkg6_tmp1-bkg6_tmp2)),total(abs(bkg6_tmp1-bkg6_tmp3)),total(abs(bkg6_tmp2-bkg6_tmp3)),total(abs(bkg6_tmp3-bkg6_tmp4))
   endelse
	endfor	
endif else begin
;  this attempts to account for the accumulated charge at high rates for sequential stopnostart events
		tofmax = 20 & tau = 250	& exptof=1.8				; cutoff at 5.5 want 5.5
		tofmax = 20 & tau = 250	& exptof=2.5				; cutoff at 5.5 want 5.5
		tofmax = 20 & tau = 250	& exptof=1.5				; cutoff at 5.5 want 5.5
		tofmax = 14 & tau = 250	& exptof=1.5				; cutoff at 5.5 want 5.5
		tofmax = 8 & tau = 250	& exptof=1.5				; cutoff at 5.5 want 5.5
		tofmax = 8 & tau = 500	& exptof=1.5				; cutoff at 5.5 want 5.5
		tofmax = 4 & tau = 500	& exptof=1.5				; cutoff at 5.5 want 5.5
		tofmax = 23 & tau = 500	& exptof=1.0				; cutoff at 5.5 want 5.5
		tofmax = 23 & tau = 1500	& exptof=1.0				; cutoff at 5.5 want 5.5
		tofmax = 10 & tau = 1500	& exptof=1.0				; cutoff at 5.5 want 5.5
		tofmax = 13 & tau = 1500	& exptof=1.0				; +/-3
		tofmax = 15 & tau = 1500	& exptof=1.0				; +/-5
		tofmax = 16 & tau = 1500	& exptof=1.0				; +/-6
		tofmax = 16 & tau = 1500	& exptof=1.0				; +/-6
	for m=0,63 do begin		; energy loop
;	mval = max(total(total(total(d1_full_64E,4),3),2),max_mind)
		tof_tmp = reform(tof_arr[fix(m/2),*])			; tof_arr[32,64]
		twt_tmp = reform(twt_arr[fix(m/2),*])
		dtof1 = (1./b_ns)*replicate(1.,64)#(tof_offset+tof_tmp-twt_tmp/2.)
		dtof2 = (1./b_ns)*replicate(1.,64)#(tof_offset+tof_tmp+twt_tmp/2.)
		tof  = (1./b_ns)*(tof_offset+tof_tmp)#replicate(1.,64)
		tof22  = (1./b_ns)*(tof_offset+tof_tmp+twt_tmp/2.)#replicate(1.,64)
	  for n=0,15 do begin
; 1.34, 1.94
		def_dead= 1.54
		def_dead= 1.64
		def_dead= 1.60
;		def_dead= 1.50
;		def_dead= 1.00
		tofmax0 = tofmax
if 0 then begin				; high coincidence rate
;		tofmax0 = 16
		tofmax1 = tofmax0-6.*(dead[m,n]/def_dead)
		tofmax2 = tofmax0
		tofmax3 = tofmax0+6.*(dead[m,n]/def_dead)
		tau  = 1500
;		tau2 = tau*(dead[m,n]/def_dead)^2.
		tau2 = tau*(dead[m,n]/def_dead)^2
endif else begin
		tofmax0 = 14.		; for low coincidence rates
		tofmax1 = tofmax0-3.	; was 3.
		tofmax2 = tofmax0
		tofmax3 = tofmax0+3.
		tau2 = 600

;		tofmax4 = tofmax0-4.	; was 3.
;		tofmax5 = tofmax0
;		tofmax6 = tofmax0+4.
;		def_dead= .60
;		tau3 = 1000*((dead[m,n]-1.)/def_dead)^2.

; this attempt involved using ratio_startnostop_valid to correct
		tofmax0 = 13.		; this works for high rates with no bkg5, but is 2x too big at low rates
		tofmax4 = tofmax0-3.	; was 3.
		tofmax5 = tofmax0
		tofmax6 = tofmax0+3.
		tau3 = 200.*rate_valid_ratio[m,n]	; rate_valid_ratio corrects for variations in ratio_startnostop_valid with rate

; this attempt uses ratio_start_rate
		tofmax0 = 13.		; 
		tofmax4 = tofmax0-3.	; was 3.
		tofmax5 = tofmax0
		tofmax6 = tofmax0+3.
		tau3 = 400.		; tbd

; this attempt uses ratio_start_rate
		tofmax0 = 15.		; 
		tofmax4 = tofmax0-5.	; was 3.
		tofmax5 = tofmax0
		tofmax6 = tofmax0+5.
;		tau3 = 200.		; tbd
		tau3 = 500.		; for bkg=5, works for high rate but factor of 3 low at low rate 
;		tau3 = 250.		; tbd

endelse
		ttt1a = tau2*(0.>alog( tofmax1/(((tof+tofmax1)<dtof1>(tof22))-tof ) ))^exptof 
		ttt1b = tau2*(0.>alog( tofmax1/(((tof+tofmax1)<dtof2>(tof22))-tof ) ))^exptof 
		ttt2a = tau2*(0.>alog( tofmax2/(((tof+tofmax2)<dtof1>(tof22))-tof ) ))^exptof 
		ttt2b = tau2*(0.>alog( tofmax2/(((tof+tofmax2)<dtof2>(tof22))-tof ) ))^exptof 
		ttt3a = tau2*(0.>alog( tofmax3/(((tof+tofmax3)<dtof1>(tof22))-tof ) ))^exptof 
		ttt3b = tau2*(0.>alog( tofmax3/(((tof+tofmax3)<dtof2>(tof22))-tof ) ))^exptof 
		ttt4a = tau3*(0.>alog( tofmax4/(((tof+tofmax4)<dtof1>(tof22))-tof ) ))^exptof 
		ttt4b = tau3*(0.>alog( tofmax4/(((tof+tofmax4)<dtof2>(tof22))-tof ) ))^exptof 
;		ttt5a = tau3*(0.>alog( tofmax5/(((tof+tofmax5)<dtof1>(tof22))-tof ) ))^exptof 
;		ttt5b = tau3*(0.>alog( tofmax5/(((tof+tofmax5)<dtof2>(tof22))-tof ) ))^exptof 
		ttt5a = 2.*tau3*(0.>alog( tofmax5/(((tof+tofmax5)<dtof1>(tof22))-tof ) ))^exptof 
		ttt5b = 2.*tau3*(0.>alog( tofmax5/(((tof+tofmax5)<dtof2>(tof22))-tof ) ))^exptof 
		ttt6a = tau3*(0.>alog( tofmax6/(((tof+tofmax6)<dtof1>(tof22))-tof ) ))^exptof 
		ttt6b = tau3*(0.>alog( tofmax6/(((tof+tofmax6)<dtof2>(tof22))-tof ) ))^exptof 
;		dttt = (ttt1-ttt2)>0.	
;		dttt = (((ttt1a-ttt1b)+(ttt2a-ttt2b)+(ttt3a-ttt3b))/3.)	
		dttt = (((ttt4a-ttt4b)+(ttt5a-ttt5b)+(ttt6a-ttt6b))/3.)	
;		dttt = (((ttt1a-ttt1b)+(ttt2a-ttt2b)+(ttt3a-ttt3b))/3.) + (((ttt4a-ttt4b)+(ttt5a-ttt5b)+(ttt6a-ttt6b))/3.)	
		for s=0,63 do begin				; this makes sure these events are not thrown away but shifted down in mass
			ind = where(dttt[s,*] gt 0.,count)
			if count ge 1 then dttt[s,ind[0]-1] = -1.*total(dttt[s,*])
;			ind = where(dttt[*,s] gt 0.,count)
;			if count ge 1 then dttt[ind[0]-1,s] = -1.*total(dttt[s,*])
		endfor	

;		startnostop_rate = 1.e-9*ratio_startnostop_valid*valid[m,n]*replicate(1.,16*64)/int_time
		startnostop_rate = 1.e-9*0.5*start_to_rate[m,n]*rate[m,n]*replicate(1.,16*64)

; 	this multi iteration is needed for high rates
		d1_tmp1 = reform(startnostop_rate*d1_full_64E[m,n,*,*],16,64)			; this calculates bkg6 for all masses
		bkg6_tmp1 = reform(dttt##d1_tmp1,16,64)					
		d1_tmp2 = reform(startnostop_rate*(d1_full_64E[m,n,*,*]-bkg6_tmp1),16,64)			
		bkg6_tmp2 = reform(dttt##d1_tmp2,16,64)
		d1_tmp3 = reform(startnostop_rate*(d1_full_64E[m,n,*,*]-bkg6_tmp2),16,64)					
		bkg6_tmp3 = reform(dttt##d1_tmp3,16,64)
		d1_tmp4 = reform(startnostop_rate*(d1_full_64E[m,n,*,*]-bkg6_tmp3),16,64)					
		bkg6_tmp4 = reform(dttt##d1_tmp4,16,64)
		bkg6_full[m,n,*,*] = bkg6_tmp4
	  endfor

;	if m eq max_mind then print,total(abs(bkg6_tmp1)),total(abs(bkg6_tmp2)),total(abs(bkg6_tmp3)),total(abs(bkg6_tmp4)),$
;		total(abs(bkg6_tmp1-bkg6_tmp2)),total(abs(bkg6_tmp1-bkg6_tmp3)),total(abs(bkg6_tmp2-bkg6_tmp3)),total(abs(bkg6_tmp3-bkg6_tmp4))

	endfor	

;tmp101 = total(total(bkg6_full>0,4),3)
;tmp102 = total(total(d1_full_64E,4),3)
;max_val = max(tmp101,ind_tmp)
;print,tmp101[ind_tmp],valid[ind_tmp],tmp102[ind_tmp]

endelse	
    endif else bkg6_full=0.

;***************************************************************************************
; bkg7 is background from coincident events where different ions (or radiation) produce start and stop
; bkg7 is proportional to total rate squared 
; bkg7 depends on the anode distribution of events - which differs for beams and penetrating radiation
; bkg7 is due to a "start_no_stop" followed by a stop from a second ion or penetrating background
; bkg7 will be underestimated when variations in solar wind flux occur during a sample period.  

; create arrays that are 64E x 16D x 16A x 8M

; TBD the following could be used to correct energy dependent efficiencies of high mass and molecular ions
	sp_mass_eff_e0 = [1.0,1.4,1.0,1.0,1.0,1.2,0.9,0.6]			; educated guess of mass stop eff at <100eV 
	sp_mass_eff_e1 = [1.0,1.2,1.0,1.0,1.0,1.2,0.9,0.6]			; educated guess of mass stop eff at 100eV-1keV 
	sp_mass_eff_e2 = [1.0,1.0,1.0,1.0,1.0,1.2,0.9,0.6]			; educated guess of mass stop eff at 1keV 
	sp_mass_eff_e3 = [1.0,1.0,1.0,1.0,1.0,1.4,1.7,1.7]			; educated guess of mass stop eff at >15keV 
;	sp_mass_eff = ????


   if not keyword_set(no_bkg7) then begin

	st_eff = replicate(tmp1,64,16,16)					; this could be anode dependent
	sp_eff = replicate(tmp2,64,16,16)					; this could be anode dependent
	qu_eff = replicate(tmp8,64,16)						; 64Ex16D
	qu_eff[*,*] = .8							; not sure which is better
	qu_eff[*,*] = tmp9>tmp8							; not sure which of these better
	qu_eff_full = reform(reform(qu_eff,64l*16)#replicate(1.,16),64,16,16)	; might want anode dependence in the future if it can be determined by inflight calib

;	64Ex16Dx16Ax64M 
	twt_full     = transpose(reform(reform(replicate(1.,2)#reform(twt_arr,32l*64),64*64)#replicate(1.,16*16),64,64,16,16),[0,2,3,1])

;	assume d1_full cnts reflects the Anode distribution of Rate
; 	anode_full: energy-def-anode distribution of counts, normalized at each energy-def
	d1_nomass    = total(d1_full,4)
	anode_tmp    = d1_nomass/(reform(reform(total(d1_nomass,3),32l*16)#replicate(1.,16),32,16,16)+.000000001)
	anode_full   = reform(replicate(1.,2)#reform(anode_tmp,32l*16*16),64,16,16)
	
; 	rate_full : energy-def-anode rate, assumes d1 anode distribution is the same as rate anode distribution
	rate_full    = reform(reform(rate,64*16)#replicate(1.,16),64,16,16)*anode_full

	dead_full    = reform(reform(dead,64*16)#replicate(1.,16),64,16,16)


; 	droop represents the decrease in start efficiency from the avg efficiency: st_eff=tmp1

;   droop_arr captures the non-rate dependent droop due to signal loss from delay line resistance - end anodes have more droop
;	droop_an_arr = replicate(1.,64*16)#droop_rate				; 
;	droop_arr    = total(droop_an_arr*anode_full,2)#replicate(1.,16)	; 
	droop_arr    = reform(replicate(1.,64*16)#droop_rate,64,16,16)		; st_no_sp has an anode dependent eff independent of droop

;	original bkg7 code for c6
;	avg_droop = replicate(total(valid*droop)/total(valid),64,16)
;	droop_eff_corr = avg_droop/droop

; 	assume anode dependence of droop is proportional anode rate; this line was replaced with the below 4 lines
;		d1_nomass64E    =               reform(replicate(.5,2)#reform(d1_nomass,32l*16*16),64,16,16)     

;   d1_nomass64E is used in calculated anode_norm
;   the >AAA prevents small values in droop_full, which would produce unphysically large droop_eff_corr values
;   AAA controls when the non-linear droop_eff_corr starts to be significant, may want something other than 1.
		en32to64 = total(cnts_c0,2)/(reform(replicate(1.,2)#total(reform(total(cnts_c0,2),2,32),1),64)+.00000001)
		en32to64_arr = reform(en32to64 # replicate(1.,16*16),64,16,16)
		d1_nomass64E    = en32to64_arr*(reform(replicate(1.,2)#reform(d1_nomass,32l*16*16),64,16,16)) >AAA	

;	print,total(en32to64),total(d1_nomass),total(d1_nomass64E),total(cnts_c0) 		; for testing

;  droop_eff_corr captures the non-linear part of droop where the highest flux anode will droop the most
;  this is imperically guessed with AAA controling the threshold for non-linear behavior

; this version does not wait the anodes properly
	anode_norm   = d1_nomass64E*reform(reform(total(d1_nomass64E,3)/(total(d1_nomass64E^2,3)+.000000001),64l*16)#replicate(1,16),64,16,16)

;	anode_norm   = (d1_nomass64E^2)/reform(reform(total(d1_nomass64E^2,3),64l*16)#replicate(1,16),64,16,16)

	droop_full   = reform(reform(droop,64l*16)#replicate(1.,16),64,16,16)*anode_norm
;	avg_droop    = total(d1_nomass64E*droop_full)/total(d1_nomass64E)
	avg_droop    = reform(reform(total(d1_nomass64E*droop_full,3)/total(d1_nomass64E,3),64*16)#replicate(1.,16),64,16,16)
	droop_eff_corr = avg_droop/droop_full
	droop_eff_corr[*]=1.
	droop_arr[*]=1.


	sp_eff_arr = reform(replicate(1.,64*16)#sp_an_eff,64,16,16)		; stop anode eff - stop anodes don't experience droop, but eff varies with anode
	ab_rate_ratio = replicate(tmp6/(tmp4+.001),64,16,16)			; A&B/RST
	st_no_sp_eff = (1.- sp_eff_arr)*ab_rate_ratio

	corr_rate = 1. + (1.-st_eff)*(1.-sp_eff)
	st_no_sp_eff = st_no_sp_eff*droop_eff_corr*droop_arr			; start eff is proportional to droop

	rbkg1 = st_no_sp_eff*rate_full*corr_rate*qu_eff_full*integ_t* (sp_eff_arr*rate_full*corr_rate*dead_full*1d*1.e-9/b_ns)

	rbkg2 = scale*reform(reform(rbkg1,64l*16*16)#replicate(1.,64),64,16,16,64)

	bkg7_full = rbkg2 * twt_full					; total bkg per tof bin

;		for testing
		if 0 then begin
		if (jj mod 4) eq 0 then begin
			print,minmax(anode_norm)
			print,minmax(droop_full)
		endif
		endif

   endif else bkg7_full=0.

;***************************************************************************************
;***************************************************************************************
; bkg8 are from Start-no-Stop followed closely by a TOF=0 event caused by Start-Stop cross-talk (these events are most easily see in apid DB)
; bkg8 is proportional to total rate squared, contaminates the low mass channels, and is generally caused by high rate O2+ at periapsis 

   if not keyword_set(no_bkg8) then begin
	bkg8_full = fltarr(64,16,16,64)
		ratio_startnostop_valid = 1.17				; (1-.65*.75)/(.65*.75*.9)
		crosstalk_to_valid_ratio = 1.7				; this is a measure of cross talk events
		b_ns  = 5.844 
		tof_offset=21						; check
		int_time = 0.00380
;		tofmax = 20 & tau = 250 & exptof=1.8			; these are the same as bkg6 except tof=0 and dtof=0
		tofmax = 25 & tau = 175 & exptof=1.8			; these should be tuned for anode 7 - the primary source of this background
		tofmax = 25 & tau = 120 & exptof=1.8			; tuned for anode 7 20191204 
	for m=0,63 do begin		; energy loop
		tof_tmp = reform(tof_arr[fix(m/2),*])
		twt_tmp = reform(twt_arr[fix(m/2),*])
		dtof1 = (1./b_ns)*replicate(1.,64)#(tof_offset+tof_tmp-twt_tmp/2.)
		dtof2 = (1./b_ns)*replicate(1.,64)#(tof_offset+tof_tmp+twt_tmp/2.)
		ttt1 = tau*(alog( tofmax/(tofmax<dtof1) ))^exptof 
		ttt2 = tau*(alog( tofmax/(tofmax<dtof2) ))^exptof
		dttt = (ttt1-ttt2)>0.		
		startnostop_rate = 1.e-9*ratio_startnostop_valid*reform(reform(valid[m,*])#replicate(1.,16*64),16,16,64)/int_time
		d1_tmp = reform(startnostop_rate*d1_full_64E(m,*,*,*),16*16,64)
		bkg8_full[m,*,*,*] = reform(dttt##d1_tmp,16,16,64)*crosstalk_to_valid_ratio
	endfor		
    endif else bkg8_full=0.


;***************************************************************************************
; bkg9 is background created by delayed Start signal from scattered molecular ion fragment that failed to produce a start at the carbon foil
; bkg9 is proportional to valid event rate - depicted in figures 33 & 34 of Mcfadden et al. 2015
; currently this background is only valid for low energy ions where the TOF bins are independent of energy

; bkg9a is primarily a low mass shoulder on the O2+ and CO2+ mass peaks seen in ground calibrations (N2+) and at Mars (O2+).
; 	This low mass shoulder background is only observed with molecules 
;	The low mass shoulder TOF relative to the primary molecule varies only with primary energy, not with TOF voltage.
;	This background is therefore internally produced in the TOF analyzer from a delayed START or early STOP signal, the latter having no known mechanism.
; 	These events are believed created by molecules hitting the start carbon foil grid (or suppression grid) and dissociating into two fragments
;	Neither fragment produces a normal START, while one fragment produces a normal STOP and the other a delayed START 
;	The delayed START is from a large angle scattering hitting an interal surface near the START foil and producing a secondary electron
;	Both O2+ and CO2+  exhibit a shoulder about 2 orders of magnitude below the nominal peak extending a factor of 2 in TOF below the peak

;	In addition, CO2+ exhibits a second peak about a factor of 17 below the primary peak at a TOF that peaks about a factor of 1.4 below the primary peak (see fig. 34 mcfadden et al.)
;	This secondary CO2+ peak is not observed in flight data at low energies and was likely a CO+ fragment from dissociation of keV ion molecules at the analyzer exit

; the following produces the low mass shoulder primarily from O2+, with contributions from other molecules such as CO2+, CO+, N2+, etc 

	bkg9_full = fltarr(32,16,16,64)

   if not keyword_set(no_bkg9) then begin

;	d1_full 	32E,16D,16A,64M
;	dead_full(E,D,A) - need to add mass
;	instead of cnts use d1_full
;	mass = reform(mvn_c6_dat.mass_arr[swp_ind,*,*]) - offset_arr

;	offset_an_arr 	= replicate(1.,64)#mass_offset				; 64 should be 32, but works anyway
;	offset_an	= reform(total(reform(anode,2,32,16),1),32,16)
;	offset_arr 	= total(offset_an_arr*offset_an,2)#replicate(1.,64)
;	mass = reform(mvn_c6_dat.mass_arr[swp_ind,*,*]) - offset_arr

	mass = reform(mvn_c6_dat.mass_arr[swp_ind,*,*])
	mass_full = transpose(reform(reform(mass,32*64)#replicate(1.,16*16),32,64,16,16),[0,2,3,1])
	twt_full = transpose(reform(reform(twt_arr,32*64)#replicate(1.,16*16),32,64,16,16),[0,2,3,1])
;	dead_full2 = reform(reform(dead_full,32l*16*16)#replicate(1.,64),32,16,16,64)

	dead_tmp = total(reform(dead*valid,2,32,16),1)/(total(reform(valid,2,32,16),1)+.000001)
	dead_full2    = reform(reform(dead_tmp,32*16)#replicate(1.,16*64),32,16,16,64)

	cxd_full1 = d1_full*dead_full2/twt_full*(mass_full gt 0 and mass_full lt 5.)				; primarily H+ - should have lower backscatter than heavies
	cxd_full2 = d1_full*dead_full2/twt_full*(mass_full gt 10 and mass_full le 24.)				; primarily O+ - should also contribute to bkg9b-9d
	cxd_full3 = d1_full*dead_full2/twt_full*(mass_full gt 24 and mass_full lt 100.)

;	bkg9a_full = .0020*twt_full*reform(tof2##reform(cxd_full3,32l*16*16,64),32,16,16,64)
;	bkg9a_full = .0025*twt_full*reform(tof2##reform(cxd_full3,32l*16*16,64),32,16,16,64)
	bkg9a_full = .0027*twt_full*reform(tof2##reform(cxd_full3,32l*16*16,64),32,16,16,64)			; 0.0027 emperically determined, only molecules contribute


	bkg9_full = bkg9a_full

   endif 

;***************************************************************************************
; bkg10 is TBD (disabled) bkg10 is internally produced background from sputtered ions generated by backscattered start-foil secondary electrons 
;	only heavy ions produce this peak - O+, O2+, CO2+ and O+ -- observed in ground calibrations
;	the primary sputtered ion peak is C++
;	however C+ is likely produced but it blends into the bkg9 shoulder so is not well estimated and can be considered part of bkg9
;	counts at C+++,C++++ are also observed but are much smaller and do not seem consistent with bkg8 background mechanism
;	a sputtered O+ might be generated, but is difficult to quantify since the bkg9 shoulder dominates
;	sputtered O++ is not generated, but actual ionospheric O++ is observed, or charge exchanged O+-->O++ during high O+ fluxes
;	the amplitude of the peaks depends on attentuator state
;	the TOF position of the peaks (especially C++) depends slightly on attuator state 
;	the size of this peak probably depends on the amount of disolved CO2 on surfaces and therefore may be anode and time dependent
;		since this peak is only important at periapsis, anode dependence can be ignored
;		since these peaks see a steady dose of CO2 at periapsis, they may not vary much in time or reflect dissolved surface gas prior to launch
;		time dependence should be checked after a deep dip 

; bkg10a is an ion peak at M/Q ~ 12, sputtered C+, similar to bkg10b 

; bkg10b is an ion peak at M/Q ~ 6-7, sputtered C++, which shifts with TOF voltage but not with primary molecule energy (see fig. 34a mcfadden et al.)
;	This peak is believed to be C++ caused by backscattered secondary electrons produced at the START carbon foil by primary molecules
;	These electrons are accelerated to 15keV and strike the ESA sputtering C++ from surfaces (probably residual CO2 gas on the surfaces).   
;	Current values for the center of these mass peaks is for anode 7 and the code may require tuning for other anodes

; bkg10c is an ion peak at M/Q ~ 6, assumed sputtered C+++,  similar to bkg10b 
; bkg10d is an ion peak at M/Q ~ 4, assumed sputtered C++++, similar to bkg10b 

; the following produces the low mass shoulder primarily from O2+, with contributions from other molecules such as CO2+, CO+, N2+, etc 

	bkg10_full = fltarr(32,16,16,64)

   if not keyword_set(no_bkg10) then begin

;	d1_full 	32E,16D,16A,64M
;	dead_full(E,D,A) - need to add mass
;	instead of cnts use d1_full
;	mass = reform(mvn_c6_dat.mass_arr[swp_ind,*,*]) - offset_arr

;	offset_an_arr 	= replicate(1.,64)#mass_offset				; 64 should be 32, but works anyway
;	offset_an	= reform(total(reform(anode,2,32,16),1),32,16)
;	offset_arr 	= total(offset_an_arr*offset_an,2)#replicate(1.,64)
;	mass = reform(mvn_c6_dat.mass_arr[swp_ind,*,*]) - offset_arr

	mass = reform(mvn_c6_dat.mass_arr[swp_ind,*,*])
	mass_full = transpose(reform(reform(mass,32*64)#replicate(1.,16*16),32,64,16,16),[0,2,3,1])
	twt_full = transpose(reform(reform(twt_arr,32*64)#replicate(1.,16*16),32,64,16,16),[0,2,3,1])
;	dead_full2 = reform(reform(dead_full,32l*16*16)#replicate(1.,64),32,16,16,64)

	dead_tmp = total(reform(dead*valid,2,32,16),1)/(total(reform(valid,2,32,16),1)+.000001)
	dead_full2    = reform(reform(dead_tmp,32*16)#replicate(1.,16*64),32,16,16,64)

	cxd_full1 = d1_full*dead_full2/twt_full*(mass_full gt 0 and mass_full lt 5.)				; primarily H+ - doesn't seem to generate this sputtered peak
	cxd_full2 = d1_full*dead_full2/twt_full*(mass_full gt 10 and mass_full le 24.)				; primarily O+ - should also contribute to bkg10a-10d
	cxd_full3 = d1_full*dead_full2/twt_full*(mass_full gt 24 and mass_full lt 100.)				; primarily O2+ and CO2+

; the following needs to be tuned for O+ and H+ backscattered secondaries - originally tuned for O2+ (cxd_full3)
; the constants may also change with time and anode, depending on ram direction at periapsis and therefore surface density of CO2 at ESA exit 

; attenuator variation in sputtered ion TOF center (assumed due to variations in where sputtered ions enter TOF)
	sputter_off = (att eq 3)*0.0 + (att eq 2)*0.0 + (att eq 1)*0.6 + (att eq 0)*0.6
; attenuator variation in sputtered ion TOF gain (assumed due to variations in where sputtered ions enter TOF)
	sputter_gain = 1. + (att eq 1)*0.3 + (att eq 0)*0.3
		
; sputtered C+ 
	bkg10a_full = 0
	bkg10a_full = sputter_gain*0.0010 *twt_full*reform((replicate(1.,64)#exp(-(tof1-43.5+sputter_off)^2/3.5^2))##reform(cxd_full3,32l*16*16,64),32,16,16,64)		; C+ sputtered from molecules
;	bkg10a_full = sputter_gain*0.0005 *twt_full*reform((replicate(1.,64)#exp(-(tof1-43.5+sputter_off)^2/3.5^2))##reform(cxd_full2,32l*16*16,64),32,16,16,64) + bkg10a_full	; C+ sputtered from O+
;	bkg10a_full = sputter_gain*0.00000*twt_full*reform((replicate(1.,64)#exp(-(tof1-43.5+sputter_off)^2/3.5^2))##reform(cxd_full1,32l*16*16,64),32,16,16,64) + bkg10a_full	; C+ sputtered from H+

; sputtered C++ 
; 0.0018  determined for 20160401
	bkg10b_full = sputter_gain*0.0018 *twt_full*reform((replicate(1.,64)#exp(-(tof1-31.7+sputter_off)^2/1.7^2))##reform(cxd_full3,32l*16*16,64),32,16,16,64)			; C++ sputtered from molecules
	bkg10b_full = sputter_gain*0.0009 *twt_full*reform((replicate(1.,64)#exp(-(tof1-31.5+sputter_off)^2/1.7^2))##reform(cxd_full2,32l*16*16,64),32,16,16,64) + bkg10b_full	; C++ sputtered from O+
;	bkg10b_full = sputter_gain*0.00000*twt_full*reform((replicate(1.,64)#exp(-(tof1-31.5+sputter_off)^2/1.7^2))##reform(cxd_full1,32l*16*16,64),32,16,16,64) + bkg10b_full	; C++ sputtered H+

; sputtered C+++   
; 0.0003  determined for 20160401
	bkg10c_full = 0
	bkg10c_full = sputter_gain*0.0003 *twt_full*reform((replicate(1.,64)#exp(-(tof1-25.0+sputter_off)^2/1.7^2))##reform(cxd_full3,32l*16*16,64),32,16,16,64)			; C+++ sputtered from molecules
	bkg10c_full = sputter_gain*0.00015*twt_full*reform((replicate(1.,64)#exp(-(tof1-23.0+sputter_off)^2/1.7^2))##reform(cxd_full2,32l*16*16,64),32,16,16,64) + bkg10c_full		; C+++ sputtered from O+
;	bkg10c_full = sputter_gain*0.00000*twt_full*reform((replicate(1.,64)#exp(-(tof1-23.0+sputter_off)^2/1.7^2))##reform(cxd_full1,32l*16*16,64),32,16,16,64) + bkg10c_full	; C+++ sputtered from H+

; sputtered C++++ 
; 0.0002  determined for 20160401
	bkg10d_full = 0
	bkg10d_full = sputter_gain*0.0002 *twt_full*reform((replicate(1.,64)#exp(-(tof1-21.0+sputter_off)^2/1.7^2))##reform(cxd_full3,32l*16*16,64),32,16,16,64)			; C++++ sputtered from molecules
	bkg10d_full = sputter_gain*0.0001 *twt_full*reform((replicate(1.,64)#exp(-(tof1-21.0+sputter_off)^2/1.7^2))##reform(cxd_full2,32l*16*16,64),32,16,16,64) + bkg10d_full		; C++++ sputtered from O+
;	bkg10d_full = sputter_gain*0.00000*twt_full*reform((replicate(1.,64)#exp(-(tof1-21.0+sputter_off)^2/1.7^2))##reform(cxd_full1,32l*16*16,64),32,16,16,64) + bkg10d_full		; C++++ sputtered from H+

	bkg10_full = bkg10a_full + bkg10b_full + bkg10c_full + bkg10d_full

   endif 

;***************************************************************************************
;***************************************************************************************
;***************************************************************************************
;***************************************************************************************

;   add background to all data except c0

	bkg_tmp = fltarr(32,16,16,64)
	if not keyword_set(no_bkg5) and total(bkg5_full) gt 0  then bkg_tmp=bkg_tmp + total(reform(bkg5_full,2,32,16,16,64),1)
	if not keyword_set(no_bkg6) and total(bkg6_full) gt 0  then bkg_tmp=bkg_tmp + total(reform(bkg6_full,2,32,16,16,64),1)
	if not keyword_set(no_bkg7) and total(bkg7_full) gt 0  then bkg_tmp=bkg_tmp + total(reform(bkg7_full,2,32,16,16,64),1)
	if not keyword_set(no_bkg8) and total(bkg8_full) gt 0  then bkg_tmp=bkg_tmp + total(reform(bkg8_full,2,32,16,16,64),1)
	if not keyword_set(no_bkg9) and total(bkg9_full) gt 0  then bkg_tmp=bkg_tmp + bkg9_full
	if not keyword_set(no_bkg10) and total(bkg10_full) gt 0  then bkg_tmp=bkg_tmp + bkg10_full

	bkg_full = bkg4 + tf_bkg4*bkg_tmp			; this assures no double counting of background

	if not keyword_set(no_update) then begin
		; note that if "add" keyword was NOT set, then bkg was zeroed out at beginning of program

		mvn_c6_dat.bkg[ii,*,*] =  mvn_c6_dat.bkg[ii,*,*] + total(total(bkg_full,3),2)

	     if size(mvn_d1_dat,/type) eq 8 then begin
		if time_c6 gt mvn_d1_dat.time[ii_d1] and time_c6 lt mvn_d1_dat.end_time[ii_d1] then $
			mvn_d1_dat.bkg[ii_d1,*,*,*] = mvn_d1_dat.bkg[ii_d1,*,*,*] + total(total(reform(bkg_full,32,4,4,16,8,8),5),2)
	     endif

		if time_c6 gt mvn_d0_dat.time[ii_d0] and time_c6 lt mvn_d0_dat.end_time[ii_d0] then $
			mvn_d0_dat.bkg[ii_d0,*,*,*] = mvn_d0_dat.bkg[ii_d0,*,*,*] + total(total(reform(bkg_full,32,4,4,16,8,8),5),2)

		if time_c6 gt time_c8-1. and time_c6 lt time_c8+1. then $
			mvn_c8_dat.bkg[ii_c8,*,*] = mvn_c8_dat.bkg[ii_c8,*,*] + total(total(bkg_full,4),3)

		if time_c6 gt time_ca-1. and time_c6 lt time_ca+1. then $
			mvn_ca_dat.bkg[ii_ca,*,*] = mvn_ca_dat.bkg[ii_ca,*,*] + reform(total(total(total(reform(bkg_full,2,16,4,4,16,64),6),3),1),16,64)

	endif

;***************************************************************************************
; treat c0 differently because it has 64 energies

   if not keyword_set(no_update) then begin
		; note that if "add" keyword was NOT set, then bkg was zeroed out at beginning of program

	bkg_tmp2 = fltarr(64,16,16,64)
	c0_cnt = reform(mvn_c0_dat.data[ii_c0,*,*])

;	bkg5_32to64 = (reform(c0_cnt[*,0])/(reform(replicate(1.,2)#total(reform(c0_cnt[*,0],2,32),1),64)+.000000001))#replicate(1.,2)
;	bkg9_32to64 = (reform(c0_cnt[*,1])/(reform(replicate(1.,2)#total(reform(c0_cnt[*,1],2,32),1),64)+.000000001))#replicate(1.,2)

;	bkg5_32to64 = reform((reform(c0_cnt[*,0])/(reform(replicate(1.,2)#total(reform(c0_cnt[*,0],2,32),1),64)+.000000001))#replicate(1.,16l*16*64),64,16,16,64)
	bkg5_32to64 = reform((reform(c0_cnt[*,0])/(reform(replicate(1.,2)#total(reform(c0_cnt[*,0],2,32),1),64)+.000000001))#replicate(1.,16l*16*64),64,16,16,64)
	bkg9_32to64 = reform((reform(c0_cnt[*,1])/(reform(replicate(1.,2)#total(reform(c0_cnt[*,1],2,32),1),64)+.000000001))#replicate(1.,16l*16*64),64,16,16,64)
	bkg10_32to64 = bkg9_32to64

;	if not keyword_set(no_bkg5) then bkg_tmp2 = bkg_tmp2 + bkg5_32to64*reform(replicate(1.,2)#reform(total(total(total(reform(tf_bkg4*bkg5_full,32,16,16,32,2),4),3),2),64),64,2)
;	if not keyword_set(no_bkg7) then bkg_tmp2 = bkg_tmp2 + reform(replicate(1,2)#reform(tf_bkg4,32l*16*16*64),64,16,16,64)*bkg7_full
;	if not keyword_set(no_bkg8) then bkg_tmp2 = bkg_tmp2 + bkg8_32to64*reform(replicate(1.,2)#reform(total(total(total(reform(tf_bkg4*bkg8_full,32,16,16,32,2),4),3),2),64),64,2)
;	if not keyword_set(no_bkg9) then bkg_tmp2 = bkg_tmp2 + bkg9_32to64*reform(replicate(1.,2)#reform(total(total(total(reform(tf_bkg4*bkg9_full,32,16,16,32,2),4),3),2),64),64,2)

;	if not keyword_set(no_bkg5) then bkg_tmp2 = bkg_tmp2 + bkg5_32to64 * reform(replicate(1.,2)#reform(tf_bkg4*bkg5_full,32l*16*16*64),64,16,16,64)

	if not keyword_set(no_bkg5) and total(bkg5_full) gt 0 then bkg_tmp2 = bkg_tmp2 + reform(replicate(1,2)#reform(tf_bkg4,32l*16*16*64),64,16,16,64)*bkg5_full
	if not keyword_set(no_bkg6) and total(bkg6_full) gt 0 then bkg_tmp2 = bkg_tmp2 + reform(replicate(1,2)#reform(tf_bkg4,32l*16*16*64),64,16,16,64)*bkg6_full
	if not keyword_set(no_bkg7) and total(bkg7_full) gt 0 then bkg_tmp2 = bkg_tmp2 + reform(replicate(1,2)#reform(tf_bkg4,32l*16*16*64),64,16,16,64)*bkg7_full
	if not keyword_set(no_bkg8) and total(bkg8_full) gt 0 then bkg_tmp2 = bkg_tmp2 + reform(replicate(1,2)#reform(tf_bkg4,32l*16*16*64),64,16,16,64)*bkg8_full
	if not keyword_set(no_bkg9) and total(bkg9_full) gt 0 then bkg_tmp2 = bkg_tmp2 + bkg9_32to64 * reform(replicate(1.,2)#reform(tf_bkg4*bkg9_full,32l*16*16*64),64,16,16,64)
	if not keyword_set(no_bkg10) and total(bkg10_full) gt 0 then bkg_tmp2 = bkg_tmp2 + bkg10_32to64 * reform(replicate(1.,2)#reform(tf_bkg4*bkg10_full,32l*16*16*64),64,16,16,64)

	bkg4_32to64 = reform((reform(c0_cnt[*,1])/(reform(replicate(1.,2)#total(reform(c0_cnt[*,1],2,32),1),64)+.000000001))#replicate(1.,16l*16*64),64,16,16,64)
	bkg4_64E = bkg4_32to64 * reform(replicate(1.,2)#reform(bkg4,32l*16*16*64),64,16,16,64)

	bkg_full2 = bkg4_64E + bkg_tmp2			; this assures no double counting of background

	if time_c6 gt time_c0-1. and time_c6 lt time_c0+1. then $
		mvn_c0_dat.bkg[ii_c0,*,*] = mvn_c0_dat.bkg[ii_c0,*,*] + total(total(total(reform(bkg_full2,64,16,16,32,2),4),3),2)

   endif

;*************************************************************************************************

; the following was for testing and is turned off
if 0 and keyword_set(save_d1_full) then begin

; the following arrays were too large to work with
;	d1_full_arr[jj,*,*,*,*] = d1_full
;	d1_full_bkg[jj,*,*,*,*] = bkg_full
;	d1_full_tim[jj] = (mvn_c6_dat.time[ii]+mvn_c6_dat.end_time[ii])/2.

	d1_full_ct1[jj,*,*,*] = total(d1_full,2)				; (32Ex16Ax64M)
	d1_full_bg1[jj,*,*,*] = total(bkg_full,2)				; (32Ex16Ax64M)
	d1_full_ct2[jj,*,*,*,*] = total(reform(d1_full,32,16,16,8,8),4)		; (32Ex16Dx16Ax8M)
	d1_full_bg2[jj,*,*,*,*] = total(reform(bkg_full,32,16,16,8,8),4)	; (32Ex16Dx16Ax8M)
	d1_full_tim[jj] = mvn_c6_dat.time[ii]
	d1_full_end[jj] = mvn_c6_dat.end_time[ii]
	d1_full_nrg[jj,*,*] = mvn_c6_dat.energy[swp_ind,*,*]
	d1_full_mas[jj,*,*] = mvn_c6_dat.mass_arr[swp_ind,*,*]
	d1_full_tof[jj,*,*] = tof_arr
	d1_full_twt[jj,*,*] = twt_arr
endif


endfor
; this is the end of the ii for loop

;*************************************************************************************************

; diagnostics for testing bkg6 - turned off

if 0 then begin
	store_data,'start_to_rate_min',data={x:mtime_c6[ind_jj],y:st_to_rate_min}
		options,'start_to_rate_min',ylog=1,yrange=[.1,1]
	store_data,'rate_max',data={x:mtime_c6[ind_jj],y:rate_max}
		options,'rate_max',ylog=1,yrange=[1.e5,1.e6]
	store_data,'st_to_rate_diff',data={x:mtime_c6[ind_jj],y:st_to_rate_diff}
		options,'start_to_rate_diff',ylog=1,yrange=[.01,1]
	store_data,'rate0_arr',data={x:mtime_c6[ind_jj],y:rate0_arr}
		options,'rate0_arr',ylog=1,yrange=[1.e4,1.e6]
 	store_data,'n_iter_arr',data={x:mtime_c6[ind_jj],y:n_iter_arr}
		options,'n_iter_arr',ylog=0,yrange=[0,16]
 
	if keyword_set(save_d1_full) then begin

; 	the following array was too large to work with
;		d1_full_str = {time:d1_full_tim,energy:d1_full_nrg,cnts:d1_full_arr,bkg:d1_full_bkg,mass:d1_full_mas,tof:d1_full_tof,twt:d1_full_twt}

		d1_full_str = {time:d1_full_tim,end_time:d1_full_end,energy:d1_full_nrg,ct1:d1_full_ct1,bg1:d1_full_bg1,$
			ct2:d1_full_ct2,bg2:d1_full_bg2,mass:d1_full_mas,tof:d1_full_tof,twt:d1_full_twt}
		save,d1_full_str,filename='C:\data\maven\idlsave\d1_full_str\'+yrmodatm+'_idlsave_d1_full_str'
	endif

endif

print,'bkg7  n_d1_4s,n_d1_16s,n_ca_4s=',n_d1_4s,n_d1_16s,n_ca_4s
print,'Run time mvn_sta_bkg_load:',systime(1)-starttime
	
end