spectrum.f

      integer nbin_t,bindim=1000
      real*8 minbin_t,maxbin_t,dxbin_t,x_axis_observable,t_histo(1:bindim),
     &       t_dnsgrd(1:bindim),t_errgrd(1:bindim)


      minbin_t=1
      maxbin_t=101
      nbin_t=50
      dxbin_t=(maxbin_t-minbin_t)/nbin_t

        do-loop over all particles in all events:

            index=int((x_axis_observable - minbin_t)/dxbin_t+1.5)
            if(index.lt.1.or.index.gt.nbin_t) goto 102
            t_histo(index)=t_histo(index)+1.0    % here you could increment by weight instead of by 1
            t_dnsgrd(index)=t_dnsgrd(index)+1
            t_errgrd(index)=t_errgrd(index)+1.0*1.0  % here the increment would be weight^2
 102        continue   % end of do-loop


c now normalization:
      do 31 i=1,nbin_t
         t_histo(i)=t_histo(i)/(dxbin_t)   % divide by bin-width to make independent of bin-width
         dns=t_dnsgrd(i)+1.d-8    % make sure for error bar you do not divide by zero
c     for dN/d_something
         t_errgrd(i)=t_histo(i)/sqrt(dns)  % error bar for this kind of histogram
 31   continue