; use: 
function old_mean_hist, obin2, mh
  dim = size(mh, /dimensions)
  me = fltarr(dim[0])
  for i=0,dim[0]-1 do begin
     me[i] = total(obin2*mh[i,*])
     if total(mh[i,*]) ne 0 then me[i,*] /= total(mh[i,*])  
  endfor
return, me
end

function mean_hist, x,y, min, max, nbins, weight=z,meanx=meanx
; set weight to one if it was not given
if min gt max then return, 0
if N_elements(z) ne N_elements(x) then begin 
   z = x
   z[*] = 1.
endif
h = histogram(x,nbins=nbins,MIN=min,MAX=max, reverse_indices=r)
my    = fltarr(nbins-1)
meanx = my 
for i=1L,nbins-1 do begin
   if r[i-1] ne r[i] then begin
      b=r[r[i-1]:(r[i]-1)]
      meanx[i-1] = mean(x[b])
      my[i-1] = total(y[b]*z[b])/total(z[b])
   endif
endfor
ind = where(h ne 0, count)
if count gt 0 then meanx = meanx[ind]
if count gt 0 then my = my[ind]

return, my
end

pro construct_histogram, data, ALL=all
@common_blocks.inc
  dims = LONG(total(all_grid_info[0].end_index gt 0)) ; number of dimensions
  if N_elements(all) ne 1 then all=0
  fields = [hi.field1, hi.field2, hi.weight]
  notnone = where(fields ne 'none' and fields ne '')
  if notnone[0] ge 0 then begin
     validfields = fields[notnone]
     Nnotnone = N_ELEMENTS(notnone)
  endif else begin
     print, 'no valid fields'
     return
  endelse
  isnone = fields eq 'none'
  Nnone = TOTAL(isnone) 
  isp = fltarr(Nnotnone)
  for i=0,Nnotnone-1 do isp[i] = FIX(is_particle_field(validfields[i]))
  
  if N_elements(uniq(isp[sort(isp)])) gt 1 then begin
     print, 'you mixed particle and grid data requesting this histogram... '
     print, isp
     return
  endif

  isparticles = isp[0]

  unitfactor = [get_unit(fields[0], ulabel=u1),get_unit(fields[1], ulabel=u2),get_unit(fields[2], ulabel=u3)]
  unitlabel = [u1,u2,u3]

  if isparticles then begin
     pfields = INTARR(N_elements(validfields))
     for i=0,N_elements(validfields)-1 do $
        pfields[i] = (where(*parti.names eq validfields[i]))[0]
     read_particle_data, fields=pfields
  endif

  
  titlestring1=validfields[0]+' ['+unitlabel[0]+']'
  if Nnotnone gt 1 then titlestring2=validfields[1]+' ['+unitlabel[1]+']' else titlestring2='dN '
;  titlestring1= (shouldbe_logged(validfields[0])) ? 'log '+titlestring1 : titlestring1
;  if Nnotnone gt 1 then titlestring2= (shouldbe_logged(validfields[1])) ? 'log '+titlestring2 : titlestring2


