/***************************************************************************
 * --------------------------------------------------------------------------
 * Template 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
 ***************************************************************************/

#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 1
Float_t SawTooth (Float_t x, Float_t p); // exercise 2

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

  // Create a new ROOT file.
  TFile hfile("pdfs.root","RECREATE","pdfs");

  // Book histograms for exercises 25.1 to 25.3 and 26.
  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 *h300 = new TH1F("h300","Cauchy distribution",  100,-10.0,10.0);
  TH1F *h400 = new TH1F("h400","Sample mean of Cauchy distribution", 
			100,-10.0,10.0);
  //------------------------------------------------------------------------

  // Exercise 1
  // (a) Fill histogram h100 with sawtooth distribution using transformation 
  //     method.



  // (b) Fill histogram 110 with sawtooth distribution using the 
  //     acceptance/rejection method.
  AcceptReject(h110, 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.


  // Exercise 4 Fill histo 400 with sample mean of n Cauchy distributions


  //------------------------------------------------------------------------

  // 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
     p      : parameter value to be used in function call of f
   * xmin   : left border of interval
   * xmax   : right border of interval
   * NumRand: number of random numbers to generate 
   */

  // Find the maximum of fmax in (xmin,xmax) by stepping through the
  // interval.
  Float_t fmax=0.0;  // maximum value of function f in (xmin, xmax)
  Int_t NumBins = histo->GetNbinsX(); // Get the number of bins of histogram.
  Float_t step = (xmax-xmin)/(Float_t)NumBins;  //step width
  for (Float_t xp=xmin; xp<=xmax; xp+=step){
    
  }
  std::cout << "Maximum of function is: "<< fmax << std::endl;


  for (int i=0; i<NumRand; i++) {

  } 

}

//-------------------------------------------------------------------------------
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 


  return rtVal;
}