#!/usr/bin/perl

# short program to find the moments of a histogram.  Stdin has two columns: x y
# Written by Joshua Spodek
# Copyright (c) 1998 Joshua Spodek
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 2
# of the License, or (at your option) any later version.
# 
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
# 
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

# set the defaults
$xcolumn=1;
$ycolumn=2;

# Read the command line arguments
while ($_=$ARGV[0]) {
    shift;
    /^-xmin/    && ($xmin=@ARGV[0],shift,next);
    /^-xmax/    && ($xmax=@ARGV[0],shift,next);
    /^-binsize/ && ($binsize=$ARGV[0],shift,next);
    /^-hew/     && ($HEW=1,next);
    /^-xcolumn/ && ($xcolumn=$ARGV[0],shift,next);
    /^-ycolumn/ && ($ycolumn=$ARGV[0],shift,next);
    print STDERR "unknown argument: $_\n";
}
if (defined($xmin) && defined($xmax)) {
    if ($xmin > $xmax) {die "lower limit must be less than upper limit\n"}
}

# Load in the data and calculate the number of counts and center
$i=0;
$num_counts = 0;
$xy_sum = 0;
while (<>) {
    next unless (/^[\d-\s]/);
    $x = (split(' '))[$xcolumn-1];
    $y = (split(' '))[$ycolumn-1];
    next if (defined($xmin) && ($x < $xmin));
    last if (defined($xmax) && ($x > $xmax));
    $num_counts += $y;
    $xy_sum     += $x * $y;
    $x[$i]=$x;
    $y[$i]=$y;
    $i++
}
$num_points=$i;
$binsize = $x[1] - $x[0];
if ($binsize == 0) {die "Error in histstats: binsize = 0!\n"}
if ($num_counts == 0) {die "Error in histstats: No counts!\n"}
$center = $xy_sum / (1.0*$num_counts);

# Now calculate the other moments of the distribution
$x2y_sum=0;
$x3y_sum=0;
$x4y_sum=0;
for ( $i=0; $i<$num_points; $i++ ) {
    $x2y_sum += $y[$i]*($x[$i] - $center)**2;
    $x3y_sum += $y[$i]*($x[$i] - $center)**3;
    $x4y_sum += $y[$i]*($x[$i] - $center)**4;
}

$second_moment = $x2y_sum / (1.0*$num_counts);
$third_moment  = $x3y_sum / (1.0*$num_counts);
$fourth_moment = $x4y_sum / (1.0*$num_counts);

$sigma    = sqrt($second_moment);
$skewness = $third_moment / $second_moment**(3/2);
$kurtosis = $fourth_moment / $second_moment**2;

# Now find the half energy width (really the half count width)
if (defined($HEW)) {
    for ( $i=1; $i<=$#x/2; $i++ ) {
	$counts = 0;
	for ( $j=1; $j<=$#x; $j++ ) {
	    if ( abs($x[$j] - $center) < $i*$binsize && $x[$j] ) {
		$counts+=$y[$j];
	    }
	}
	$HEW = 2*$i*$binsize;
	last if ($counts >= $num_counts/2);
    }
}

printf( "%6.4f %6.4f %6.4f %6.4f %6.4f\n",
       $center, $sigma, $num_counts, $binsize, $HEW )

#print STDERR "skewness = $skewness\n";
#print STDERR "kurtosis = $kurtosis\n";