;make histogram
  if isparticles then begin
     typeind = where(*parti.names eq 'particle_type')
     ind = lindgen(N_elements((*parti.d[pfields[0]])))
     if parti.type gt 0 then begin
        th = histogram((*parti.d[typeind[0]])[ind])
        print, N_elements(th), ' different types of particles found.'
        if typeind[0] gt -1 then begin
           nind = where((*parti.d[typeind[0]])[ind] eq parti.type)
           if nind[0] gt -1 then ind = ind[nind] else print, 'no matching particle types found. Using all instead.'
        end else print, 'no field with name "particle_type" read'
     endif
     x = reform((*parti.d[pfields[0]])[ind])*unitfactor[0]
     if Nnotnone gt 1 then y = double(reform((*parti.d[pfields[1]])[ind])*unitfactor[1])
     if Nnotnone gt 2 then z = double(reform((*parti.d[pfields[2]])[ind]))
  endif else begin ; field data:
     if all then oj->r3D,data_fields=validfields[0:Nnotnone-1] $
     else oj->r3D, center=center, cube=[slice_size[2]-slice_size[0], cube_dim], $
                   data_fields=validfields[0:Nnotnone-1] 
     x = double(reform((oj->data3d(validfields[0])))*unitfactor[0])
     if Nnotnone gt 1 then y = double(reform((oj->data3d(validfields[1])))*unitfactor[1])
     if Nnotnone gt 2 then z = double(reform((oj->data3d(validfields[2]))))
  endelse

  min1 = 0 & min2 =0 & max1=1 & max2=1
  min1 = (min(x, max=max1)) 
  max1 = max1 
  if hi.max1 ne  1e30 then max1 = hi.max1
  if hi.min1 ne -1e30 then min1 = hi.min1

  if shouldbe_logged(validfields[0]) then begin
     min1 = (min1 gt 0.) ? alog10(min1) : 0
     max1 = (max1 gt 0.) ? alog10(max1) : 0
     x = alog10(temporary(x))
     x = temporary(x) < max1 > min1
  endif
  if Nnotnone gt 1 then begin
     if hi.min2 eq -1e30 then min2 = min(y) else min2 = hi.min2
     if hi.max2 eq  1e30 then max2 = max(y) else max2 = hi.max2
     if ((notnone[1] eq 1) and shouldbe_logged(validfields[1])) then begin
        min2 = (min2 gt 0.) ? alog10(min2) : 0
        max2 = (max2 gt 0.) ? alog10(max2) : 0
        y = alog10(temporary(y))
        y = temporary(y) < max2 > min2
     endif
  endif
  
  obin1 = [min1,max1]
  obin2 = [min2,max2]

print ,'show hist', nnotnone
print, 'ranges:'
print, validfields
print , min1, max1
print,  min2, max2

  case Nnotnone of 
     1: BEGIN
        mh = hist1d(x, $
                        min=min1, max=max1, $
                        binsize=(max1-min1)/(hi.nbins1-1), obin=obin1,DENSITY=dens)
        if (hi.cumoption eq 1) then mh = total(mh, /cumulative)
        if (hi.cumoption eq 2) then mh = reverse(total(reverse(mh), /cumulative))
; for cumulative ones normalize so largest value is 1.
        if ((not useunits) and (hi.cumoption gt 0)) then mh = 1.D*mh/max(mh) 

        obin1 = obin1- 0.5*(max1-min1)/N_elements(obin1)
        if shouldbe_logged(validfields[0]) then begin
           obin1=10.^obin1
           min1 = 10.D^min1
           max1 = 10.D^max1
        endif
        if hi.min2 eq -1e30 then min2 = min(mh) else min2 = hi.min2
        if hi.max2 eq  1e30 then max2 = max(mh) else max2 = hi.max2

;        print ,'show hist for one', nnotnone, obin1, mh
        plot, obin1, mh, xlog=shouldbe_logged(validfields[0]), $
              ylog=hi.logoption, psym=10, $
              xtitle=titlestring1, ytitle='dN', $
              xrange=[min1,max1], $
              xstyle=1, ystyle=1, yrange=[min2,max2], $
              charsize=2, thick=2, $
              /noclip, pos=[0.15,0.15,0.98,0.98], $
              ticklen=.02


     END
     2: BEGIN
        if notnone[1] eq 2 then begin
           mh = hist1d(x, y, $
                       min=min1, max=max1, $
                       binsize=(max1-min1)/hi.nbins1, obin=obin1,DENSITY=dens)
           if (hi.cumoption eq 1) then mh = total(mh, /cumulative)
           if (hi.cumoption eq 2) then mh = reverse(total(reverse(mh), /cumulative))
; for cumulative ones normalize so largest value is 1.
           if ((not useunits) and (hi.cumoption gt 0)) then mh /= max(mh) 

           if shouldbe_logged(validfields[0]) then obin1=10.^obin1
           if hi.min2 eq -1e30 then min2 = min(mh) else min2 = hi.min2
           if hi.max2 eq  1e30 then max2 = max(mh) else max2 = hi.max2
           
