#!/usr/bin/perl -w
use strict;
use Data::Dumper;
# use POSIX qw(ceil floor);
use GD;
use GD::Text;
use GD::Graph;
# note this is Andrew's version...
use GD::Graph::histogram;
use PDL;
use PDL::Fit::Gaussian;

my ($pixx,$pixy,$infile)=@ARGV;

# print "$pixx $pixy $infile\n";

local $| = 1;

# It's not necessary to define the event hash, but it helps maintain the code...
my %event = (
   numr		=> 0.0,
   dnom		=> 0.0,
   tsr		=> 0,		# time since rese
   pre		=> 0,		# baseline
   post		=> 0,
   ph0		=> 0,		# post - pre
   ph1		=> 0.0,		# offset-corrected
   E		=> 0.0,
);

our ($m,$b,$r, @res); 		# linear regression fit params

our @bin = ();
our @freq = ();
our @rawpeaks = ();
our @rawcnts = ();
my @fitpeaks = ();
my @fwhms = ();

my @dacsets = (5,27,49,71,93);	# assume equal spacing between external DAC settings

my @pixeldata = ();
# Event data organized by starting-cap as follows: scap[starting-cap].event[ii].numr ...
my @scap;
my @capcnts;

my @nev;
my $ii = 0;
my $jj = 0;
my $kk = 0;
for ($ii=0;$ii<32;$ii++) {
   for ($jj=0;$jj<32;$jj++) {
      $nev[$ii][$jj] = 0;
   }
}

# print "Reading input data...";
open INPUT, "<$infile" or die "ERROR - Can't open $infile\n";
my $sc = 0;
my $x = 0;
my $y = 0;
my @fields;
my $evtot = 0;

while( <INPUT> ) {
   chomp;
   @fields = split " ", $_;
   $x = $fields[3];
   $y = $fields[4];
   $pixeldata[$x][$y][$nev[$x][$y]] = $_;
   $nev[$x][$y]++;
   $evtot++;
}
close INPUT;
# print "Done\n";


#for ($ii=0;$ii<32;$ii++) {
#   for ($jj=0;$jj<32;$jj++) {
#      print "$nev[$ii][$jj] ";
#   }
#}

my $minx=0;
my $miny=0;
my $maxx=31;
my $maxy=31;
if ($pixx != -1) {
   $minx=$pixx;
   $miny=$pixy;
   $maxx=$pixx;
   $maxy=$pixy;
}

