import ROOT
import math


class Generator():
	def __init__(self):
		self.rand = ROOT.TRandom()

	def get_rv(self):
		x = (x_up - x_down) * self.rand.Uniform() + x_down
		y = (y_up - y_down) * self.rand.Uniform() + y_down
		return x, y


class Integrator():
	def __init__(self, func, rv_gen):
		self.func = func
		self.rv_gen = rv_gen
		self.series_count = 0
		self.n_Eval = 0.0
		self.n_Eval2 = 0.0

	def iterate(self, num_iterations):
		self.series_count += 1

		n_Ef = 0.0
		n_Ef2 = 0.0
		for i in range(num_iterations):
			x, y = self.rv_gen.get_rv()
			#print x, y
			func_val = self.func.Eval(x, y)
			n_Ef += func_val
			n_Ef2 += func_val**2

		Ef = n_Ef / num_iterations
		self.n_Eval += Ef
		self.n_Eval2 += Ef**2
		Ef2 = n_Ef2 / num_iterations
		sd = ( Ef2 - Ef**2 )**0.5
		return Ef, sd

	def get_ev_sd(self):
		Eval = self.n_Eval / self.series_count
		Eval2 = self.n_Eval2 / self.series_count
		sd = ( Eval2 - Eval**2 )**0.5
		return Eval, sd


x_down = 0.0
x_up = 0.5
y_down = 1.0
y_up = 3.0
num_series = 10
num_iterations_per_series = 1000
func = ROOT.TF2("func", "exp( x^2 + y^2 ) / ( 1 + x * y**0.5 / 2 )", x_down, x_up, y_down, y_up)


rvgen = Generator()
integrator = Integrator(func, rvgen)
for series in range(num_series):
	val, sd = integrator.iterate(num_iterations_per_series)
	accval, accsd = integrator.get_ev_sd()
	print "series:", series, "\tval:", val, "\tsd:", sd, "\taccval:", accval, "\taccsd:", accsd


raw_input("Press ENTER to exit.")