;           print ,'show hist for two', nnotnone, obin1, mh
           plot, obin1, (mh) , xlog=shouldbe_logged(validfields[0]), $
                 ylog=hi.logoption, psym=10, $
                 xtitle=titlestring1, ytitle='dN', $
                 xrange=[min(obin1),max(obin1)], $
                 yrange=[min2,max2], $
                 xstyle=1, ystyle=1, $
                 charsize=2, thick=2
           meanx = obin1
           meanh = mh
        endif else begin
           if ((dims eq 1) or (hi.cumoption eq 3)) then begin ; requested plot instead of histogram
              obin1 = x
              if shouldbe_logged(validfields[0]) then begin
                 obin1=10.^obin1
                 min1 = 10.^min1
                 max1 = 10.^max1                 
              endif
              mh = y
              if shouldbe_logged(validfields[1]) then begin
                 mh=10.^mh
                 min2 = 10.^min2
                 max2 = 10.^max2
              endif
              plot, obin1 > min1, mh , xlog=shouldbe_logged(validfields[0]), $
                    ylog=shouldbe_logged(validfields[1]), psym=4, $
                    xtitle=titlestring1, ytitle=titlestring2, $
                    xrange=[min1,max1], $
                    yrange=[min2,max2], $
                    xstyle=1, ystyle=1, $
                    charsize=2, thick=1,SYMSIZE=0.5
           endif else begin ; otherwise make a histogram
              mh = hist2d(x, $
                          y, $
                          min1=min1, max1=max1, $
                          min2=min2, max2=max2, $
                          binsize1=(max1-min1)/(hi.nbins1), $
                          binsize2=(max2-min2)/(hi.nbins2), $
                          obin1=obin1,obin2=obin2,DENSITY=dens)
              meanh = mean_hist(x, y, min1,max1,hi.nbins1, meanx=meanx)
           endelse
           END
     END
     3: BEGIN
        mh = hist2d(x, $
                    y, $
                    z, $
                    min1=min1, max1=max1, $
                    min2=min2, max2=max2, $
                    binsize1=(max1-min1)/(hi.nbins1), $
                    binsize2=(max2-min2)/(hi.nbins2), $
                    obin1=obin1,obin2=obin2,DENSITY=dens)
        meanh = mean_hist(x, y, min1,max1,hi.nbins1, weight=z,meanx=meanx)
     END
     ELSE: print, 'somthings wrong in show_histogram ..'
 end

  mh = double(mh)
  print, 'size(mh):', size(mh)
  tmh = total(mh,/double)
  print, 'used ', N_elements(x), ' points'
  print, 'total(histogram): ', tmh
  
  if (Nnotnone ge 2 and dims gt 1) then begin
     if hi.option eq 0 then begin  ; convert to fraction of total per bin 
        if not useunits then begin 
           mh = mh/(tmh)        ; not using units make it show a fraction
        endif
     endif else begin           ; average value rather than fraction of total
        ind = where(dens ne 0, count)
; if average is asked but no field don't average:
        if hi.weight ne 'none' and count gt 0 then  mh[ind] /= DOUBLE(dens[ind])
     endelse
  endif
  
  print, 'total(histogram): ', total(mh, /double)
  if N_elements(z) gt 0 then print, 'total(weight):', total(z, /double)
  
  
  data = mh
  mid = min(data, max=maxdata)
  ind = where(data gt mid)
  if ind[0] ge 0 then sec = min(data[ind]) else sec = 1e-40
  print, 'min, secmin and max in histogram:', mid, sec, maxdata
  data = shouldbe_logged(hi.weight) ? alog10(data > sec) : data > sec
  
  hi.titlestring1 = titlestring1
  hi.titlestring2 = titlestring2
  hi.xrange = [min(obin1),max(obin1)]
  hi.yrange  =[min(obin2),max(obin2)]

  x_stuetz = x
  if Nnotnone gt 1 then y_stuetz = y
  if Nnotnone gt 2 then z_stuetz = z

  if shouldbe_logged(validfields[0]) then x_stuetz = 10.^x
  if (Nnotnone gt 1) then if (notnone[1] eq 1) then  if shouldbe_logged(validfields[1]) then y_stuetz = 10.^y
  if Nnotnone gt 2 then  if shouldbe_logged(validfields[2]) then z_stuetz = 10.^z
  
  return
end