# For each pixel...
for ($ii=$minx;$ii<=$maxx;$ii++) {		# each X
   for ($jj=$miny;$jj<=$maxy;$jj++) {		# each Y
      @scap = ();
      @capcnts = (0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0);

      # For each event in pixel (ii,jj), sort by starting-cap, and pick out numr,dnom,pre,post
#      print "$nev[$ii][$jj]\n";
      for($kk=0;$kk<$nev[$ii][$jj];$kk++) {
         @fields = split " ",  $pixeldata[$ii][$jj][$kk];
#	 print $pixeldata[$ii][$jj][$kk];

         # Cut on time-of-rise denominator
         if ($fields[2]>20.0) {
            $sc = $fields[0];
            $scap[$sc]{event}[$capcnts[$sc]]{numr} = $fields[1];
            $scap[$sc]{event}[$capcnts[$sc]]{dnom} = $fields[2];
            $scap[$sc]{event}[$capcnts[$sc]]{pre} = $fields[5];
            $scap[$sc]{event}[$capcnts[$sc]]{post} = $fields[6];
            $scap[$sc]{event}[$capcnts[$sc]]{ph0} = $fields[6] - $fields[5];
            $scap[$sc]{event}[$capcnts[$sc]]{tsr} = $fields[7];
            $capcnts[$sc]++;
	 }
      }


      #for ($ll=0;$ll<$capcnts[1];$ll++) {
      #   print "$scap[1]{event}[$ll]{pre} $scap[1]{event}[$ll]{post} $scap[1]{event}[$ll]{numr} $scap[1]{event}[$ll]{dnom}\n";
      #}

      # Look for pixels with too few events per starting-cap
      my $stop = 0;
      my $cap = 0;
      for ($cap=0;$cap<16;$cap++) {
         if ($capcnts[$cap] < 50) {
	    $stop = 1;
	 }
      }
      if ($stop) {
	 printf(STDERR "%2d %2d ",$ii,$jj);
	 print STDERR "ERROR too few events: ";
	 for ($cap=0;$cap<16;$cap++) {
	    print STDERR "$capcnts[$cap] ";
         }
         print STDERR "\n";	   
         next;
      }

      my @capdata = ();
      my $ll = 0;
      for ($cap=0;$cap<16;$cap++) {
         for ($ll=0;$ll<$capcnts[$cap];$ll++) {
            $capdata[$cap][$ll] = $scap[$cap]{event}[$ll]{ph0};
         }
      }      

      # First-pass thru data for each starting-cap, checking for enough statistics


#      my $min = 0;
#      my $width = 1500;
#      my $max = int(9000/$width); # channel_max/width
#      for ($cap=0;$cap<16;$cap++) {

#         @bin = ();
#         @freq = ();
#         my_histo ($width, $min, $max, @{$capdata[$cap]});

#         $stop = 0;
#         for ($ll = 0; $ll < 5; $ll++) {
#	    if ($freq[$ll] < 5) {
#	       $stop = $ll + 1;
#	       last;
#	    }
#	 }
#	 if ($stop) {
##            for (my $i = $min; $i <= $max; $i++) {
##               print "$bin[$i]   $freq[$i]\n";
##            }
#
#	    printf(STDERR "%2d %2d ",$ii,$jj);
#	    printf(STDERR "WARNING - not enough data for cap %2d: %2d %2d %2d %2d %2d\n",
#	    	$cap, $freq[0], $freq[1], $freq[2], $freq[3], $freq[4]);
##	    last;
#         }
#      }

      # Second-pass, more accurate peak-find. If 5 peaks found, compute mean of the 5 peaks
      $stop = 0;
      $ll = 0;
      my $min = 0;
      my $width = 10;
      my $max = int(8000/$width); # channel_max/width
      my @rpmean = ();
      my @p_offset = ();
      my @offset = ();
      for ($cap=0;$cap<16;$cap++) {

         @bin = ();
         @freq = ();
         my_histo ($width, $min, $max, @{$capdata[$cap]});

         @rawpeaks = ();
         @rawcnts = ();
         findpeaks($width,-3,5,3,2,$min,$max);
      
	 my $tmp = 0;
	 if (($#rawpeaks == 5)&&(($rawpeaks[5]-$rawpeaks[4]) <100)&&($rawpeaks[5]>6500)) {
	    printf(STDERR "%2d %2d ",$ii,$jj);
	    print STDERR "WARNING - found 6 peaks for cap $cap: ";
            for ($ll = 0; $ll <= $#rawpeaks; $ll++) {
               printf(STDERR "%7.2f ",$rawpeaks[$ll]);
            }
            print STDERR " - ignoring 5th peak\n";
	    splice(@rawpeaks, 4, 1);
	 }
	 if (($#rawpeaks < 2 ) || ($#rawpeaks > 4)) {
	    $tmp = $#rawpeaks +1;
#            for (my $i = $min; $i <= $max; $i++) {
#               print "$bin[$i]   $freq[$i]\n";
#            }
	    printf(STDERR "%2d %2d ",$ii,$jj);
	    print STDERR "ERROR - found $tmp peaks for cap $cap: ";
            for ($ll = 0; $ll <= $#rawpeaks; $ll++) {
               printf(STDERR "%7.2f ",$rawpeaks[$ll]);
            }
            print STDERR "\n";
	    $stop = 1;
	    last;
         }
	 if ($#rawpeaks != 4 ) {
	    $tmp = $#rawpeaks +1;
#            for (my $i = $min; $i <= $max; $i++) {
#               print "$bin[$i]   $freq[$i]\n";
#            }
	    printf(STDERR "%2d %2d ",$ii,$jj);
	    print STDERR "WARNING - found $tmp peaks for cap $cap: ";
            for ($ll = 0; $ll <= $#rawpeaks; $ll++) {
               printf(STDERR "%7.2f ",$rawpeaks[$ll]);
            }
            print STDERR "\n";
         }
         for ($ll = 0; $ll < 5; $ll++) {
            $rpmean[$cap][$ll] = -1;
	 }
         for ($ll = $#rawpeaks; $ll >= 0; $ll--) {
	    if (($rawpeaks[$ll]>350)&&($rawpeaks[$ll]<450)&&($rpmean[$cap][0] == -1)) {
               $rpmean[$cap][0] = $rawpeaks[$ll];
	    }
	    if (($rawpeaks[$ll]>1700)&&($rawpeaks[$ll]<2300)&&($rpmean[$cap][1] == -1)) {
               $rpmean[$cap][1] = $rawpeaks[$ll];
	    }
	    if (($rawpeaks[$ll]>3300)&&($rawpeaks[$ll]<4100)&&($rpmean[$cap][2] == -1)) {
               $rpmean[$cap][2] = $rawpeaks[$ll];
	    }
	    if (($rawpeaks[$ll]>4800)&&($rawpeaks[$ll]<5800)&&($rpmean[$cap][3] == -1)) {
               $rpmean[$cap][3] = $rawpeaks[$ll];
	    }
	    if (($rawpeaks[$ll]>6300)&&($rawpeaks[$ll]<7600)&&($rpmean[$cap][4] == -1)) {
               $rpmean[$cap][4] = $rawpeaks[$ll];
	    }
	 }
      }

      if ($stop) {
         next;
      }

#      print Dumper(@rpmean);

      # Average the 16 peak-means
      my @rpmm = ();
      my $ccnt = 0;
      for ($ll = 0; $ll <5; $ll++) {
         $rpmm[$ll] = 0;
	 $ccnt = 0;
         for ($cap=0;$cap<16;$cap++) {
            if ( $rpmean[$cap][$ll] >0) {
	       $rpmm[$ll] += $rpmean[$cap][$ll];
	       $ccnt++;
	    }
         }
	 if ($ccnt==0) {
#           print STDERR Dumper(@rpmean);
	    printf(STDERR "%2d %2d ",$ii,$jj);
	    print STDERR "ERROR - found NO peak $ll for any cap\n";
	    $stop = 1;
	    last;
	 }
	 
         $rpmm[$ll] /= $ccnt;
         for ($cap=0;$cap<16;$cap++) {
            if ( $rpmean[$cap][$ll] >0) {
	       $p_offset[$cap][$ll] = $rpmean[$cap][$ll] - $rpmm[$ll];
	    } else {
	       $p_offset[$cap][$ll] = -100000;
	    }
	 }
      }

      if ($stop) {
         next;
      }
#      print Dumper(@rpmm);
#      print Dumper(@p_offset);

      # Compute the overall offset for each starting-cap
      for ($cap=0;$cap<16;$cap++) {
         $offset[$cap] = 0;
	 $ccnt = 0;
#	 printf(STDERR "%2d   ",$cap);
         for ($ll = 0; $ll <5; $ll++) {
	    if ($p_offset[$cap][$ll] > -99999) {
               $offset[$cap] += $p_offset[$cap][$ll];
	       $ccnt++;
	    }
#	    printf(STDERR "%4.1f ",$p_offset[$cap][$ll]);
	 }
	 $offset[$cap] = $offset[$cap]/$ccnt;   
#         printf(STDERR "%4.1f\n",$offset[$cap]);
      }

      # For each event for this pixel, compute the energy-estimate E 
      my @Elist = ();
      my @bline = ();
      my @tsr = ();
      for ($cap=0;$cap<16;$cap++) {
         for ($ll=0;$ll<$capcnts[$cap];$ll++) {
            $scap[$cap]{event}[$ll]{ph1} = $scap[$cap]{event}[$ll]{post} - $scap[$cap]{event}[$ll]{pre} - $offset[$cap];
            $scap[$cap]{event}[$ll]{E} = $scap[$cap]{event}[$ll]{ph1} *
            	(1.0 + (2.15E-3*$scap[$cap]{event}[$ll]{numr}/$scap[$cap]{event}[$ll]{dnom}));
            push @Elist, $scap[$cap]{event}[$ll]{E};
            push @bline, $scap[$cap]{event}[$ll]{pre};
            push @tsr, $scap[$cap]{event}[$ll]{tsr};
	 }
      }

#      for ($ll = 0; $ll <= $#Elist; $ll++) {
#         printf(STDOUT "%7.2f\n",$Elist[$ll]);
#      }

      # Histogram E-values for this pixel, and find the 5 peaks
      @bin = ();
      @freq = ();
      $min = 0;
      $width = 4;
      $max = int(8000/$width); # channel_max/width
         
      my_histo ($width, $min, $max, @Elist);
#      for (my $i = $min; $i <= $max; $i++) {
#         print "$i $bin[$i]   $freq[$i]\n";
#      }
      @rawpeaks = ();
      @rawcnts = ();
      my $ampl = int($#Elist/50);
#      print "ampl= $ampl\n";
      my $pthresh = -15;
      if ($ampl < 50) { $pthresh = -10;}
      findpeaks($width,$pthresh,$ampl,4,3,$min,$max);
      my $tmp = 0.0;

#      for ($ll = 0; $ll <= $#rawpeaks; $ll++) {
#         printf(STDOUT "%7.2f\n",$rawpeaks[$ll]);
#      }

      # Check for shoulder on low-side of highest peak
      if (($#rawpeaks == 5)&&(($rawpeaks[5]-$rawpeaks[4]) <100)&&($rawpeaks[5]>6500)&&($rawcnts[5]/$rawcnts[4] > 5 )) {
	 printf(STDERR "%2d %2d ",$ii,$jj);
	 print STDERR "WARNING - found 6 peaks in E-binning: ";
         for ($ll = 0; $ll <= $#rawpeaks; $ll++) {
            printf(STDERR "%7.2f ",$rawpeaks[$ll]);
         }
         print STDERR " - ignoring 5th peak\n";
	 splice(@rawpeaks, 4, 1);
      }

      # Bail-out for this pixel if npeaks !=5
      if ( $#rawpeaks !=4) {
	 $tmp = $#rawpeaks +1;
	 printf(STDERR "%2d %2d ",$ii,$jj);
	 print STDERR "ERROR: Found $tmp peaks in E-binning: ";
         for ($ll = 0; $ll <= $#rawpeaks; $ll++) {
            printf(STDERR "%7.2f ",$rawpeaks[$ll]);
         }
         print STDERR "\n";
	 next;
#         for (my $i = $min; $i <= $max; $i++) {
#            print "$bin[$i]   $freq[$i]\n";
#         }
      }

      # refine peak estimation with a gaussian fit

      @fitpeaks = ();
      @fwhms = ();
      for ($ll = 0; $ll <= $#rawpeaks; $ll++) {
         my $b_idx = int($rawpeaks[$ll]/$width);
         my @xdat = @bin[($b_idx-7)..($b_idx+7)];
         my @ydat = @freq[($b_idx-7)..($b_idx+7)];

#      for (my $i = 0; $i <= $#xdat; $i++) {
#         print "$xdat[$i]   $ydat[$i]\n";
#      }
      
         my $x_pdl = pdl(\@xdat);
         my $y_pdl = pdl(\@ydat);
         my ($cen, $pk, $fwhm, $back, $err, $fit) = fitgauss1d($x_pdl, $y_pdl);
#         print "gfit $ll: $cen,$pk,$fwhm\n";
	 $fitpeaks[$ll] = $cen + ($width/2.0);
	 $fwhms[$ll] = $fwhm;
	 if ($fwhm > 40) {
            printf(STDERR "%2d %2d ",$ii,$jj);
	    print STDERR "WARNING - FWHM $ll too wide: $fwhm\n";
	    $fwhms[$ll]=38.0;
#	    print STDERR Dumper(@bin[($b_idx-10)..($b_idx+10)]);
#	    print STDERR "\n";
#	    print STDERR Dumper(@freq[($b_idx-10)..($b_idx+10)]);
#	    print STDERR "\n";
	 }
      }

      # fit a straight line to the 5 peaks
      linfit(5,\@dacsets,\@fitpeaks);

      # print out all relevent results
      printf("%2d %2d ",$ii,$jj);
#      for ($ll = 0; $ll <5; $ll++) {
#         printf("%8.3f ",$rawpeaks[$ll]);
#      }
      for ($ll = 0; $ll <5; $ll++) {
         printf("%8.3f ",$fitpeaks[$ll]);
      }
      for ($ll = 0; $ll <5; $ll++) {
         printf("%8.3f ",$fwhms[$ll]);
      }

      printf("% 9.6e  % 9.6e  % 9.6e  ",$m,$b,$r);
      for ($ll = 0; $ll <5; $ll++) {      
         printf("% 6.4e  ",$res[$ll]/$fitpeaks[4]);
      }
      for ($cap=0;$cap<16;$cap++) {
         printf("% 4.2f  ",$offset[$cap]);
      }
      print "\n";      


      # fit a straight line to the baseline vs time-since-reset data
      linfit($#tsr +1,\@tsr,\@bline);
#      printf("% 9.6e  % 9.6e  % 9.6e  \n",$m,$b,$r);
      my @tsr2 = ();
      my @bline2 = ();
      for ($ll = 0; $ll <=$#tsr; $ll++) {      
         if ($res[$ll] > 200) {
#	    printf("%d % 6.4e\n",$ll,$res[$ll]);
	 } else {
            push @tsr2, $tsr[$ll];	    
            push @bline2, $bline[$ll];	    
	 }
      }
#      print "\n";      

      linfit($#tsr2 +1,\@tsr2,\@bline2);
#      printf("% 9.6e  % 9.6e  % 9.6e  \n",$m,$b,$r);
      for ($ll = 0; $ll <=$#tsr; $ll++) {      
         if ($res[$ll] > 200) {
#	    printf("%d % 6.4e\n",$ll,$res[$ll]);
	 }
      }
   }
}

sub findpeaks {
 my ($bin_width,$pthresh,$ampl,$llim,$ulim,$min,$max) = @_;
 my $odiff = 0;
 my $diff = 0;
 my $ii = 0;
 my $jj = 0;
 for ($ii = 6; $ii <= ($max - $min -4); $ii++) {
   my $diff = $freq[$ii] - $freq[$ii -1];
   if ( (($diff <= $pthresh) && ($odiff >= 0)) || (($diff <0) && ($freq[$ii -1] >=$ampl)) ) {
      my $peak = 0;
      my $tot = 0;
      for ($jj=$ii-$llim;$jj<=$ii+$ulim;$jj++) {
         $peak += ($freq[$jj] * $bin[$jj]);
	 $tot += $freq[$jj];
      }
#      if ($tot <15) {
#	 print "Sparse data: ";
#         for ($jj=$ii-$llim;$jj<=$ii+$ulim;$jj++) {
#	    print "$freq[$jj] ";
#         }
#	 print "\n";
#      }
      if ($tot > 10) {
         $peak = $bin_width/2.0 + $peak/$tot -0.0;
         if (! $rawpeaks[$#rawpeaks]) {
            push @rawpeaks, $peak;
            push @rawcnts, $tot;
         } else {
            if ($peak > (50.0 + $rawpeaks[$#rawpeaks])) {
               push @rawpeaks, $peak;
               push @rawcnts, $tot;
            }
         }
      }
   }
   $odiff = $diff;
 }
}

sub round {
    my($number) = shift;
    return int($number + .5 * ($number <=> 0));
}

sub my_histo {
  my ($bin_width, $min, $max, @list) = @_;

  # This calculates the frequencies for all available bins in the data set
  my %histogram;
  $histogram{ceil(($_ + 1) / $bin_width) -1}++ for @list;

#  my $max;
#  my $min;

  # Calculate min and max
#  while ( my ($key, $value) = each(%histogram) )
#  {
#    $max = $key if !defined($min) || $key > $max;
#    $min = $key if !defined($min) || $key < $min;
#  }


  for (my $i = $min; $i <= $max; $i++)
  {
    $bin[$i]       = sprintf("% 10d", ($i) * $bin_width);
    $freq[$i] = $histogram{$i} || 0;

#    $freq = "#" x $freq;
#    print $bin." ".$freq."\n";
  }

#  print "===============================\n\n";
#  print "    Width: ".$bin_width."\n";
#  print "    Range: ".$min."-".$max."\n\n";
}


# ---------------------------------------------------------------------------------------------------------------------------------
#   sqr()  -  Return the square of a number.
# ---------------------------------------------------------------------------------------------------------------------------------

sub sqr {
$_[0] * $_[0];
}


# ---------------------------------------------------------------------------------------------------------------------------------
#   linear regression analysis for a set of data given as (x,y) pairs.
# ---------------------------------------------------------------------------------------------------------------------------------
sub linfit {

my($n, $xref, $yref) = @_;

my @x = @{$xref};
my @y = @{$yref};

our ($m,$b,$r, @res);

my $sumx = 0;
my $sumx2 = 0;
my $sumxy = 0;
my $sumy = 0;
my $sumy2 = 0;
my $ii = 0;

for ($ii=0;$ii<$n;$ii++) {                                                          # loop for all data points
   $sumx  += $x[$ii];                                                               # compute sum of x
   $sumx2 += $x[$ii] * $x[$ii];                                                     # compute sum of x**2
   $sumxy += $x[$ii] * $y[$ii];                                                     # compute sum of x * y
   $sumy  += $y[$ii];                                                               # compute sum of y
   $sumy2 += $y[$ii] * $y[$ii];                                                     # compute sum of y**2
}


# ---------------------------------------------------------------------------------------------------------------------------------
#   Compute and print results.
# ---------------------------------------------------------------------------------------------------------------------------------

$m = ($n * $sumxy  -  $sumx * $sumy) / ($n * $sumx2 - sqr($sumx));                  # compute slope
$b = ($sumy * $sumx2  -  $sumx * $sumxy) / ($n * $sumx2  -  sqr($sumx));            # compute y-intercept
$r = ($sumxy - $sumx * $sumy / $n) /                                                # compute correlation coefficient
         sqrt(($sumx2 - sqr($sumx)/$n) * ($sumy2 - sqr($sumy)/$n));

#printf "\nSlope        m = %15.6e\n", $m;                                           # print results
#printf "y-intercept  b = %15.6e\n", $b;
#printf "Correlation  r = %15.6e\n", $r;

my @fy;
for ($ii=0;$ii<$n;$ii++) {                                                      
   $fy[$ii] = $m*$x[$ii] + $b;
#   $res[$ii] = ($y[$ii]-$fy[$ii])/$y[$n-1];
   $res[$ii] = ($y[$ii]-$fy[$ii]);
#   printf("%8.4f  %8.4f  %.4e\n",$y[$ii],$fy[$ii],($y[$ii]-$fy[$ii])/$y[$n-1]);
}

}