function mvn_sep_sim_rotate_by_180,data1          ; simulate the B side
  data2 = data1
  data2.pos[0] = -data1.pos[0]
  data2.pos[1] = -data1.pos[1]
  data2.dir[0] = -data1.dir[0]
  data2.dir[1] = -data1.dir[1]
  data2.edep[*,0] = data1.edep[*,1]
  data2.edep[*,1] = data1.edep[*,0]
  ;data2.a      = data1.b
  ;data2.b      = data1.a
  return,data2
end


function mvn_sep_inst_response_retrieve,pathnames,age_limit=age_limit
  pdunn = file_retrieve(/struct)
  if ~keyword_set(age_limit) then age_limit=900 ; 3600*4
  pdunn.min_age_limit=age_limit
  pdunn.local_data_dir = '~/data/pdunn/'
  pdunn.remote_data_dir = 'http://sprg.ssl.berkeley.edu/~davin/pdunn/'
  pdunn.archive_ext = '.arc'
  pdunn.archive_dir = 'archive/'
  pdunn.ignore_filesize = 1
  files = file_retrieve(pathnames,_extra=pdunn)
  return,files
end



pro mvn_sep_read_mult_sim_files,simstat,data,desc=desc,pathnames=pathnames,type=type,dosymm=dosymm
  last_pathname=''
  str_element,simstat,'pathnames',last_pathname
  if array_equal(pathnames,last_pathname) then begin
    dprint,dlevel=1, 'Pathnames: ',pathnames[0],' Match. Returning.'
    return
  endif
  simstat=0
  files = mvn_sep_inst_response_retrieve(pathnames)
  str_element,/add,simstat,'desc',desc
  str_element,/add,simstat,'files',files
  str_element,/add,simstat,'pathnames',pathnames
  dprint,'Loading: ',files
  fmt = {event:0L,einc:0.,pos:[0.,0.,0.],dir:[0.,0.,0.],edep:[[0.,0.,0.],[0.,0.,0.]],E_tot:0.}   ; [F T O] order in these files
  ;#    Event  |  (keV)  |  Pos_x     Pos_y     Pos_z    |  Dir_x     Dir_y     Dir_z    |     AF        AT        AO        BF        BT        BO     Total
  npart=0
  data=0
  ;e_inc_num=0.
  ;if keyword_set(xbinsize) then begin
  ;endif

  for i=0,n_elements(files)-1 do begin
    file = files[i]
    d = read_asc(file,simstat,format=fmt)
    ns = max(d.event)     ; estimate of number of test particles
    descfile = str_sub(file,'.dat','.txt')
    if file_test(/regular,descfile) then printdat,read_asc(descfile,simstat)
    ; str_element,simstat,'numberofparticles',ns
    dprint,dlevel=2,'Loaded',n_elements(d),' data samples out of ',ns,' in ',file
    npart += ns
    append_array,data,d
    ;  if keyword_set(xbinsize) then begin
    ;      e_inc_num +=
    ;  endif
  endfor

  data.edep = reverse(data.edep,1)   ; switch order to O T F
  data.dir = shift(data.dir,1,0)     ; rotate through axes
  data.pos = shift(data.pos,1,0)     ;

  if n_elements(type) eq 1 then str_element,/add,data,'type',replicate(type,n_elements(data))
  ;str_element,/add,data,'det_ed',reform([data.a,data.b],3,2,n_elements(data))
  tp = n_elements(type) eq 1 ? type : 0
  ;if ~keyword_set(npart) then npart = max(data.event)
  sim_area = 0.
  str_element,simstat,'sim_area',sim_area
  if 0 then begin
    w = where(data.event gt shift(data.event,-1),nroll)              ; account for multiple files
    npart = total(/preserve,data[w].event)
  endif else nroll =1

  str_element,/add,simstat,'npart',npart
  str_element,/add,simstat,'sim_area',sim_area * nroll     ; account for multiple files - assume all have same area.
  if keyword_set(dosymm) then begin
    data=[data,mvn_sep_sim_rotate_by_180(data)]
    simstat.sim_area *= 2
    simstat.npart *=2
  endif

  ;str_element,/add,simstat,'data',data
end


function reform2,vec
  if n_elements(vec) eq 1 then return, vec[0]
  return, reform(vec)
end



;+
;FUNCTION:  crossp_trans(a,b)
;INPUT:
; a,b:  real(3,n) vector arrays dimension (3,n) or (3)
;PURPOSE:
; performs cross product on arrays
;CREATED BY:
; Davin Larson
;-
function crossp_trans,a,b
  if n_params() ne 2 then message, ' Wrong format, Use: crossp2(a,b)'
  dim_a = size(/dimen,a)
  dim_b = size(/dimen,b)
  dim = n_elements(dim_a) eq 2 ? dim_a : dim_b   ; crude but works
  c=replicate(a[0]*b[0] , dim)
  ;   printdat,dim_a,dim_b,c
  c[0,*]= reform2(a[1,*]) * reform2(b[2,*]) - reform2(a[2,*]) * reform2(b[1,*])
  c[1,*]= reform2(a[2,*]) * reform2(b[0,*]) - reform2(a[0,*]) * reform2(b[2,*])
  c[2,*]= reform2(a[0,*]) * reform2(b[1,*]) - reform2(a[1,*]) * reform2(b[0,*])
  return,c
end



