#!/usr/bin/env python

# ---------------------------------------------------------------------------------- #
# 
#  Python/ROOT macro for illustrating how to generate random numbers following a
#  specific PDF using uniformly distributed random numbers. Both the Accept-Reject
#  (Von Neumann) and transformation methods are used.
#
#  Author: Troels C. Petersen (NBI)
#  Email:  petersen@nbi.dk
#  Date:   11th of September 2014
# 
# ---------------------------------------------------------------------------------- #

from ROOT import *        # Import ROOT libraries (to use ROOT)
import math               # Import MATH libraries (for extended MATH functions)

gStyle.SetOptStat("emr")
gStyle.SetOptFit(1111)

SavePlots = False         # Determining if plots are saved or not

r = TRandom3()
r.SetSeed(1)

Npoints = 10000



#----------------------------------------------------------------------------------- #
# Problem 1: Produce random points following f(x) = exp(-x/3), x in [0,inf]
#----------------------------------------------------------------------------------- #

Hist_exp = TH1F("Hist_exp", ";x (exponentially distributed);Frequency", 100, 0.0, 20.0)

for i in range (Npoints) :
    x = -3.0 * log(r.Uniform())
    Hist_exp.Fill(x)

canvas1 = TCanvas("canvas1","", 50, 50, 600, 400)
Hist_exp.SetLineColor(kBlue)
Hist_exp.SetLineWidth(2)
Hist_exp.Draw()

# Normalization is Npoints * binwidth = 2000:
f_exp = TF1("f_exp", "2000 * 1.0/3.0 * exp(-x/3.0)", 0, 20.0)
f_exp.Draw("same")

if (SavePlots) : canvas1.SaveAs("Hist_exp.pdf")




#----------------------------------------------------------------------------------- #
# Problem 2: Produce random points following f(x) = x*cos(x), x in [0,pi/2]
#----------------------------------------------------------------------------------- #

pi = TMath.Pi()

Hist_xcos = TH1F("Hist_xcos", ";x (exponentially distributed);Frequency", 100, 0.0, pi/2.0)

Nhit = 0
for i in range (Npoints) :
    x = pi/2.0 * r.Uniform()     # Range that f(x) is defined/wanted in
    y = 1.0 * r.Uniform()        # Upper bound for function values (a better bound exists!)
    if (y < x*cos(x)) :          # Test if x-value is accepted or not
        Nhit += 1
        Hist_xcos.Fill(x)

integral = float(Nhit) / float(Npoints) * pi/2.0 * 1.0
eintegral = -1.0                 # Calculate this yourself!
print "  Integral of f(x) = x*cos(x), x in [0,pi/2] is: %7.4f +- %6.4f"%(integral, eintegral)

canvas2 = TCanvas("canvas2","", 80, 80, 600, 400)
Hist_xcos.SetLineColor(kBlue)
Hist_xcos.SetLineWidth(2)
Hist_xcos.Draw()

# Normalization here is Npoints * binwidth = pi*50:
f_xcos = TF1("f_xcos", "3.1416 * 50 * 0.570796*x*cos(x)", 0, pi/2.0)
f_xcos.Draw("same")

if (SavePlots) : canvas2.SaveAs("Hist_xcos.pdf")



#----------------------------------------------------------------------------------- #
# Problem for you: Try to combine the two yourself:
#----------------------------------------------------------------------------------- #
#
# Estimate the integral of exp(-x/3.0)*cos(x)^2 in the interval [0,inf] using a
# combination of Transformation and Hit-and-Miss method.

Nhit = 0
for i in range ( Npoints ) :
      
    # Well, the inside of the loop is up to you!
    Nhit += 1


integral  = -1.0     # Well, put your value here
eintegral = -1.0     # ...
print "  Integral of f(x) = exp(-x/3)*cos(x)^2, x in [0,inf] is: %7.4f +- %6.4f"%(integral,eintegral)





raw_input( ' ... Press enter to exit ... ' )

# ----------------------------------------------------------------------------------- #
# 
# The purpose of this exercise is to combine Transformation and Hit-and-Miss methods
# to evaluate an integral (here only in 1D, to make it easy/illustrative).
# 
# Questions:
# ----------
#  1) Find the integral of exp(-x/3.0)*cos(x)^2 in the interval [0,inf] using a
#     combination of the Transformation and Hit-and-Miss methods.
# 
# ----------------------------------------------------------------------------------- #