#!/usr/bin/perl
# 
# Use the chisq method to determine the probability that 
#    distribution of real streaks is same as distribution
#    of simulated streaks.
#
# MWR 7/1/2009

$debug = 1;

$start_bin = 25;
$end_bin = 35;
$d_bin = 5;

# this has the distribution of real streaks
$real_file = "./real_streaks_dist.dat";
# and this has the distribution of simulated streaks
$sim_file = "./sim_streaks_dist.dat";


# read in data 
$num_real = 0;
open(INFILE, "$real_file") || die("can't open file $real_file");
while (<INFILE>) {
  $line = $_;
  chomp($line);
  if ($line =~ /^#/) {
    next;
  }

  @words = split(/\s+/, " " . $line);
  $real_length_array[$num_real] = $words[1];
  $real_count_array[$num_real] = $words[2];
  if ($debug > 0) {
    printf " real elem %3d has length %6d count %10d \n",
             $num_real, $real_length_array[$num_real], $real_count_array[$num_real];
  }
    
  $num_real++;

}

$num_sim = 0;
open(INFILE, "$sim_file") || die("can't open file $sim_file");
while (<INFILE>) {
  $line = $_;
  chomp($line);
  if ($line =~ /^#/) {
    next;
  }

  @words = split(/\s+/, " " . $line);
  $sim_length_array[$num_sim] = $words[1];
  $sim_count_array[$num_sim] = $words[2];
  if ($debug > 0) {
    printf " sim elem %3d has length %6d count %10d \n",
             $num_sim, $sim_length_array[$num_sim], $sim_count_array[$num_sim];
  }
    
  $num_sim++;

}


# now do the chi-sq calculation
$chisq = 0;
for ($length = $start_bin; $length <= $end_bin; $length += $d_bin) {
  if ($debug > 0) {
    printf " now length %5d \n", $length;
  }


  # place the number of real streaks of this length in "p"
  $p = -1;
  for ($i_real = 0; $i_real < $num_real; $i_real++) {
    if ($real_length_array[$i_real] == $length) {
      $p = $real_count_array[$i_real];
      if ($debug > 0) {
        printf " real length %3d  count %10d \n", $length, $p;
      }
      last;
    }
  }
  if ($p == -1) {
    printf " failed to find real array element with length %5d \n", $length;
    exit(1);
  }

  # place the number of sim streaks of this length in "q"
  $q = -1;
  for ($i_sim = 0; $i_sim < $num_sim; $i_sim++) {
    if ($sim_length_array[$i_sim] == $length) {
      $q = $sim_count_array[$i_sim];
      if ($debug > 0) {
        printf " sim  length %3d  count %10d \n", $length, $q;
      }
      last;
    }
  }
  if ($q == -1) {
    printf " failed to find sim  array element with length %5d \n", $length;
    exit(1);
  }

  # normalize the simulated counts -- divide by 1000
  $q /= 1000;
  if ($debug > 0) {
    printf "  p %10.1f q %10.1f \n", $p, $q;
  }

  # now compute the contribution of this bin to the chisq
  $top = ($p - $q)*($p - $q);
  $bot = $p + $q;
  $this = $top / $bot;
  $chisq += $this;
  if ($debug > 0) {
    printf "   top %8.3f bot %8.3f  this %8.3f  chisq now %12.5e \n",
         $top, $bot, $this, $chisq;
  }





     


}



exit 0;