pro mvn_sep_response_simdata_rand,type,data=data,simstat=simstat,seed=seed,window=win
  if not keyword_set(type) then type = 0
  if keyword_set(simstat) then return
  dprint,'Creating simulated data distribution'
  d={event:0L,einc:0.,pos:[0.,0.,0.],dir:[0.,0.,0.],A:[0.,0.,0.],B:[0.,0.,0.],e_tot:0., type:0}
  d = {event:0L,einc:0.,pos:[0.,0.,0.],dir:[0.,0.,0.],edep:[[0.,0.,0.],[0.,0.,0.]],E_tot:0.}   ; [F T O] order in these files
  n = 20000000L
  n = 10000000L
  ;n = 1000000L
  log10_sr = [0,5d]
  str_element,/add,simstat,'sim_energy_range', 10d ^ log10_sr
  str_element,/add,simstat,'sim_energy_log', 1
  str_element,/add,simstat,'sphere_radius', 80.  ; mm
  str_element,/add,simstat,'sim_area', simstat.sphere_radius^2 *!dpi  ; mm^2
  str_element,/add,simstat,'npart',n
  str_element,/add,simstat,'desc','Sim'
  str_element,/add,simstat,'Particle_name','simpart'
  str_element,/add,simstat,'xbinsize',.1d
  ;str_element,/add,simstat,'area', 1.
  data = replicate(d,n)
  data.event = lindgen(n)
  data.einc = 10d^(log10_sr[0] + randomu(seed,n) * (log10_sr[1]-log10_sr[0]))
  ; Random distribution of starting points over surface of sphere
  phi = randomu(seed,n) * 360d
  theta = asin(randomu(seed,n) * 2-1) *180/!dpi
  sphere_to_cart,simstat.sphere_radius,theta,phi,vec=vec
  data.pos = transpose(vec)
  ; Get Random directions
  phi = randomu(seed,n) * 360d
  oldway=0
  if oldway then begin
    theta = asin(randomu(seed,n) * 2-1) *180/!dpi     ;  uniform distribution of directons
  endif else begin
    theta = asin( sqrt( randomu(seed,n)  )   ) *180/!dpi  ; non-uniform - Better for surface of sphere
  endelse
  sphere_to_cart,1,theta,phi,vec=vec
  data.dir = transpose(vec)
  ; Must rotate direction distribution into position of starting point on sphere
  if oldway then begin
    w = where( total(data.pos * data.dir,1) gt 0, nw)
    if nw ne 0 then data[w].dir = -data[w].dir              ; force inward trajectories  only
  endif else begin
    q = get_quaternion(-data.pos,/last_index)
    q[0,*] = -q[0,*]    ; inverse rotation
    temp = quaternion_rotation(data.dir,q,/last_index)
    data.dir = temp
  endelse
  if n_elements(win) ne 0 then begin
    wi,win,/show,wsize=[600,800]
    xyz_to_polar,vec,theta=th2,phi=ph2,mag=r,/ph_0_360
    !p.multi=[0,0,4]
    printdat,phi,theta
    plot,xb,histbins(phi,xb),xtitle='PHI'
    if keyword_set(ph2) then oplot,color=2,xb,histbins(phi,xb)
    plot,xb,histbins(theta,xb,/shift),xtitle='THETA'
    if keyword_set(th2) then oplot,color=6,xb,histbins(theta,xb)
    ;  plot,xb,histbins(r,xb,range=[0,2])
    psym=10
    plot,xb,histbins(vec[*,0],xb,/shift),yrange=n*[0,2.]/n_elements(xb),psym=psym
    oplot,xb,histbins(vec[*,1],xb,/shift),color=2,psym=psym
    oplot,xb,histbins(vec[*,2],xb,/shift),color=6,psym=psym

    plot,xb,histbins( total( data.pos *data.dir,1), xb),xtitle='POS dot DIR'
    !p.multi=0
    printdat,xb
  endif
  ;data.type = type
  emin = ([10.,0.,200.])[type+1]
  emin = 0
  data.edep[2,*] =  reform([1,1] # (data.einc-emin) > 0.,[1,2,n])
end


function mvn_sep_adc_calibration
  ;message,'obsolete.  Contained within mvn_sep_lut2map.pro'
  adc_scale =  [[[ 43.77, 38.49, 41.13 ] ,  $  ;1A          O T F
    [ 41.97, 40.29, 42.28 ]] ,  $  ;1B
    [[ 40.25, 44.08, 43.90 ] ,  $  ;2A
    [ 43.22, 43.97, 41.96 ]]]   ;  2B
  adc_scale = adc_scale / 59.5
  return,adc_scale
end

;
;function mvn_sep_response_flux_to_bin_cr,b,omega,pnum,response=r,parameter=par  ;,omega=omega
;if n_elements(omega) eq 0 then omega = [0,1]
;;n_omega = n_elements(omega)
;;n_einc  = n_elements(r.e_inc)
;;n_bins
;;types = [0,1]
;flux =  func(r.e_inc,omega,param=par)
;flux_dim = size(/dimen,flux)
;resp = r.bin3
;resp_dim = size(/dimen,resp)
;
;de_inc = (r.xlog) ? (r.e_inc * r.xbinsize  * alog(10)) : r.xbinsize
;de_inc = de_inc # replicate(1,flux_dim[1])
;de_flux = de_inc * flux
;
;;printdat,de_inc,de_flux
;resp = reform(resp,resp_dim[0]*resp_dim[1],resp_dim[2])
;de_flux = reform(de_flux,flux_dim[0] * flux_dim[1])
;;printdat,resp,de_flux
;CR = resp ## de_flux
;CR *=  (r.area_cm2 / r.nd)
;if keyword_set(b) then CR=CR[b]
;return,reform(CR)
;end
;


;function particle_flux,energy,omega,parameter=par,response=r,dflux50=dflux50,name=name
;if n_elements(omega) eq 0 then omega=[0,1]
;if ~keyword_set(par) || keyword_set(r) then begin
;    if ~keyword_set(r) then r=0
;    if ~keyword_set(name) then name = ''
;    if ~keyword_set(dflux50) then dflux50 = 100.
;    nrg = 10.^[1.5,2,2.5,3,3.5,4,4.5,6.]
;    flx = dflux50 * (nrg/50.)^(-3)
;    flux = spline_fit3(nrg,nrg,flx,/xlog,/ylog,par =pflux)
;    par = {func:'particle_flux', $
;           name:name, $
;           spec:replicate(pflux,n_elements(omega)) , $
;           response:r,  $
;           units_name:'flux'}
;    par.spec[1].ys -= .1
;    if n_params() eq 0 then return,par
;    printdat,par
;endif
;fluxes=fltarr(n_elements(energy),n_elements(omega))
;for i=0,n_elements(omega)-1 do begin
;  fluxes[*,i] = spline_fit3(energy,param=par.spec[omega[i]])
;endfor
;return,fluxes
;end



;
;function all_flux,energy,omega,parameter=par,Electron_response=re,Proton_response=rp,species=species
;if ~keyword_set(species) then species=''
;if n_elements(omega) eq 0 then omega=[0,1]
;if ~keyword_set(par) || keyword_set(re) || keyword_set(rp) then begin
;    electron = particle_flux(response=re,name='Electron',dflux50=.5)
;    proton   = particle_flux(response=rp,name='Proton',dflux50=10.)
;    par ={func:'all_flux', pts:{electron:electron, Proton:proton}, units:'CR' }
;    if n_params() eq 0 then return,par
;endif
;dt = size(/type,energy)
;if (dt eq 4) || (dt eq 5) then begin
;   el_flux = func(energy,omega,param=par.pts.electron)
;   pr_flux = func(energy,omega,param=par.pts.proton)
;   if species eq 'electron' then return,el_flux
;   if species eq 'proton' then return,pr_flux
;   return,el_flux + pr_flux
;endif
;;if n_params() eq 0 then energy=par.pts.electron.response.e_inc
;
;if 1 || par.units eq 'CR' then begin
;   cr_e = mvn_sep_response_flux_to_bin_CR(param=par.pts.electron,response=par.pts.electron.response)
;   if species eq 'electron' then return,cr_e
;   cr_p = mvn_sep_response_flux_to_bin_CR(param=par.pts.proton  ,response=par.pts.proton.response)
;   if species eq 'proton' then return,cr_p
;   cr_t = cr_p+ cr_e
;endif
;
;if dt eq 2  or dt eq 1 or dt eq 3 then return, cr_t[energy]    ; energy is the index
;
;return,cr_t
;
;end




;
;
;function particle_flux,energy,omega,type,nrg=nrg,flx=flx,parameter=par
;if keyword_set(nrg) and keyword_set(flx) then begin
;    newflux = spline_fit3(energy,nrg,flx,/xlog,/ylog,par=electron)
;    newflux = spline_fit3(energy,nrg,flx,/xlog,/ylog,par=ion)
;    par = {func:'particle_flux',$
;           electron:electron,       $
;           ion: ion, $
;           units_name:'Flux'}
;endif
;elec_flux = spline_fit3(energy,param=par.electron)
;ion_flux = spline_fit3(energy,param=par.ion)
;return,elec_flux
;end
;

;  Multiple matrix plots
pro mvn_sep_response_matrix_plots,r,window=win,single=single
  if keyword_set(win) then     wi,win,/show,wsize=[1100,850]
  labels = strsplit('XXX O T OT F FO FT FTO Total',/extract)
  zrange = minmax(r.g4,/pos)
  xrange = r.xbinrange
  yrange = r.ybinrange
  options,lim,xlog=1,/ylog,xrange=xrange,/ystyle,/xstyle,yrange=yrange,xmargin=[10,10],/zlog,zrange=zrange,/no_interp,xtitle='Energy incident (keV)',ytitle='Enery Deposited (keV)'
  if not keyword_set(single) then !p.multi = [0,4,4]
  if not keyword_set(ok1) then ok1 = 1
  wnum=0
  atten = (['Open','Closed'])[r.attenuator]
  SEP  = (['???','SEP1','SEP2'])[r.sensornum]
  title = r.desc+' '+r.particle_name+' '+SEP+' '+atten+' '

  for side =0,1 do begin
    slabel = side ? 'B' : 'A'
    for fto = 1,8 do begin
      if fto eq 8 then G2 = total(r.g4[*,*,1:7,side],3) $
      else G2 = r.g4[*,*,fto,side]
      dprint,dlevel=3,side,fto,total(g2)
      options,lim,title = title+slabel+'_'+labels[fto]
      if keyword_set(single) then begin
        if single ne fto then continue
        if side ne 0 then continue
      endif
      specplot,r.e_inc,r.e_meas,G2,limit=lim
      oplot,dgen(),dgen(),linestyle=1
    endfor
  endfor
  !p.multi=0
end


;  Multiple matrix plots
pro mvn_sep_response_matrix_plots_lin,r,window=win,single=single
  if keyword_set(win) then     wi,win,/show,wsize=[1100,850]
  labels = strsplit('XXX O T OT F FO FT FTO Total',/extract)
  zrange = minmax(r.g4,/pos)
  xrange = r.xbinrange
  yrange = r.ybinrange
  xrange = [0.,200]
  yrange = [0.,200]
  options,lim,xlog=0,ylog=0,xrange=xrange,/ystyle,/xstyle,yrange=yrange,xmargin=[10,10],/zlog,zrange=zrange,/no_interp,xtitle='Energy incident (keV)',ytitle='Enery Deposited (keV)'
  if not keyword_set(single) then !p.multi = [0,4,4]
  if not keyword_set(ok1) then ok1 = 1
  wnum=0
  for side =0,1 do begin
    slabel = side ? 'B' : 'A'
    for fto = 1,8 do begin
      if fto eq 8 then G2 = total(r.g4[*,*,1:7,side],3) $
      else G2 = r.g4[*,*,fto,side]
      dprint,dlevel=3,side,fto,total(g2)
      options,lim,title = slabel+'_'+labels[fto]
      if keyword_set(single) then begin
        if single ne fto then continue
        if side ne 0 then continue
      endif
      specplot,r.e_inc,r.e_meas,G2,limit=lim
      oplot,dgen(),dgen(),linestyle=1
    endfor
  endfor
  !p.multi=0
end



pro mvn_sep_response_plot_gf,r,window=win,ylog=ylog,xrange=xrange  ;,face=face
  ;            x O  T OT  F  OF  FT FTO   Total
  colors = [0,2, 4, 1, 6,  5,  3, 0,   5]
  linestyle = [0,2]
  yrange = [0,2]
  if ~keyword_set(face) then face=0
  face_str = (['Aft','Both','Front'])[face+1]
  atten = (['Open','Closed'])[r.attenuator]
  SEP  = (['???','SEP1','SEP2'])[r.sensornum]
  title = r.desc+' '+r.particle_name+' '+SEP+' '+atten+' '+face_str


  if keyword_set(win) then     wi,win,/show,wsize=[600,600]
  if keyword_set(ylog) then yrange=[.0001,100]
  str_element,r,'fdesc',subtitle
  einc = r.e_inc
  ;   str_element,r,'xbinrange',xrange
  if not keyword_set(xrange) then xrange = minmax(einc)
  for side=0,1 do begin
    ls = linestyle[side]
    G = total( total(r.G4[*,*,*,side],3), 2)
    if ~keyword_set(fst) then  plot,[1,2],/nodata,xrange=xrange,/xlog,/xstyle,thick=2,yrange=yrange,ylog=ylog,xtitle='Incident Energy (keV)',ytitle='Geometric Area (cm^2)',title=title,subtitle=subtitle
    fst=1
    ;     oplot , einc,G   ,color = colors[8],linestyle =ls   ; total
    for fto = 1,7 do begin
      G = total(r.g4[*,*,fto,side],2)
      oplot,einc,G > yrange[0]/3,color=colors[fto],linestyle=ls,thick=1
    endfor
  endfor
end

pro mvn_sep_response_plot,r,window=win,ylog=ylog,fluxfunc=fluxfunc,fst=fst
  colors = [0,2,4,1,6,5,3,0,5]
  linestyle = [0,2]
  yrange = [0,2]
  if keyword_set(win) then     wi,win,/show,wsize=[600,800]
  if keyword_set(ylog) then yrange=[.0001,1000]
  einc = r.e_inc
  emeas = r.e_meas
  flux = keyword_set(fluxfunc) ? func(param=fluxfunc,einc) : replicate(1,n_elements(einc))
  de_inc = (r.xlog) ? (einc * r.xbinsize  * alog(10)) : r.xbinsize
  deflux = de_inc * flux
  de_meas = (r.ylog) ? (emeas * r.ybinsize *alog(10)) : r.ybinsize
  dprint,dlevel=1,r.desc,r.xbinsize,r.ybinsize
  for side=0,1 do begin
    slabel = side ? 'B' : 'A'
    ls = linestyle[side]
    CR= reform( total(r.G4[*,*,1:7,side],3) ## deflux ) / de_meas
    dprint,'Total Counts in ',slabel,':',total(cr * de_meas)
    if ~keyword_set(fst) then  plot,[1,2],/nodata,xrange=minmax([emeas,einc]),/xlog,yrange=yrange,ylog=ylog,xtitle='Energy (keV)',ytitle='Count Rate per keV',title=r.desc
    fst=1
    oplot , emeas,CR   ,color = colors[8],linestyle =ls,psym=10
    for fto = 1,7 do begin
      CR = reform(r.g4[*,*,fto,side] ## deflux) /de_meas
      oplot,emeas,CR,color=colors[fto],linestyle=ls,psym=10
      dprint,dlevel=3,side,fto,total(CR * de_meas)
    endfor
  endfor
end


; ADC bin Response  MATRIX
pro mvn_sep_response_bin_matrix_plot_old,r,window=win,face=face
  if keyword_set(win) then     wi,win,/show,wsize=[500,800],icon=0
  xrange = r.xbinrange
  if n_elements(face) eq 0 then face=0
  face_str = (['Aft','Both','Front'])[face+1]

  title= r.desc+' '+r.particle_name+' ('+r.mapname+') '+face_str
  zrange = minmax(float(r.bin3[*,*,0:255]) ,/pos)
  str_element,r,'fdesc',subtitle
  options,lim,xlog=1,xrange=xrange,yrange=[-1,260],/xstyle,/ystyle,xmargin=[10,10],/zlog,zrange=zrange,/no_interp,xtitle='Incident Energy (keV)',ytitle='Bin Number',title=title,subtitle=subtitle
  ;if not keyword_set(ok1) then ok1 = 1
  if keyword_set(face) then z = reform( r.bin3[*, face lt 0, *] ) else z = total(/pres,r.bin3,2)
  specplot,r.e_inc,r.bin_val,z,limit=lim
  ;bmap = mvn_sep_lut2map(mapnum=r.mapnum)
  bmap = mvn_sep_lut2map(mapnum=r.mapnum,sensor=r.sensornum)
  ;bmap = r.bmap
  for tid=0,1 do begin
    for fto=1,7 do begin
      w = where(bmap.tid eq tid and bmap.fto eq fto,nw)
      if nw gt 0 then begin
        b = bmap[w].bin
        bmap0 = bmap[w[0]]
        oplot,b*0.+ xrange[0]*1.5,b,psym=bmap0.psym,symsize=.5,color=bmap0.color
        xyouts,xrange[0]*1.5,average(b),' '+bmap0.name,color=bmap0.color
      endif
    endfor
  endfor

end


; ADC bin Response  MATRIX (transposed)
pro mvn_sep_response_bin_matrix_plot,r,window=win,face=face,transpose=transpose,overplot=overplot,energy_range=ei_range,zlog=zlog
  if n_elements(zlog) eq 0 then zlog=1
  if ~keyword_set(ei_range) then ei_range = minmax(r.e_inc)
  bin_range = [-2,260]
  if n_elements(face) eq 0 then face=0
  face_str = (['Aft','Both','Front'])[face+1]
  atten_str = (['Open','Closed'])[r.attenuator]
  SEP  = (['???','SEP1','SEP2'])[r.sensornum]
  title= r.desc+' '+r.particle_name+' ('+r.mapname+' '+atten_str+' '+sep+') '+face_str
  resp_matrix = float(r.bin3[*,*,0:255] )   * (r.sim_area /100 / r.nd * 3.14)
  zrange = minmax(resp_matrix ,/pos)
  if keyword_set(face) then z = reform( resp_matrix[*, face lt 0, *] ) else z = total(/pres,resp_matrix,2)
  str_element,r,'fdesc',subtitle
  if keyword_set(transpose) then begin
    options,lim,ylog=1,xrange=bin_range,yrange=ei_range,/xstyle,/ystyle,xmargin=[10,10],zlog=zlog,zrange=zrange,/no_interp,ytitle='Incident Energy (keV)',xtitle='Bin Number',title=title
    if keyword_set(win) then     wi,win,/show,wsize=[1300,800],icon=0
    x = indgen(256)
    y = r.e_inc
    z = transpose(z)
  endif else begin
    options,lim,xlog=1,xrange=ei_range,yrange=bin_range,/xstyle,/ystyle,xmargin=[10,10],zlog=zlog,zrange=zrange,/no_interp,xtitle='Incident Energy (keV)',ytitle='Bin Number',title=title,subtitle=subtitle
    if keyword_set(win) then     wi,win,/show,wsize=[500,800],icon=0
    y = indgen(256)
    x = r.e_inc
  endelse
  ;if not keyword_set(ok1) then ok1 = 1
  specplot,x,y,z,limit=lim
  overplot=get_plot_state()
  bmap = mvn_sep_get_bmap(r.mapnum,r.sensornum)
  ;bmap = r.bmap
  labpos1 = 5.
  labpos2 = 8.
  for tid=0,1 do begin
    for fto=1,7 do begin
      w = where(bmap.tid eq tid and bmap.fto eq fto,nw)
      if nw gt 0 then begin
        b = bmap[w].bin
        bmap0 = bmap[w[0]]
        if keyword_set(transpose) then begin
          oplot,b,b*0.+ labpos1,psym=bmap0.psym,symsize=.5,color=bmap0.color
          xyouts,average(b),labpos2,' '+bmap0.name,color=bmap0.color  ,align=.5
        endif else begin
          oplot,b*0.+ ei_range[0]*1.5,b,psym=bmap0.psym,symsize=.5,color=bmap0.color
          xyouts,ei_range[0]*1.5,average(b),' '+bmap0.name,color=bmap0.color
        endelse
      endif
    endfor
  endfor

end


pro mvn_sep_response_each_bin,r,bins=bins0,window=win,ylog=ylog,omega=omega
  if n_elements(omega) eq 0 then omega=0
  if keyword_set(win) then     wi,win,/show,wsize=[1200,500],icon=0
  ;bmap = mvn_sep_lut2map(mapnum=r.mapnum)
  bmap =r.bmap
  einc = r.e_inc
  wght = 1/einc^2
  yrange = minmax(r.bin3,pos=ylog)
  title= r.desc+' ('+r.mapname+')'
  plot,/nodata,minmax(einc,/pos),yrange,/xlog,ylog=ylog,xtitle='Incident Energy (keV)',ytitle='GF',title=title
  bins = keyword_set(bins0) ? bins0 : indgen(256)
  for i=0,n_elements(bins)-1 do begin
    bin=bins[i]
    dgf = r.bin3[*,omega,bin]
    oplot,einc,dgf > yrange[0]/2.,color=bmap[bin].color,psym=-bmap[bin].psym,symsize=1
    if keyword_set(bins0) then begin
      dgf_max = max(dgf,emx)
      einc_avg = average(dgf*einc*wght)/average(dgf*wght)
      einc_max = einc[emx]
      txt = '!c'+bmap[bin].name+'!c'+strtrim(bin,2)
      xyouts,einc_avg,yrange[1],txt,align=.5,color=bmap[bin].color,charsize=1
    endif
  endfor

end

pro mvn_sep_response_each_bin_GF,r,bins=bins0,sbins=sbins,window=win,ylog=ylog,overplot=overplot,omega=omega
  if n_elements(omega) eq 0 then omega=0
  if keyword_set(win) then     wi,win,/show,wsize=[1200,600],icon=0

  bmap = r.bmap
  e_inc = r.e_inc
  de_inc = (r.xlog) ? (r.e_inc * r.xbinsize  * alog(10)) : r.xbinsize
  wght = 1/e_inc^2
  yrange = minmax(r.bin3,pos=ylog)
  yrange = [.8,1e6]
  atten_str = (['Open','Closed'])[r.attenuator]
  SEP  = (['???','SEP1','SEP2'])[r.sensornum]
  ;title= r.desc+' ('+r.mapname+')'
  title= r.desc+' '+r.particle_name+' ('+r.mapname+' '+atten_str+' '+sep+') '
  if ~keyword_set(overplot) then plot,/nodata,minmax(e_inc,/pos),yrange,ystyle=1,/xlog,ylog=ylog,xtitle='Incident Energy (keV)',ytitle='GF',title=title
  overplot = 1
  bins = keyword_set(bins0) ? bins0 : indgen(256)
  for i=0,n_elements(bins)-1 do begin
    bin=bins[i]
    dgf = r.bin3[*,omega,bin]
    oplot,e_inc,dgf > yrange[0]/2.,color=bmap[bin].color,psym=-bmap[bin].psym,symsize=1
    dgf_max = max(dgf,emx)
    einc_avg = average(dgf*e_inc*wght)/average(dgf*wght)
    einc_max = e_inc[emx]
    txt = bmap[bin].name+'!c'+strtrim(bin,2)
    GF = total(dgf )  ;* de_inc)
    oplot,[einc_avg],[gf],/psym
    if keyword_set(bins0) then begin
      xyouts,einc_avg,gf,txt,align=.5,color=bmap[bin].color,charsize=1
    endif
  endfor
  for i=0,n_elements(sbins)-1 do begin
    bin=sbins[i]
    dgf = r.bin2[*,bin]
    oplot,e_inc,dgf > yrange[0]/2.,color=bmap[bin].color,psym=-bmap[bin].psym,symsize=1
    dgf_max = max(dgf,emx)
    einc_avg = average(dgf*e_inc*wght)/average(dgf*wght)
    einc_max = e_inc[emx]
    txt = bmap[bin].name+'!c'+strtrim(bin,2)
    GF = total(dgf )  ;* de_inc)
    oplot,[einc_avg],[gf],/psym
    ;   if keyword_set(bins0) then begin
    xyouts,einc_avg,gf,txt,align=.5,color=bmap[bin].color,charsize=1
    ;   endif
  endfor

end




function mvn_sep_response_spectra,r,fluxfunc
  message,'Obsolete'
  ;bmap =  mvn_sep_lut2map(lut=r.lut)
  ;bmap =  mvn_sep_lut2map(mapnum=r.mapnum)
  bmap = r.bmap
  scl = [ [[0.],[0.]] ,1/ r.adc_scale ]
  remap = [0,1,2,1,3,0,3,2]
  ftot_scale = scl[remap,*]           ; X O T OT F FO FT FTO
  ftot_scale[[3,6],*] *=2  ; OT and FT events divided by two
  ftot_scale[7,*] *= 4     ; FTO events divided by four
  escale = ftot_scale
  bmap.x = average(bmap.adc,1) * escale[bmap.fto,bmap.tid]
  bmap.dx = bmap.num * escale[bmap.fto,bmap.tid]

end

pro mvn_sep_inst_bin_response,simstat,data,new_seed=new_seed,noise_level=noise_level,mapnum=mapnum,bmap=bmap
  ;common mvn_sep_inst_bin_response_com , seed
  if size(/type,simstat) ne 8  then begin
    undefine,data
    return
  endif
  if n_elements(new_seed) ne 0 then seed = new_seed
  str_element,/add,simstat,'seed',new_seed
  n = n_elements(data)
  if n le 1 then begin
    dprint,'Must have at least 2 successful events'
    return
  endif
  ;str_element,simstat,'sensornum',sensornum
  sensornum = simstat.sensornum
  if ~keyword_set(mapnum) then str_element,simstat,'mapnum',mapnum
  if ~keyword_set(noise_level) then str_element,simstat,'noise_level',noise_level
  if n_elements(noise_level) ne 1 then noise_level=1.
  ;if n_elements(sensornum) ne 1 then sensornum=1
  if n_elements(mapnum) ne 1 then mapnum = 8
  noise_rms = noise_level * [[2.,3.,2.],[2.,3.,2.]]   ; noise O T F in kev
  adc_scale = mvn_sep_adc_calibration()
  ;adc_scale = adc_scale[*,*,sensornum]
  threshold = noise_rms * 5   ; 5 sigma threshold
  shft = [1,1,1,2,1,1,2,4]
  if keyword_set(mapnum) then begin
    lut = mvn_sep_create_lut(mapnum=mapnum)
    mapname = mvn_sep_mapnum_to_mapname(mapnum)
    bmap = mvn_sep_lut2map(lut=lut,sensor=sensornum)
    lut = fix( reform(lut,4096,2,8) )
    lut[*,*,0] = 256                 ;  not triggered (not detected)
    lut[*,*,5] = 257                 ;  FO event  (not handled correctly yet)
  endif
  one_n = replicate(1,n)
  str_element,/add,data,'fto',bytarr(2,n)
  str_element,/add,data,'em',fltarr(2,n)
  str_element,/add,data,'bin',intarr(2,n)
  for side=0,1 do begin
    ;   slabel = side ? 'B' : 'A'
    ;   ec3 = side ? data.b : data.a                          ; collected (deposited) energy in each of 3 detectors
    ec3 = data.edep[*,side]                                     ; energy deposited in each of 3 detectors for that side
    noise3 =  (noise_rms[*,side] # one_n) * randomn(seed,3,n)     ; generate noise in kev
    em3 = ec3 + noise3
    fto3 = em3 gt (threshold[*,side] # one_n)           ; determine FTO pattern  based on # above  threshold
    em3 = em3 * fto3                                             ; clear untriggered channels
    em  = total(em3,1)                                           ; total energy in all 3 triggered channels
    scl3 = (adc_scale[*,side,sensornum-1] # one_n)
    adc3 = long(em3 * scl3)  <  4095
    ftocode = reform(fto3 ## [1,2,4])                        ; This does not properly account for FO events!!
    adc = total(/pres,adc3,1) / shft[ftocode]    ; FT and OT  adc values are divided by 2,  FTO are divided by 4
    adc_bin = lut[adc,side * one_n,ftocode]
    data.fto[side] = ftocode
    data.em[side] = em
    data.bin[side] = adc_bin
    ;   if keyword_set(add) then begin
    ;     str_element,/add,data,slabel+'_em',em
    ;     str_element,/add,data,slabel+'_fto',ftocode
    ;     str_element,/add,data,slabel+'_adc',adc
    ;     str_element,/add,data,slabel+'_bin',adc_bin
    ;   endif
  endfor
  str_element,/add,simstat,'noise_level',noise_level
  str_element,/add,simstat,'mapnum',mapnum
  str_element,/add,simstat,'mapname',mapname
  ;str_element,/add,simstat,'lut',lut
  ;str_element,/add,simstat,'bmap',bmap

end



;  this function returns a true or false depending on if projection of particle direction will
;  pass through a rectangle centered at pos with edges of size width.
function mvn_sep_response_projection_filter,data,pos,width
  a = (where(width eq 0,nw))[0]   ; must be 2 for now
  if nw ne 1 then message,'At least one of width dimensions must be zero'
  if a ne 0 then message, 'Not ready yet'
  one = replicate(1,n_elements(data))
  pos_ = pos # one
  pdist= reform( ( (pos_[a] # one)- data.pos[a]) / data.dir[a]  )          ; projected distance to a-plane
  proj = data.pos + ([1,1,1] # pdist) * data.dir                  ; intersection point of dir line with rectangle's plane
  d  = sqrt(total( (pos_ - proj)^2 ,1) )                       ;  distance to center of rectangle
  ok = pdist gt 0                                                 ; use only particles moving toward rectangle
  ok = ok and proj[1,*] lt pos[1]+width[1]/2.  and proj[1,*] gt pos[1]-width[1]/2.
  ok = ok and (proj[2,*] lt pos[2]+width[2]/2.  and proj[2,*] gt pos[2]-width[2]/2.)
  ;  ok =ok and ((A_rad lt center) or (B_rad lt center))
  return,ok
end



function mvn_sep_response_data_filter,simstat,data,filter=filter, $
  ;     plus_xdir=plus_xdir,minus_xdir=minus_xdir, $
  a_side=a_side,b_side=b_side,o_det=o_det,f_det=f_det,center=center,erange=erange, $
  derange=derange, $
  detname=detname, $
  xdir=xdir,ypos=ypos, $
  fdesc=fdesc, $
  dir = dir, angle = angle,  $
  col_af=col_af,det_af=det_af,impact_bmin=impact_bmin,impact_pos=impact_pos
  filter=0
  str_element,/add,filter,'xdir',xdir
  str_element,/add,filter,'ypos',ypos
  str_element,/add,filter,'dir',dir
  str_element,/add,filter,'angle',angle
  str_element,/add,filter,'a_side',a_side
  str_element,/add,filter,'b_side',b_side
  str_element,/add,filter,'o_det',o_det
  str_element,/add,filter,'f_det',f_det
  str_element,/add,filter,'center',center
  str_element,/add,filter,'erange',erange
  str_element,/add,filter,'derange',derange
  str_element,/add,filter,'col_af',col_af
  str_element,/add,filter,'det_af',det_af
  str_element,/add,filter,'impact_bmin',impact_bmin
  str_element,/add,filter,'impact_pos',impact_pos
  str_element,/add,filter,'detname',detname

  printdat,filter,/val,output=s
  s[0]=' '
  fdesc = strjoin(strcompress(s,/remove_all),', ')
  str_element,/add,filter,'fdesc',fdesc
  n = n_elements(data)
  if n lt 1 then begin
    dprint,'At least 1 data points are required'
    return,-1
  endif
  ok = data.e_tot ge 0                   ; all true
  one = replicate(1,n_elements(data))
  if keyword_set(center) then begin
    det_width=[0,14.4,8.4]
    one = replicate(1,n_elements(data))
    A_pos = [-10,25.,-12.5] # one                              ; note: Z position is a guess.
    A_dist= ( (A_pos[0] # one)- data.pos[0]) / data.dir[2]
    A_proj = data.pos + ([1,1,1] # A_dist) * data.dir
    A_rad  = sqrt(total( (A_pos -A_proj)^2 ,1) )
    B_pos = A_pos *([-1,-1,1] # one)
    B_dist= ( (B_pos[2] # one)- data.pos[2]) / data.dir[2]
    B_proj = data.pos + ([1,1,1]# B_dist) * data.dir
    B_rad  = sqrt(total( (B_pos -B_proj)^2 ,1) )
    ok =ok and ((A_rad lt center) or (B_rad lt center))
  endif
  det_pos = [10.,-25.,-12.5]
  det_width = [0,14.4,8.4]
  ;det_width = [0,10d,10d]
  ;det_pos = [0,0,0]

  col_pos =  det_pos + [30.,0,0]
  col_width = 2 * det_width
  if keyword_set(det_Af) then  ok = ok and mvn_sep_response_projection_filter(data,det_pos,det_width)
  if keyword_set(col_Af) then  ok = ok and mvn_sep_response_projection_filter(data,col_pos,col_width)
  if keyword_set(impact_bmin) then begin
    dpos = (impact_pos # one) -data.pos
    c = crossp_trans(dpos,data.dir)
    b = sqrt(total( c^2, 1 ))
    ok = ok and (b le impact_bmin)
  endif
  if keyword_set(XDIR) then ok = ok and (data.dir[0] * xdir ge 0)
  ;if keyword_set(plus_xdir) then ok = ok and (data.dir[0] gt 0)
  ;if keyword_set(minus_xdir) then ok = ok and (data.dir[0] lt 0)
  if keyword_set(YPOS) then ok = ok and (data.pos[1] * YPOS ge 0)
  angt = 10d   ; 10 degrees
  if keyword_set(angle) then angt = angle
  if keyword_set(DIR)  then begin
    ang = acos( total( data.dir * (dir # one ),1 )/sqrt(total(dir^2.))  )  *180 /!dpi
    ok = ok and (ang lt angt)
  endif
  if keyword_set(a_side) then ok = ok and (data.pos[1] lt 0)
  if keyword_set(b_side) then ok = ok and (data.pos[1] gt 0)
  if keyword_set(F_det)  then ok = ok and (data.pos[0]*data.pos[1] lt 0)
  if keyword_set(O_det)  then ok = ok and (data.pos[0]*data.pos[1] gt 0)
  if keyword_set(erange) then ok = ok and (data.einc ge erange[0] and data.einc lt erange[1])
  if keyword_set(derange) then ok = ok and (data.e_tot ge derange[0] and data.e_tot lt derange[1])
  if keyword_set(detname) then begin
    ;   str_element,simstat,'bmap',bmap
    bmap = simstat.bmap
    ;   if ~keyword_set(bmap) then   bmap = mvn_sep_lut2map(lut=simstat.lut)
    bins = where( strmatch(bmap.name,detname) ,nbins)
    ok1 = 0
    for side=0,1 do for b=0,nbins-1 do ok1 = ok1 or (data.bin[side] eq bins[b])
    ok = ok and ok1
  endif
  w_ok = where(ok,n_ok)
  npart =0
  str_element,simstat,'npart',npart
  s[0]=strjoin(strtrim([n_ok,n,npart],1),'/')
  fdesc = strjoin(strcompress(s,/remove_all),', ')
  str_element,/add,filter,'fdesc',fdesc


  return,ok
end




function mvn_sep_inst_response,simstat,data0,mapnum=mapnum ,noise_level=noise_level,filter=filter ,bmap=bmap
  if n_elements(data0) le 1 then begin
    dprint,'Must have at least 2 successful events'
    return,0
  endif

  mvn_sep_inst_bin_response,simstat,data0,mapnum=mapnum,noise_level=noise_level  ,bmap=bmap

  w= where( mvn_sep_response_data_filter(simstat,data0,_extra=filter,filter=out_filter),nw)
  if nw le 1 then begin
    dprint,'Filter leaves no data. Must have at least 2 successful events'
    ;   return,0
  endif
  if nw lt 10 then begin
    dprint,'Very few samples.',nw
  endif
  if nw ne 0 then data=data0[w]

  ;str_element,simstat,'nsuccess',nw
  ;str_element,simstat,'type',type
  ;str_element,simstat,'sensornum',sensornum
  str_element,simstat,'npart',npart
  str_element,simstat,'sim_energy_range',simrange
  str_element,simstat,'sim_energy_log',simlog
  str_element,simstat,'n_omega',n_omega
  ;str_element,simstat,'noise_level',noise_level
  str_element,simstat,'desc',desc
  str_element,simstat,'sim_area',area_mm2
  str_element,simstat,'xbinsize',xbinsize
  str_element,simstat,'ybinsize',ybinsize

  str_element,out_filter,'fdesc',fdesc

  area_cm2 = area_mm2/100.
  ;if n_elements(noise_level) ne 1 then noise_level=1.
  ;if n_elements(sensornum) ne 1 then sensornum=0
  ;if n_elements(type) ne 1 then type=0               ;  -1: electrons,  1:protons,   2:???
  if ~keyword_set(xbinsize) then xbinsize= .025
  if ~keyword_set(ybinsize) then ybinsize= xbinsize

  srange =  simlog ? alog10(simrange) : simrange
  brange=  minmax([1.,simrange] )        ;[1.,100e3]

  xlog=1
  ylog=1
  xbinrange = brange
  ybinrange = brange
  ndata = n_elements(data)
  ND = npart * xbinsize / (srange[1] - srange[0])
  nx_einc = long((xlog ? alog10(xbinrange[1]/xbinrange[0]) : xbinrange[1]-xbinrange[0]) / xbinsize)
  ny_emeas= long((ylog ? alog10(ybinrange[1]/ybinrange[0]) : ybinrange[1]-ybinrange[0]) / ybinsize)
  xs0 = xlog ? alog(xbinrange[0]) : xbinrange[0]
  ys0 = ylog ? alog(ybinrange[0]) : ybinrange[0]
  ny_bins = 258L   ; max(lut)+1
  if ~keyword_set(n_omega) then n_omega = 2L
  G4=fltarr(nx_einc,ny_emeas,8,2)
  adcbin_hist = lonarr(nx_einc,n_omega,ny_bins)
  adcbin_index= lindgen(nx_einc,n_omega,ny_bins)
  bin_val = indgen(ny_bins)
  ei_val = ( (indgen(nx_einc)+.5d) *xbinsize) + xs0
  if xlog then ei_val = 10.d^ ei_val
  em_val = ( (indgen(ny_emeas)+.5d) *ybinsize) + ys0
  if ylog then em_val = 10.d^ em_val
  if ndata ne 0 then begin
    einc = data.einc
    one_n = replicate(1,ndata)
    for side=0,1 do begin
      ftocode = data.fto[side]
      adc_bin = data.bin[side]
      em      = data.em[side]
      einc_bin = long(  (xlog ? alog10(einc/xbinrange[0]) : einc-xbinrange[0]) / xbinsize )
      omega_bin = data.dir[0] gt 0                                ;  separate into two 'angle' bins based on the x direction
      ind = adcbin_index[einc_bin,omega_bin,adc_bin]
      h = histogram(ind,binsize=1,min=0,max= nx_einc*ny_bins*n_omega-1)
      adcbin_hist += h
      for fto = 0,7 do begin
        ;       if fto eq 0 then ok = ftocode gt 0 else $
        ok = fto eq ftocode
        w = where(ok,nw)
        if nw eq 0 then begin
          dprint,dlevel=3,side,fto, ' Not encountered ',desc
          continue
        endif
        G2 = histbins2d(einc[w],em[w],ei_val2,em_val2,xbinsize=xbinsize,ybinsize=ybinsize,xrange=xbinrange,yrange=ybinrange,xlog=xlog,ylog=ylog)
        G2 = G2/ND * area_cm2 ;*!dpi; *2* 4 ; normalize to area  (cm^2) * 4pi ster
        if ~keyword_set(g4) then G4 = replicate(0.,[size(/dimen,G2),8,2] )
        G4[*,*,fto,side] = g2
      endfor
    endfor
  endif


  response= simstat
  str_element,/add,response,'bmap',bmap
  str_element,/add,response,'ndata',ndata
  str_element,/add,response,'nd',nd
  str_element,/add,response,'xlog',xlog
  str_element,/add,response,'ylog',ylog
  str_element,/add,response,'xbinsize',xbinsize
  str_element,/add,response,'ybinsize',ybinsize
  str_element,/add,response,'xbinrange',xbinrange
  str_element,/add,response,'ybinrange',ybinrange
  str_element,/add,response,'e_inc',ei_val
  str_element,/add,response,'e_meas',em_val
  str_element,/add,response,'G4',g4
  str_element,/add,response,'bin3',adcbin_hist
  str_element,/add,response,'GB3' , adcbin_hist *  (response.sim_area /100 / response.nd * 3.14)
  str_element,/add,response,'bin_val',bin_val
  peakeinc = mvn_sep_inst_response_peakeinc(response,width=30)
  str_element,/add,response,'peakeinc',peakeinc
  str_element,/add,response,'fdesc',fdesc
  str_element,/add,response,'filter',out_filter

  return,response
end


pro mvn_sep_response_omega_plot,window=win,data,_extra=ex,filter=filter,binscale=binscale,simstat=simstat,posflag=posflag
  if n_elements(data) le 1 then begin
    dprint,'At least 2 data points are required'
    ;   return
  endif
  if keyword_set(win) then     wi,win,/show,wsize=round([800,400]*1.),icon=0
  ;printdat,ex
  ok = mvn_sep_response_data_filter(simstat,data,_extra=ex,filter=filter,fdesc=fdesc)

  title='Angular Direction'
  desc = '-'
  particle = '-'
  if keyword_set(posflag) then title='Position Angle'
  str_element,simstat,'title',title
  str_element,simstat,'desc',desc
  str_element,simstat,'particle_name',particle
  title = title+' '+desc+' '+particle
  subtitle = fdesc

  ;if keyword_set(simstat) then begin
  ;   npart=-1
  ;   str_element,simstat,'npart',npart
  ;   n = total(ok)
  ;   dprint,n ,' of ',npart
  ;endif
  w = where(ok,nw)
  erase
  ;if nw eq 0 then return
  xrange = [0,360]
  yrange = [-1,1.]
  !p.multi = [0,2,2]
  options,lim,xtitle='Phi (degrees)',ytitle='sin(theta) ',title=title,xmargin=[10,10],xrange=xrange,yrange=yrange,/xstyle,/ystyle,subtitle=fdesc
  if nw ne 0 then begin
    dir = keyword_set(posflag) ? data[w].pos : data[w].dir
    xyz_to_polar,transpose(dir),theta=th,phi=ph,mag=r,/ph_0_360
    sth = sind(th)
  endif
  plot,_extra=lim,[0,1],/nodata
  if keyword_set(binscale) then  begin    ;binscale =1
    xbinsize = binscale / !x.s[1] / !d.x_size
    ybinsize = binscale / !y.s[1] / !d.y_size
  endif
  zv = histbins2d(ph,sth,xv,yv,xbinsize=xbinsize,ybinsize=ybinsize,xrange=xrange,yrange=yrange)
  ;printdat,xbinsize,ybinsize
  ;printdat,yv,minmax(yv)
  options,lim,/zlog,zrange=minmax(/pos,[zv[*],1.,10])    ,/no_interp
  specplot,xv,yv,zv,limit=lim,/overplot
  ;   if keyword_set(binscale) && binscale ge 10 then oplot,_extra=lim,ph,sind(th),psym=3;,/nodata

  plot,total(zv,1),yv,yrange=yrange,xmargin=[10,10] ,/ystyle  ;,xrange=minmax(tot
  plot,xv,total(zv,2),xrange=xrange,xmargin=[10,10],/xstyle

  ;directions
  if nw ne 0 then begin
    if keyword_set(posflag) then   range=minmax(dir) else range=[-1,1]
    plot,xb,histbins(dir,xb,range=range),/nodata
    nb = n_elements(xb)
    oplot,xb,histbins(dir[0,*],xb,range=range,/shift,nbins=nb),color=2,psym=10
    oplot,xb,histbins(dir[1,*],xb,range=range,/shift,nbins=nb),color=4,psym=10
    oplot,xb,histbins(dir[2,*],xb,range=range,/shift,nbins=nb),color=6,psym=10
  endif

  if 0 then begin
    specplot,xv,yv,zv,limit=lim;,/overplot
    if binscale ge 10 then oplot,_extra=lim,ph,sind(th),psym=3;,/nodata
  endif
  !p.multi=0



end



pro mvn_sep_response_aperture_plot,window=win,data,_extra=ex,filter=filter,binscale=binscale,simstat=simstat
  ok = mvn_sep_response_data_filter(simstat,data,_extra=ex,filter=filter)
  w = where(ok,nw)
  ;if nw le 1 then begin
  ;  dprint,'Not enough data'
  ;  return
  ;endif
  title='Start Position'
  desc = ''
  particle = ''
  str_element,simstat,'title',title
  str_element,simstat,'desc',desc
  str_element,simstat,'particle_name',particle
  title = title+' '+desc+' '+particle
  str_element,filter,'fdesc',subtitle
  if keyword_set(win) then     wi,win,/show,wsize=round([650,330]*1.),icon=0
  ;w = where(ok,nw)
  ;dprint,nw
  ;printdat,ex,filter
  xrange = minmax(data.pos[1])
  yrange = minmax(data.pos[2])
  options,lim,xtitle='Y position (mm)',ytitle='Z position (mm)',title=title,xmargin=[10,10],/isotropic,xrange=xrange,yrange=yrange,/xstyle,/ystyle,subtitle=subtitle
  plot,_extra=lim,[0,1],/nodata
  if nw eq 0 then begin
    dprint, 'No selected data!'
    ;     erase
    return
  endif
  ;    xp = data[w].pos[0]
  ;    yp = data[w].pos[1]
  if keyword_set(binscale) then begin                          ; binscale =1 gives one bin per pixel
    xbinsize = binscale / !x.s[1] / !d.x_size
    ybinsize = binscale / !y.s[1] / !d.y_size
  endif
  zv = histbins2d(data[w].pos[1],data[w].pos[2],xv,yv,xbinsize=xbinsize,ybinsize=ybinsize)
  options,lim,/zlog,zrange=minmax(/pos,zv)    ,/no_interp
  specplot,xv,yv,zv,limit=lim
  ;   printdat,xv,yv,zv,xbinsize,ybinsize
end



pro makeallpngs,name
  device,window_state=ws
  for i =0,20 do  if ws[i] then  makepng,name+'_'+strtrim(i,2),window=i

end

pro mvn_sep_response_plots,simstat,data,filter=f,window=win             ;,mapnum=mapnum,noise_level=noise_level,seed=seed
  if ~ keyword_set(simstat) || ~ keyword_set(data) then return
  if not keyword_set(win) then win=1
  binscale=3
  mvn_sep_response_aperture_plot,data,simstat=simstat,window= win++,_extra=f ,binscale=binscale
  mvn_sep_response_omega_plot,data,simstat=simstat,window=win++,_extra=f ,binscale=binscale, /posflag
  mvn_sep_response_omega_plot,data,simstat=simstat,window=win++,_extra=f ,binscale=binscale

  ;we= where( mvn_sep_response_data_filter(simstat,data,_extra=f,filter=f2),nwe)
  resp = mvn_sep_inst_response(simstat,data,filter=f)
  ;printdat,resp
  if ~keyword_set(resp) then stop

  mvn_sep_response_matrix_plots,resp,window=win++
  mvn_sep_response_matrix_plots,resp,window=win++,single=4
  ;mvn_sep_response_bin_matrix_plot,resp,window=win++ ,face=0         ; both faces
  mvn_sep_response_bin_matrix_plot,resp,window=win++ ,face=-1
  mvn_sep_response_bin_matrix_plot,resp,window=win++ ,face=+1
  mvn_sep_response_plot_gf,resp,window=win++,/ylog;,xrange=[1e2,1e6]

end






; Begin program here

if 0 then $
  allfiles = mvn_sep_inst_response_retrieve('results2/results/*/mvn_sep_*.dat',age_limit=900)

;testrun = '4pi_magcorrect'
;testrun = 'Geom1'
;testrun = '4pi'
testrun = 'oxygen'
testrun = 'oxygen_dir'
;testrun = 'run1b'
;testrun = 'Geom1_front'
;testrun = 'Geom1_back'
;testrun = '4pi_sim'
testrun = '4pi_run05_sep1'
testrun = '4pi_run05_sep2'
;testrun = 0
mapnum=9

if not keyword_set(testrun) then testrun = 'run04_sep2'
undefine,resp_e0,resp_p0,resp_g0,resp_e1,resp_p1,resp_g1,bmap

if ~keyword_set(ltestrun) || ltestrun ne testrun then begin
  undefine,  simstat_e,  simstat_p ,  simstat_g ,  simstat_Ox
endif
ltestrun =testrun
case testrun of

  'run1': begin
    if ~keyword_set(data_e1) then mapnum_last=0
    mvn_sep_response_simulation1_load,-1,data=data_e1,simstat=simstat_e  ; load electrons (but only if data_e1 is not defined
    mvn_sep_response_simulation1_load, 1,data=data_p1,simstat=simstat_p   ; load ions
  end

  'run1b': begin
    mvn_sep_read_mult_sim_files,simstat_e,data_e1,pathnames=['results/mvn_sep_e-_AF.dat','results/mvn_sep_e-_AO.dat'],type=-1,/dosymm
    mvn_sep_read_mult_sim_files,simstat_p,data_p1,pathnames=['results/mvn_sep_proton_AF.dat','results/mvn_sep_proton_AO.dat'],type=+1,/dosymm
    options,simstat_e,xbinsize=0.05
    options,simstat_p,xbinsize=0.05
  end


  '4pi_sim': begin
    mvn_sep_response_simdata_rand,-1,data=data_e1,simstat=simstat_e,seed=seed  ;,window=win
  end

  'high_dens_detectors': begin
    mvn_sep_read_mult_sim_files,simstat_e,data_e1,pathnames='results2/results/KoreaDetector/run01/mvn_sep2_e-_seed01_open_AF.dat',type=-1
    ;mvn_sep_read_mult_sim_files,simstat_e,data_e1,pathnames='results2/results/KoreaDetector/run01/mvn_sep2_e-_seed01_open_AO.dat',type=-1
    ;mvn_sep_read_mult_sim_files,simstat_p,data_p1,pathnames='results2/results/KoreaDetector/run01/mvn_sep2_proton_seed01_open_AF.dat',type=-1
    ;mvn_sep_read_mult_sim_files,simstat_p,data_p1,pathnames='results2/results/KoreaDetector/run01/mvn_sep2_proton_seed01_open_AO.dat',type=-1
  end


  'Geom1_front': begin
    mvn_sep_read_mult_sim_files,simstat_e,data_e1,pathnames='results2/results/geometric_factor/mvn_sep_e-_.dat',type=-1
    mvn_sep_read_mult_sim_files,simstat_p,data_p1,pathnames='results2/results/geometric_factor/mvn_sep_proton_.dat',type=+1
    mvn_sep_read_mult_sim_files,simstat_g,data_g,pathnames='results2/results/geometric_factor/mvn_sep_gamma_open_.dat',type=0
    options,simstat_e,xbinsize=0.1
    options,simstat_p,xbinsize=0.1
    options,simstat_g,xbinsize=0.1
  end

  'Geom1_back': begin
    mvn_sep_read_mult_sim_files,simstat_e,data_e1,pathnames='results2/results/geometric_factor/mvn_sep_e-_back_side_open_.dat',type=-1
    mvn_sep_read_mult_sim_files,simstat_p,data_p1,pathnames='results2/results/geometric_factor/mvn_sep_proton_back_side_open_.dat',type=+1
    mvn_sep_read_mult_sim_files,simstat_g,data_g,pathnames='results2/results/geometric_factor/mvn_sep_gamma_back_side_open_.dat',type=0
    options,simstat_e,xbinsize=0.1
    options,simstat_p,xbinsize=0.1
    options,simstat_g,xbinsize=0.1
  end

  'Geom1_both': begin
    mvn_sep_read_mult_sim_files,simstat_e,data_e1,pathnames='results2/results/geometric_factor/'+['mvn_sep_e-_.dat','mvn_sep_e-_back_side_open_.dat'],type=-1
    mvn_sep_read_mult_sim_files,simstat_p,data_p1,pathnames='results2/results/geometric_factor/'+['mvn_sep_proton_.dat','mvn_sep_proton_back_side_open_.dat'],type=+1
    mvn_sep_read_mult_sim_files,simstat_g,data_g,pathnames='results2/results/geometric_factor/'+['mvn_sep_gamma_open_.dat','mvn_sep_gamma_back_side_open_.dat'],type=0
    options,simstat_e,xbinsize=0.1
    options,simstat_p,xbinsize=0.1
    options,simstat_g,xbinsize=0.1
  end

  'closed_back': begin
    mvn_sep_read_mult_sim_files,simstat_e,data_e1,pathnames='results2/results/geometric_factor/mvn_sep_e-_back_side_closed_.dat',type=-1
    mvn_sep_read_mult_sim_files,simstat_p,data_p1,pathnames='results2/results/geometric_factor/mvn_sep_proton_back_side_closed_.dat',type=+1
    mvn_sep_read_mult_sim_files,simstat_g,data_g,pathnames='results2/results/geometric_factor/mvn_sep_gamma_back_side_closed_.dat',type=0
    options,simstat_e,xbinsize=0.1
    options,simstat_p,xbinsize=0.1
    options,simstat_g,xbinsize=0.1
  end

  'closed_front': begin
    mvn_sep_read_mult_sim_files,simstat_e,data_e1,pathnames='results2/results/geometric_factor/mvn_sep_e-_closed_.dat',type=-1
    mvn_sep_read_mult_sim_files,simstat_p,data_p1,pathnames='results2/results/geometric_factor/mvn_sep_proton_closed_.dat',type=+1
    mvn_sep_read_mult_sim_files,simstat_g,data_g,pathnames='results2/results/geometric_factor/mvn_sep_gamma_closed_.dat',type=0
    options,simstat_e,xbinsize=0.1
    options,simstat_p,xbinsize=0.1
    options,simstat_g,xbinsize=0.1
  end

  'run2_iso':begin
    mvn_sep_read_mult_sim_files,simstat_e,data_e1,pathnames='results2/results/run2_isotropic/mvn_sep_e-_.dat',type=-1
    mvn_sep_read_mult_sim_files,simstat_p,data_p1,pathnames='results2/results/run2_isotropic/mvn_sep_proton_.dat',type=+1
    options,simstat_e,xbinsize=1.d/5
    options,simstat_p,xbinsize=1.d/5
  end

  '4pi_open':begin
    ;mvn_sep_read_mult_sim_files,simstat_e,data_e1,pathnames='results2/results/4PI/run01/mvn_sep_e-_all.dat',type=-1
    ;mvn_sep_read_mult_sim_files,simstat_e,data_e1,pathnames='results2/results/4PI/run01/mvn_sep_e-_successful.dat',type=-1
    mvn_sep_read_mult_sim_files,simstat_e,data_e1,pathnames='results2/results/4PI/run01/mvn_sep_e-_seed??_open.dat',type=-1
    mvn_sep_read_mult_sim_files,simstat_p,data_p1,pathnames='results2/results/4PI/run01/mvn_sep_proton_seed??_open.dat',type=+1
    mvn_sep_read_mult_sim_files,simstat_g,data_g ,pathnames='results2/results/4PI/run01/mvn_sep_gamma_seed??_open.dat',type=0

    ;mvn_sep_read_mult_sim_files,simstat_p,data_p1,pathnames='results2/results/run2_isotropic/mvn_sep_proton_.dat',type=+1
    options,simstat_e,xbinsize=1.d/10
    options,simstat_p,xbinsize=1.d/10
    options,simstat_g,xbinsize=1.d/10
  end

  '4pi_magcorrect':begin
    mvn_sep_read_mult_sim_files,simstat_e,data_e1,pathnames='results2/results/run2_isotropic/mvn_sep_e-magneticfield_corrected_.dat',type=-1
    mvn_sep_read_mult_sim_files,simstat_p,data_p1,pathnames='results2/results/run2_isotropic/mvn_sep_protonmagneticfield_corrected_.dat',type=+1
    options,simstat_e,xbinsize=1.d/5
    options,simstat_p,xbinsize=1.d/5
  end

  'oxygen':begin  ;
    mvn_sep_read_mult_sim_files,simstat_Ox,data_Ox,pathnames='results2/results/oxygen/mvn_sep_oxygen_.dat',type=2,dosymm=1
    str_element,/add,simstat_Ox,'desc','Oxygen Sim'   ; Needs checking!
    str_element,/add,simstat_Ox,'sensornum',2   ; Needs checking!
    str_element,/add,simstat_Ox,'sim_energy_log',1
    str_element,/add,simstat_Ox,'sim_energy_range',[10.,1e7]
    resp_O0 = mvn_sep_inst_response(simstat_Ox,data_Ox,mapnum=mapnum,bmap=bmap)
  end

  'oxygen_dir':begin  ;
    mvn_sep_read_mult_sim_files,simstat_Ox,data_Ox,pathnames='results2/results/oxygen/mvn_sep_oxygen_randomDir_.dat',type=2,dosymm=1
    str_element,/add,simstat_Ox,'desc','Ox GEOM'   ; Needs checking!
    str_element,/add,simstat_Ox,'particle_name','Oxygen'   ; Needs checking!
    str_element,/add,simstat_Ox,'desc','GEOM'   ; Needs checking!
    str_element,/add,simstat_Ox,'sensornum',2   ; Needs checking!
    str_element,/add,simstat_Ox,'sim_energy_log',1
    str_element,/add,simstat_Ox,'sim_energy_range',[10.,1e7]
    if simstat_ox.sim_area eq 0 then simstat_ox.sim_area= 2 * !pi * 15. ^2   ;  2 circles
    resp_O0 = mvn_sep_inst_response(simstat_Ox,data_Ox,mapnum=mapnum,bmap=bmap)
  end

  'geom1_all_elec':begin  ;
    mvn_sep_read_mult_sim_files,simstat_e,data_e1,pathnames='results2/results/geometric_factor/mvn_sep_e-_allEvents_.dat',type=-1
  end

  'gamma':begin  ;
    mvn_sep_read_mult_sim_files,simstat_g,data_g,pathnames='results2/results/geometric_factor/mvn_sep_gamma_back_side_open_.dat',type=0
  end


  'run03_sep1':begin
    mvn_sep_read_mult_sim_files,simstat_e,pathnames='results2/results/4PI/run03/mvn_sep1_e-_atten0_seed??_.dat',type=-1,data_e
    mvn_sep_read_mult_sim_files,simstat_p,pathnames='results2/results/4PI/run03/mvn_sep1_proton_atten0_seed??_.dat',type=+1,data_p
    mvn_sep_read_mult_sim_files,simstat_g,pathnames='results2/results/4PI/run03/mvn_sep1_gamma_atten0_seed??_.dat',type=0,data_g
    resp_e0 = mvn_sep_inst_response(simstat_e,data_e,mapnum=mapnum,bmap=bmap)
    resp_p0 = mvn_sep_inst_response(simstat_p,data_p,mapnum=mapnum,bmap=bmap)
    resp_g0 = mvn_sep_inst_response(simstat_g,data_g,mapnum=mapnum,bmap=bmap)
    mvn_sep_read_mult_sim_files,simstat_e1,data_e1,pathnames='results2/results/4PI/run03/mvn_sep1_e-_atten1_seed??_.dat',type=-1
    mvn_sep_read_mult_sim_files,simstat_p1,data_p1,pathnames='results2/results/4PI/run03/mvn_sep1_proton_atten1_seed??_.dat',type=+1
    mvn_sep_read_mult_sim_files,simstat_g1,data_g1 ,pathnames='results2/results/4PI/run03/mvn_sep1_gamma_atten1_seed??_.dat',type=0
    resp_e1 = mvn_sep_inst_response(simstat_e1,data_e1,mapnum=mapnum,bmap=bmap)
    resp_p1 = mvn_sep_inst_response(simstat_p1,data_p1,mapnum=mapnum,bmap=bmap)
    resp_g1 = mvn_sep_inst_response(simstat_g1,data_g1,mapnum=mapnum,bmap=bmap)
  end


  'run03_sep2':begin
    mvn_sep_read_mult_sim_files,simstat_e,pathnames='results2/results/4PI/run03/mvn_sep2_e-_atten0_seed??_.dat',type=-1,data_e
    mvn_sep_read_mult_sim_files,simstat_p,pathnames='results2/results/4PI/run03/mvn_sep2_proton_atten0_seed??_.dat',type=+1,data_p
    mvn_sep_read_mult_sim_files,simstat_g,pathnames='results2/results/4PI/run03/mvn_sep2_gamma_atten0_seed??_.dat',type=0,data_g
    resp_e0 = mvn_sep_inst_response(simstat_e,data_e,mapnum=mapnum,bmap=bmap)
    resp_p0 = mvn_sep_inst_response(simstat_p,data_p,mapnum=mapnum,bmap=bmap)
    resp_g0 = mvn_sep_inst_response(simstat_g,data_g,mapnum=mapnum,bmap=bmap)
    mvn_sep_read_mult_sim_files,simstat_e1,data_e1,pathnames='results2/results/4PI/run03/mvn_sep2_e-_atten1_seed??_.dat',type=-1
    mvn_sep_read_mult_sim_files,simstat_p1,data_p1,pathnames='results2/results/4PI/run03/mvn_sep2_proton_atten1_seed??_.dat',type=+1
    mvn_sep_read_mult_sim_files,simstat_g1,data_g1 ,pathnames='results2/results/4PI/run03/mvn_sep2_gamma_atten1_seed??_.dat',type=0
    resp_e1 = mvn_sep_inst_response(simstat_e1,data_e1,mapnum=mapnum,bmap=bmap)
    resp_p1 = mvn_sep_inst_response(simstat_p1,data_p1,mapnum=mapnum,bmap=bmap)
    resp_g1 = mvn_sep_inst_response(simstat_g1,data_g1,mapnum=mapnum,bmap=bmap)
  end


  'run04_sep1':begin
    mvn_sep_read_mult_sim_files,simstat_e,pathnames='results2/results/4PI/run04/mvn_sep1_e-_atten0_seed??_.dat',type=-1,data_e
    mvn_sep_read_mult_sim_files,simstat_p,pathnames='results2/results/4PI/run04/mvn_sep1_proton_atten0_seed??_.dat',type=+1,data_p
    ;mvn_sep_read_mult_sim_files,simstat_g,pathnames='results2/results/4PI/run04/mvn_sep1_gamma_atten0_seed??_.dat',type=0,data_g
    resp_e0 = mvn_sep_inst_response(simstat_e,data_e,mapnum=mapnum,bmap=bmap)
    resp_p0 = mvn_sep_inst_response(simstat_p,data_p,mapnum=mapnum,bmap=bmap)
    ;resp_g0 = mvn_sep_inst_response(simstat_g,data_g,mapnum=mapnum,bmap=bmap)
    mvn_sep_read_mult_sim_files,simstat_e1,data_e1,pathnames='results2/results/4PI/run04/mvn_sep1_e-_atten1_seed??_.dat',type=-1
    mvn_sep_read_mult_sim_files,simstat_p1,data_p1,pathnames='results2/results/4PI/run04/mvn_sep1_proton_atten1_seed??_.dat',type=+1
    ;mvn_sep_read_mult_sim_files,simstat_g1,data_g1 ,pathnames='results2/results/4PI/run04/mvn_sep1_gamma_atten1_seed??_.dat',type=0
    resp_e1 = mvn_sep_inst_response(simstat_e1,data_e1,mapnum=mapnum,bmap=bmap)
    resp_p1 = mvn_sep_inst_response(simstat_p1,data_p1,mapnum=mapnum,bmap=bmap)
    ;resp_g1 = mvn_sep_inst_response(simstat_g1,data_g1,mapnum=mapnum,bmap=bmap)
  end



  'run04_sep2':begin
    mvn_sep_read_mult_sim_files,simstat_e,pathnames='results2/results/4PI/run04/mvn_sep2_e-_atten0_seed??_.dat',type=-1,data_e
    mvn_sep_read_mult_sim_files,simstat_p,pathnames='results2/results/4PI/run04/mvn_sep2_proton_atten0_seed??_.dat',type=+1,data_p
    ;mvn_sep_read_mult_sim_files,simstat_g,pathnames='results2/results/4PI/run04/mvn_sep2_gamma_atten0_seed??_.dat',type=0,data_g
    resp_e0 = mvn_sep_inst_response(simstat_e,data_e,mapnum=mapnum,bmap=bmap)
    resp_p0 = mvn_sep_inst_response(simstat_p,data_p,mapnum=mapnum,bmap=bmap)
    ;resp_g0 = mvn_sep_inst_response(simstat_g,data_g,mapnum=mapnum,bmap=bmap)
    mvn_sep_read_mult_sim_files,simstat_e1,data_e1,pathnames='results2/results/4PI/run04/mvn_sep2_e-_atten1_seed??_.dat',type=-1
    mvn_sep_read_mult_sim_files,simstat_p1,data_p1,pathnames='results2/results/4PI/run04/mvn_sep2_proton_atten1_seed??_.dat',type=+1
    ;mvn_sep_read_mult_sim_files,simstat_g1,data_g1 ,pathnames='results2/results/4PI/run04/mvn_sep2_gamma_atten1_seed??_.dat',type=0
    resp_e1 = mvn_sep_inst_response(simstat_e1,data_e1,mapnum=mapnum,bmap=bmap)
    resp_p1 = mvn_sep_inst_response(simstat_p1,data_p1,mapnum=mapnum,bmap=bmap)
    ;resp_g1 = mvn_sep_inst_response(simstat_g1,data_g1,mapnum=mapnum,bmap=bmap)
  end


  'run02B_sep1':begin
    mvn_sep_read_mult_sim_files,simstat_e,pathnames='g4work/sep/results/run02/mvn_sep1_e-_atten0_seed??_.dat',type=-1,data_e
    mvn_sep_read_mult_sim_files,simstat_p,pathnames='g4work/sep/results/run02/mvn_sep1_proton_atten0_seed??_.dat',type=+1,data_p
    ;mvn_sep_read_mult_sim_files,simstat_g,pathnames='g4work/sep/results/run02/mvn_sep1_gamma_atten0_seed??_.dat',type=0,data_g
    resp_e0 = mvn_sep_inst_response(simstat_e,data_e,mapnum=mapnum,bmap=bmap)
    resp_p0 = mvn_sep_inst_response(simstat_p,data_p,mapnum=mapnum,bmap=bmap)
    ;resp_g0 = mvn_sep_inst_response(simstat_g,data_g,mapnum=mapnum,bmap=bmap)
    mvn_sep_read_mult_sim_files,simstat_e1,data_e1,pathnames='g4work/sep/results/run02/mvn_sep1_e-_atten1_seed??_.dat',type=-1
    mvn_sep_read_mult_sim_files,simstat_p1,data_p1,pathnames='g4work/sep/results/run02/mvn_sep1_proton_atten1_seed??_.dat',type=+1
    ;mvn_sep_read_mult_sim_files,simstat_g1,data_g1 ,pathnames='g4work/sep/results/run02/mvn_sep1_gamma_atten1_seed??_.dat',type=0
    resp_e1 = mvn_sep_inst_response(simstat_e1,data_e1,mapnum=mapnum,bmap=bmap)
    resp_p1 = mvn_sep_inst_response(simstat_p1,data_p1,mapnum=mapnum,bmap=bmap)
    ;resp_g1 = mvn_sep_inst_response(simstat_g1,data_g1,mapnum=mapnum,bmap=bmap)
  end

  'run02B_sep2':begin
    mvn_sep_read_mult_sim_files,simstat_e,pathnames='g4work/sep/results/run02/mvn_sep2_e-_atten0_seed??_.dat',type=-1,data_e
    mvn_sep_read_mult_sim_files,simstat_p,pathnames='g4work/sep/results/run02/mvn_sep2_proton_atten0_seed??_.dat',type=+1,data_p
    ;mvn_sep_read_mult_sim_files,simstat_g,pathnames='g4work/sep/results/run02/mvn_sep2_gamma_atten0_seed??_.dat',type=0,data_g
    resp_e0 = mvn_sep_inst_response(simstat_e,data_e,mapnum=mapnum,bmap=bmap)
    resp_p0 = mvn_sep_inst_response(simstat_p,data_p,mapnum=mapnum,bmap=bmap)
    ;resp_g0 = mvn_sep_inst_response(simstat_g,data_g,mapnum=mapnum,bmap=bmap)
    mvn_sep_read_mult_sim_files,simstat_e1,data_e1,pathnames='g4work/sep/results/run02/mvn_sep2_e-_atten1_seed??_.dat',type=-1
    mvn_sep_read_mult_sim_files,simstat_p1,data_p1,pathnames='g4work/sep/results/run02/mvn_sep2_proton_atten1_seed??_.dat',type=+1
    ;mvn_sep_read_mult_sim_files,simstat_g1,data_g1 ,pathnames='g4work/sep/results/run02/mvn_sep2_gamma_atten1_seed??_.dat',type=0
    resp_e1 = mvn_sep_inst_response(simstat_e1,data_e1,mapnum=mapnum,bmap=bmap)
    resp_p1 = mvn_sep_inst_response(simstat_p1,data_p1,mapnum=mapnum,bmap=bmap)
    ;resp_g1 = mvn_sep_inst_response(simstat_g1,data_g1,mapnum=mapnum,bmap=bmap)
  end


  '4pi_run05_sep1':begin
    mvn_sep_read_mult_sim_files,simstat_e,pathnames='g4work/davinMaven/results/4PI/run05/mvn_sep1_e-_atten0_seed??_.dat',type=-1,data_e
    mvn_sep_read_mult_sim_files,simstat_p,pathnames='g4work/davinMaven/results/4PI/run05/mvn_sep1_proton_atten0_seed??_.dat',type=+1,data_p
    ;mvn_sep_read_mult_sim_files,simstat_g,pathnames='g4work/davinMaven/results/4PI/run05/mvn_sep1_gamma_atten0_seed??_.dat',type=0,data_g
    str_element,/add,simstat_e,'desc','4Pi'
    str_element,/add,simstat_p,'desc','4Pi'
    ;  str_element,/add,simstat_g1,'desc','4Pi'
    resp_e0 = mvn_sep_inst_response(simstat_e,data_e,mapnum=mapnum,bmap=bmap)
    resp_p0 = mvn_sep_inst_response(simstat_p,data_p,mapnum=mapnum,bmap=bmap)
    ;resp_g0 = mvn_sep_inst_response(simstat_g,data_g,mapnum=mapnum,bmap=bmap)
    mvn_sep_read_mult_sim_files,simstat_e1,data_e1,pathnames='g4work/davinMaven/results/4PI/run05/mvn_sep1_e-_atten1_seed??_.dat',type=-1
    mvn_sep_read_mult_sim_files,simstat_p1,data_p1,pathnames='g4work/davinMaven/results/4PI/run05/mvn_sep1_proton_atten1_seed??_.dat',type=+1
    ;mvn_sep_read_mult_sim_files,simstat_g1,data_g1 ,pathnames='g4work/davinMaven/results/4PI/run05/mvn_sep1_gamma_atten1_seed??_.dat',type=0
    str_element,/add,simstat_e1,'desc','4Pi'
    str_element,/add,simstat_p1,'desc','4Pi'
    ;  str_element,/add,simstat_g1,'desc','4Pi'
    resp_e1 = mvn_sep_inst_response(simstat_e1,data_e1,mapnum=mapnum,bmap=bmap)
    resp_p1 = mvn_sep_inst_response(simstat_p1,data_p1,mapnum=mapnum,bmap=bmap)
    ;resp_g1 = mvn_sep_inst_response(simstat_g1,data_g1,mapnum=mapnum,bmap=bmap)
  end

  '4pi_run05_sep2':begin
    mvn_sep_read_mult_sim_files,simstat_e,pathnames='g4work/davinMaven/results/4PI/run05/mvn_sep2_e-_atten0_seed??_.dat',type=-1,data_e
    mvn_sep_read_mult_sim_files,simstat_p,pathnames='g4work/davinMaven/results/4PI/run05/mvn_sep2_proton_atten0_seed??_.dat',type=+1,data_p
    ;mvn_sep_read_mult_sim_files,simstat_g,pathnames='g4work/davinMaven/results/4PI/run05/mvn_sep2_gamma_atten0_seed??_.dat',type=0,data_g
    str_element,/add,simstat_e,'desc','4Pi'
    str_element,/add,simstat_p,'desc','4Pi'
    resp_e0 = mvn_sep_inst_response(simstat_e,data_e,mapnum=mapnum,bmap=bmap)
    resp_p0 = mvn_sep_inst_response(simstat_p,data_p,mapnum=mapnum,bmap=bmap)
    ;resp_g0 = mvn_sep_inst_response(simstat_g,data_g,mapnum=mapnum,bmap=bmap)
    mvn_sep_read_mult_sim_files,simstat_e1,data_e1,pathnames='g4work/davinMaven/results/4PI/run05/mvn_sep2_e-_atten1_seed??_.dat',type=-1
    mvn_sep_read_mult_sim_files,simstat_p1,data_p1,pathnames='g4work/davinMaven/results/4PI/run05/mvn_sep2_proton_atten1_seed??_.dat',type=+1
    str_element,/add,simstat_e1,'desc','4Pi'
    str_element,/add,simstat_p1,'desc','4Pi'
    ;mvn_sep_read_mult_sim_files,simstat_g1,data_g1 ,pathnames='g4work/davinMaven/results/4PI/run05/mvn_sep2_gamma_atten1_seed??_.dat',type=0
    resp_e1 = mvn_sep_inst_response(simstat_e1,data_e1,mapnum=mapnum,bmap=bmap)
    resp_p1 = mvn_sep_inst_response(simstat_p1,data_p1,mapnum=mapnum,bmap=bmap)
    ;resp_g1 = mvn_sep_inst_response(simstat_g1,data_g1,mapnum=mapnum,bmap=bmap)
  end


endcase

filename = testrun+'_response-map-'+strtrim(mapnum,2)+'.sav'
save,file=filename,resp_e0,resp_p0,resp_g0,resp_e1,resp_p1,resp_g1,resp_O0,bmap,mapnum


undefine,f
;dir =[1,0,0]
;xdir = -1
str_element,/add,f,'xdir',xdir
str_element,/add,f,'erange',erange
str_element,/add,f,'dir',dir
str_element,/add,f,'angle',angle
str_element,/add,f,'ypos',ypos
str_element,/add,f,'f_det',f_det
str_element,/add,f,'a_side',a_side
str_element,/add,f,'det_af',det_af
str_element,/add,f,'col_af',col_af
str_element,/add,f,'derange',derange
str_element,/add,f,'detname',detname
str_element,/add,f,'impact_bmin',impact_bmin
str_element,/add,f,'impact_pos',impact_pos
printdat,f,output=s,/val
fdesc = strjoin(strcompress(s,/remove_all),', ')
dprint,fdesc ;,/val

if 1 then begin
  mvn_sep_inst_bin_response,simstat_e,data_e,mapnum=mapnum,noise_level=noise_level
  mvn_sep_inst_bin_response,simstat_p,data_p,mapnum=mapnum,noise_level=noise_level
  mvn_sep_inst_bin_response,simstat_g,data_g,mapnum=mapnum,noise_level=noise_level
  ;printdat,data_p1,simstat_p

  ;ok = mvn_sep_response_data_filter(simstat_e,data_e1,_extra=f,filter=f)
  win=0
  mvn_sep_response_plots,simstat_e,data_e,window=win,filter=f
  mvn_sep_response_plots,simstat_p,data_p,window=win,filter=f
  mvn_sep_response_plots,simstat_p1,data_p1,window=win,filter=f
  mvn_sep_response_plots,simstat_g,data_g,window=win,filter=f

endif

if 0 then begin
  resp=resp_p0
  bmap =resp.bmap
  omega=0
  over=0
  bins0 = where(strmatch(bmap.name,'B-O'))
  ;mvn_sep_response_each_bin,resp_p0,bins=bins0,window=20,ylog=1,omega=omega
  mvn_sep_response_each_bin_GF,resp_p0,bins=bins0,window=20,ylog=1,omega=omega,over=over
endif


if testrun eq 'oxygen_dir' then begin
  mvn_sep_inst_bin_response,simstat_Ox,data_Ox,mapnum=mapnum,noise_level=noise_level
  win=0
  mvn_sep_response_plots,simstat_Ox,data_Ox,window=win,filter=f

endif



if 0 then begin
  erange=[0,0]
  w = where(data_e1.einc lt 1e4 and data_e1.fto[0] eq 7)
  plot,[1,1],/nodata,xrange=[1,1e4],yrange=[1,1e4],/xlog,/ylog
  ;oplot,data_e1[w].edep[0,0],data_e1[w].edep[1,0],psym=3
  color = bytescale(data_e1[w].einc,range=[500,1e4],/log)
  plots,data_e1[w].edep[2,0],data_e1[w].edep[1,0],psym=4,color=color,symsize=.3
endif


end