/***************************************************************************
 * --------------------------------------------------------------------------
 * Solution to exercise PDFs and MC-Methods
 * -> Probability density functions - Monte Carlo Methods 
 * --------------------------------------------------------------------------
 * 13.06.2005, Klaus Rabbertz
 *
 * modification log:
 * kr, 06.07.2006: Removed explicit ex11 naming
 * as, 07.06.2008: Various changes
 * gq, 04.01.2011: some adaptions 
 ***************************************************************************/

#include "TFile.h" 
#include "TH1.h"
#include "TRandom.h"
#include <stdio.h>
#include <iostream>
#include <math.h>
#include "TMath.h"

// Function to genericly implement the acceptance-rejection method. 
void AcceptReject (TH1F *histo, Float_t (*f)(Float_t x, Float_t p), Float_t p,
		   Float_t xmin, Float_t xmax, int NumRand);

// Two example functions:
Float_t f1(Float_t x, Float_t p);        // p*(1*x^2) exercise 2
Float_t SawTooth (Float_t x, Float_t p); // exercise 3

int main()
{
  Int_t NumRand=100000;  // Number of random number to be filled into
                         // a histogram.

  // Create a new ROOT file.
  TFile hfile("pdfs.root","RECREATE","pdfs");
  
  TH1F *h100 = new TH1F("h100","Sawtooth distribution (transformation)",
			100,0.0,1.0);
  TH1F *h110 = new TH1F("h110","Sawtooth distribution (accept-reject)",
			100,0.0,1.0);
  TH1F *h200 = new TH1F("h200","3/8(1+x^2)",           100,-1.0,1.0);
  TH1F *h400 = new TH1F("h300","Cauchy distribution",
			100,-10.0,10.0);
  TH1F *h500 = new TH1F("h400","Sample mean of Cauchy distribution",
			100,-10.0,10.0);
  //------------------------------------------------------------------------

  Float_t dummy=0.0;

  // Exercise 1
  // (a) Fill histogram h100 with sawtooth distribution using transformation 
  //     method.
  Float_t xmax=1.0;                
  for (Int_t i=0;i<NumRand; i++) {
    dummy = gRandom->Rndm(1); 
    h100->Fill(xmax*sqrt(dummy));
  }
  // (b) Fill histogram 210 with sawtooth distribution using the 
  //     acceptance/rejection method.
  AcceptReject(110, SawTooth, 1.0, 0.0, 1.0, NumRand);
  //------------------------------------------------------------------------

  // Exercise 2 
  // Fill histogram h200 with 3/8(1+x^2) using acceptance/rejection method.
  // Put the code for this exercise in function "AcceptReject".
  AcceptReject(h200, f1, 3.0/8.0, -1.0, 1.0, NumRand);
  //------------------------------------------------------------------------

  // Exercise 3
  // Fill histogram 300 with Cauchy-distribution using the transformation
  // method.
  for (Int_t i=0; i<NumRand; i++) {
    dummy = gRandom->Rndm(1); 
    h300->Fill(tan(TMath::Pi()*(dummy-0.5)));
  }

  // Exercise 4  Fill histo 400 with sample mean of n Cauchy distributions
  Int_t   n=10;
  Float_t mean=0;
  for  (Int_t i=0; i<NumRand; i++){
    mean = 0.0;
    for (Int_t j=0; j<n; j++) {
      dummy = gRandom->Rndm(1); 
      mean += tan(TMath::Pi()*(dummy-0.5));
    } // for j
    h400->Fill(mean/(Float_t)n);  //fill normalised sample mean into histo
  } //for i
  //------------------------------------------------------------------------

  // Save all histograms in the ROOT file and close it.
  hfile.Write();
  hfile.Close();

  return 0;
}                 

//-------------------------------------------------------------------------------
void AcceptReject (TH1F *histo, Float_t (*f)(Float_t x,Float_t p), Float_t p,
		   Float_t xmin, Float_t xmax, int NumRand){
  /* transform uniform random numbers to numbers following a distribution
   * described by f
   * Input variables: 
   * histo  : ID of histogram to fill
   * f      : function the random numbers should follow
   * xmin   : left border of interval
   * xmax   : right border of interval
   * NumRand: number of random numbers to generate 
   */


  // local variables
  Float_t fmax=0.0;  // maximum value of function f in (xmin, xmax)
  Float_t u,v;       // generated random numbers in (0,1)
  Float_t x,y;       // random numbers in (xmin, xmax,) (0, fmax)
  Int_t NumBins;     // Number of bins of histogram histo
  Float_t step, dummy;
  int flag=0;

  // Find the maximum of f in (xmin,xmax) by stepping through the
  // interval.
  NumBins = histo->GetNbinsX();  // Get the number of bins of histogram.
  step = (xmax-xmin)/(Float_t)NumBins;  //step width
  for (Float_t x=xmin; x<=xmax; x+=step){
    dummy = (*f)(x,p);
    if (dummy>fmax) {
      fmax = dummy;  // new maximum found
    }
  }

  std::cout << "Maximum of function is: "<< fmax << std::endl;

  for (int i=0; i<NumRand; i++) {
    while (flag == 0) {
      u  = gRandom->Rndm(1); //x - coordinate
      x  = xmin + u*(xmax-xmin);
      v  = gRandom->Rndm(1);
      y  = v*fmax;             //y - coordinate
      
      if (y <= f(x,p)) {
	flag =  1;
	histo->Fill(x);
	//cout << "accepted x= " << x << " f(x)= " << f(x) << " y= " << y << endl;
      } //if
    } // while 

    flag = 0;  //reset flag
  } // for

}

//-------------------------------------------------------------------------------
Float_t f1(Float_t  x, Float_t p) {
  return p*(1+x*x);
}

//-------------------------------------------------------------------------------
Float_t SawTooth(Float_t x, Float_t p) {
  // p is a parameter (here: xmax)
  Float_t rtVal=0.0; // return value 

  if (x<p) {
    rtVal = (Float_t)2.0*x/(p*p);
  }
    
  return rtVal;
}