#!/usr/bin/perl
#
#  histo2.pl   perl script for histogram binning 
#  (c) Hans-Joachim Drescher    hajo@free.fr  
#  licence: GPL Version 2
#  drescher at th dot physik dot uni dash frankfurt dot de 
#
$opt_c=1;
$opt_w=0;
$opt_skip='^\s*$';
$opt_only='*';
use Getopt::Long;
use POSIX qw(ceil floor);
GetOptions(qw(dx:f f log sc:f c:i w:i h norm:i max nd:f x0:f s:f skip:s cond:s only:s));
$opt_c--;
$opt_w--;
if ($opt_h) {
print '
   make fast histograms with this perl script
   it takes input from files or standard input. 

   examples:

   bin the second column of some file with bin-width 2.5:
           histo2.pl  -c 2 -dx 2.5 file.txt

   bin the 3rd with weight in the 4th. logaritmic and 10 bins per decade
   and normalize by bin-width and number of entries:
           histo2.pl -log -nd 10 -c 3 -w 4 -norm 11 file.txt


   histo2.pl  -c col_x -w col_weight -dx dx -log -nd i -norm inorm -x0 X0
              -f 
	      
         -c col_x       column of x values [ 1 ]
         -w col_weight  column of weight [weight is 1 if not specified]
         -dx dx         bin width (is factor in log-case)  [ 10 ]
         -nd  i         computes dx as number of bins per decade in log case
         -log           log-binning
         -x0            bin-middle of any bin [100]
         -f             print empty bins 
         -max           takes naximum in each bin

         -norm ijk     normalization [0]
            k=0   no normalization
            k=1   normalize by number of entries
            k=4   normalize by number of entries in bin
            j=1   by bin-width
            i     not yet implemented 
         -s  f     scale y by 1/factor
         -sc f     scale x by 1/factor
';
exit;
}


$x0=100;
$x0 = $opt_x0 if defined($opt_x0);
$dx=10;
$dx=$opt_dx if defined($opt_dx);
$dx=10.0**(1.0/$opt_nd) if defined($opt_nd);

if( $opt_log) {
    $lx0=log($x0);
    $ldx=log($dx);
}

$n=0;
while (<>){
    next if /^\s*\#/;
    next if /$opt_skip/;
    (@A)=split;
    if ( defined($opt_cond)) {
	($u=$opt_cond) =~ s/\$([0-9]+)/\$A\[-1+\1\]/g;
	# print "here",eval($u),"--->$u<---$opt_cond\n";
	if( ! eval($u) ) {next;}
    }
    $n++;
    if( $opt_log ){ # log case
	if( $dx <= 0 ) {
	    $a=$A[$opt_c];   
	} elsif ( $A[$opt_c] > 0.0 )  {
	    $a=exp(floor((log($A[$opt_c])-$lx0)/$ldx+0.5)*$ldx+$lx0);
	}
	#print STDERR "log:",$a,"  $A[$opt_c]\n";
    }else{    # lin case 
	if( $dx <= 0 ) {            #
	    $a=$A[$opt_c];
	} else {
	    $a=$x0+floor(($A[$opt_c]-$x0)/$dx+0.5)*$dx;
	}
    }
    $a *= $opt_sc if defined($opt_sc); 
    if ( $opt_w > 0 ) {
	if( $opt_max ) {
	    $bin{$a} = $A[$opt_w] if  $bin{$a} < $A[$opt_w] ;
	} else {
	    $bin{$a} += $A[$opt_w];     # 
	    $bin2{$a} += $A[$opt_w]*$A[$opt_w];     # 
	}
    } else {
	$bin{$a} += 1;	            #
    }
    $nbin{$a} += 1;
}

$s=0.0;$s2=0.0;$w=0.0;$ny=0.0;
foreach $x ( sort {$a<=>$b} keys(%bin)) {
    if( $opt_log ){
	$deltax=$x*(sqrt($dx)-1.0/sqrt($dx));
    } else {
	$deltax=$dx;
    }
    $y=$bin{$x};
    $y /= $deltax   if $opt_norm/10 % 10 == 1;
    $y /= $n        if $opt_norm    % 10 == 1;
    $y /= $nbin{$x} if $opt_norm    % 10 == 4;
    $y /= $opt_s    if defined($opt_s);
 
    $dy=sqrt($bin2{$x});
    $dy /= $deltax   if $opt_norm/10 % 10 == 1;
    $dy /= $n        if $opt_norm    % 10 == 1;
    $dy /= $nbin{$x} if $opt_norm    % 10 == 4;
    $dy /= $opt_s    if defined($opt_s);

    if( defined($xlast) && $opt_f ) {
	if( $opt_log ){
	    $i=int(log($x/$xlast)/$ldx+0.5);
	    for ($ii=1;$ii<$i;$ii++){
		print $xlast*$dx**$ii," 0 0\n";
	    }
	} else {
	    $i=int(($x-$xlast)/$dx+0.5);
	    for ($ii=1;$ii<$i;$ii++){
		print $xlast+$ii*$dx," 0 0\n";
	    }
	}
    }
    print "$x ",$y," ",$dy,"  ",$nbin{$x},"\n";
    $s+=$x*$y;
    $s2+=$x*$x*$y;
    $w+=$y;
    $xlast=$x;
}
if ( $n != 0 ) {
    printf STDERR "%s %e %s %e %s %e %s %e\n","total:",$n,"mean:",$s/$w," sig ",sqrt($s2/$w-($s/$w)**2),"  em ",sqrt($s2/$w-($s/$w)**2)/sqrt(1.0*$n),"\n";
}

print "\n\n";

sub cond {
    $ae=eval("$opt_cond");
    print STDERR "$_ $opt_cond -->$ae<--\n";
    return $ae;
}