#!/usr/bin/env ruby
#
# subtract.rb
#
# Copyright Jun Makino All rights reserved.
$LOAD_PATH.unshift(File.dirname($0))
$LOAD_PATH.push(ENV["ACSROOT"]+"/src/utilities") if ENV["ACSROOT"] != nil
require "clop.rb"

options_definition_string = <<-END

  Description: process TA100 spectrum data
  Long description:
    This program compares control data (background) and actual measured
    data
    (c) 2011- Jun Makino


  Short name: 		-l
  Long name:		--lower_limit_for_energy
  Value type:		float
  Default value:        0.3
  Variable name: 	emin
  Description:		Lower limit for calculating sum
  Long description:
    Lower limit for calculating the accumulated difference




  Short name: 		-u
  Long name:		--upper_limit_for_energy
  Value type:		float
  Default value:        1.5
  Variable name: 	emax
  Description:		Upper limit for calculating sum
  Long description:
    Upper limit for calculating the accumulated difference


  Short name: 		-d
  Long name:		--data_file
  Value type:		string
  Variable name: 	data_file
  Description:		Input data file
  Long description:
    The name of the input data file. The format is TA100U 
    spectrum data (text file)




  Short name: 		-b
  Long name:		--background_data_file
  Value type:		string
  Variable name: 	background_file
  Description:		Input background file
  Long description:
    The name of the input background data file. The format is TA100U 
    spectrum data (text file)





  Short name: 		-w
  Long name:		--sample_wait
  Value type:		float
  Default value:        100
  Variable name: 	sample_weight
  Description:		sample weight in gram
  Long description:
    The weight of the sample in unit of gram




  Short name: 		-k
  Long name:		--kalium_fraction
  Value type:		float
  Default value:        0
  Variable name: 	kalium_fraction
  Description:		kalium fraction in persentage
  Long description:
    kalium fraction in persentage





  Short name: 		-c
  Long name:		--ce_conversion_factor
  Value type:		float
  Default value:        110
  Variable name: 	ce_conversion_factor
  Description:		Caesium conversion factor
  Long description:
   Caesium conversion factor in unit of count/Bq





  Short name: 		-K
  Long name:		--kaliumn_conversion_factor
  Value type:		float
  Default value:        70
  Variable name: 	k_conversion_factor
  Description:		kalium  conversion factor
  Long description:
   Kalium conversion factor in unit of count/Bq







  END

$clop = parse_command_line(options_definition_string, false)

def fabs(x)
   y = x
  y = -y if y < 0
  y
end

def convert_to_bq_per_kg(count, sigma)
#  STDERR.print count, " ",sigma, "\n"
  kalium_count = $clop.kalium_fraction*$clop.sample_weight*$clop.k_conversion_factor
  net_count = count - kalium_count
  net_count_error = sigma * 3
  conversion_factor = 1000/$clop.sample_weight/$clop.ce_conversion_factor
  STDERR.print "\nCs amount = #{net_count*conversion_factor} +- "
  STDERR.print " #{net_count_error*conversion_factor} Bq/kg\n"
end

data_file=open($clop.data_file,"r") 
background_file=open($clop.background_file,"r") 
emin = $clop.emin
emax = $clop.emax
data_file.gets
background_file.gets
files=[data_file, background_file]
integration_times = files.map{|x| x.gets.chomp.split[1].to_f}
files.each{|x| 11.times{x.gets}}
sum = 0
asum = 0
bsum = 0
scale = integration_times[1]/integration_times[0]
s2 = 86400/integration_times[1]
while s=data_file.gets
  t=background_file.gets
  a=s.chomp.split
  a2=a[2].to_f
  b2=t.chomp.split[2].to_f
  diff = scale*a2 - b2
  e = a[1].to_f
  sigma=0
  if e > emin
    sum += diff 
    asum += scale*scale*a2
    bsum += b2
  end
  exit if e > emax
  sigma = asum + bsum
  print( a[1], " ", a2*scale*s2, " ", b2*s2,  " ", diff*s2, " ", sum*s2," ",
     sqrt(sigma)*s2, " ", asum/scale*s2, " ", bsum*s2, "\n")
  convert_to_bq_per_kg(sum*s2, sqrt(sigma)*s2)  if (fabs(a[1].to_f-0.8) < 0.001) 
    

end