/*
     LIPS
     shuttlerless photon transfer
       gain xy file binning soft
       
          gain_bin_v11.c : Ver 1.1
	  
          2002.9.24  H.Akitaya
	  2003.10.13 H.Akitaya

  input file : column#1 pixel(No.)
                     #2 count(ADU)
		     #3 Deviation(ADU^2)
		     #4 Deviation - readout noise^2 (ADU^2)
                     #5 gain = (#2)/(#4) (ADU/e-)

  output file : column#1: x1
                column#2: x2
                column#3: (x1+x2)/2
                column#4: Number of Data Points
                column#5: gain 
                column#6: stddev
                column#7: stddev / sqrt(#4)
		column#8: rejected data points
		
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#define VERSION "1.1 [2003.10.14 H.Akitaya]"

#define STEP 25
#define UPPER_REJECT_LIM 1.0  /* Allowed Maximum Gain Value */
#define NSREJ_UPPER 2.5       /* sigma rejection - upper limit */
#define NSREJ_LOWER 2.5       /* sigma rejection - lower limit */

main(int argc, char *argv[])
{
  FILE *fin, *fout;
  int scale, j, i, n_data, step, c_from, c_to;
  float px, rawc, dev_all, dev_c, gain, count;
  float mean_gain, tot_dev, stddev;
  float tmp_count[20000], tmp_gain[20000];
  
  int n_data_fin;
  float mean_gain_fin, stddev_fin, gain_tot,  tot_dev_fin;
  float error_sig;

/*  step = STEP;*/

  /* read parameters */
  
  if(argc != 6){
    printf("gain_bin ver.");
    printf( VERSION );
    printf( "\nusage : gain_bin (c_from) (c_to) (step) (data_fn) (output_fn)\n" );
   exit(1) ;
  }
  c_from = atoi( argv[1] );
  c_to   = atoi( argv[2] );
  step   = atoi( argv[3] );

  /* file open */
  
  if( (fin = fopen(argv[4], "r")) == NULL){
    printf("cannot open %s\n", argv[1]);
    exit(1);
  }
  fclose( fin );
  
  if( (fout = fopen(argv[5], "w")) == NULL){
    printf("cannot open %s\n", argv[2]);
    exit(1);
  }
  for(i=c_from; i<=c_to; i+=step ){

    

    /* store data into the array */
    
    mean_gain = 0.0;
    tot_dev = 0.0;
    n_data=0;         /* number of Data Points */
    
    fin = fopen(argv[4],"r");
    
    while( fscanf( fin, "%e%e%e%e%e", &px, &count, &dev_all, &dev_c, &gain) != EOF){
      if ( count > i && count <= (i+step) && gain < UPPER_REJECT_LIM){
	tmp_count[n_data] = count;
	tmp_gain[n_data] = gain;
	n_data++;
      }
      
    }
    fclose( fin );

    
    if( n_data > 2){

      /* calculate the mean */
	 
      for( j=1; j <= n_data; j++){
	mean_gain += ( tmp_gain[j-1] / n_data );
      }
      
      /* calculate the deviation and the std.dev */
	 
      for( j=1; j <= n_data; j++){
	tot_dev +=
	  ( (tmp_gain[j-1] - mean_gain) * (tmp_gain[j-1] - mean_gain) / n_data );
      }
      stddev = sqrt ( tot_dev );

      gain_tot = 0.0;
      tot_dev_fin = 0.0;
      n_data_fin = 0;

      /* n-sigma rejection */
      
      for( j=1; j <= n_data; j++) {
	if( (tmp_gain[j-1] <= (mean_gain + NSREJ_UPPER * stddev) ) &&
	   (tmp_gain[j-1] >= (mean_gain - NSREJ_LOWER * stddev ) ) ) {
	  gain_tot += tmp_gain[j-1];
	  n_data_fin++;
	}
      }
      mean_gain_fin = gain_tot / n_data_fin;
      for( j=1; j <= n_data; j++){
	if( (tmp_gain[j-1] <= (mean_gain + NSREJ_UPPER * stddev) ) &&
	   (tmp_gain[j-1] >= (mean_gain - NSREJ_LOWER * stddev ) ) ) {
	  tot_dev_fin +=
	    ( (tmp_gain[j-1] - mean_gain_fin) * (tmp_gain[j-1] - mean_gain_fin) / n_data_fin );
	}
      }
      stddev_fin = sqrt ( tot_dev_fin );

      if ( n_data_fin > 1 ){
	error_sig = stddev_fin / ( sqrt( n_data_fin - 1 ) );
      } else {
	error_sig = 0.0;
      }

      /* print out results into STDOUT and the file */
      
      printf( "%i %i %f %i %f %f %f %i\n",
	     i, i+step, (i+step/2.0), n_data_fin, mean_gain_fin,
	     stddev_fin , error_sig, ( n_data - n_data_fin ) );
      fprintf( fout, "%i %i %f %i %f %f %f %i\n",
	      i, i+step, (i+step/2.0), n_data_fin, mean_gain_fin,
	     stddev_fin , error_sig, ( n_data - n_data_fin ) );
    }
    
  }
  fclose( fout );
}