/* In this ROOT function we generate a distribution according to sin(x) between 0 and pi

To run do:
root 
.L generate.C+ 
boxgenerate(10000)
anagenerate(10000)
*/

// include C++ STL headers 
#include <iostream>

using namespace std;

// ROOT library obejcts
#include <TF1.h> // 1d function class
#include <TH1.h> // 1d histogram classes
#include <TH2.h> // 2d histogram classes
#include <TRandom.h> // random generators
#include <TStyle.h>  // style object
#include <TMath.h>   // math functions
#include <TCanvas.h> // canvas object

void boxgenerate(Int_t nEvents); // box method
void anagenerate(Int_t nEvents); // analytical method
void rootfuncgenerate(Int_t nEvents); // ROOT method (a bit dangerous since we don't know what happens!)

//________________________________________________________________________
void boxgenerate(Int_t nEvents) 
{
  cout << "Generating " << nEvents << " events" << endl << endl;

  // create histogram objects
  TH2F* hThrown = new TH2F("hThrown", "Input random variables", 
			   100, 0, TMath::Pi(), 100, 0, 1);
  TH2F* hAccepted = new TH2F("hAccepted", "Accepted random variables", 
			     100, 0, TMath::Pi(), 100, 0, 1);
  TH1F* hSin = new TH1F("hSin", "Box generated sin(x) distribution", 
			100, 0, TMath::Pi());
  
  // create random generator object
  TRandom* randGen = new TRandom();

  // make a loop for the number of events
  for(Int_t n = 0; n < nEvents; n++) {
    
    if((n+1)%1000==0)
      cout << "event " << n+1 << endl;
    
    Float_t x = randGen->Rndm()*TMath::Pi(); // flat random between 0 and pi
    Float_t y = randGen->Rndm(); // flat random between 0 and 1

    hThrown->Fill(x, y); // fill 2d histogram with input values

    if(y <= TMath::Sin(x)) { // condition
      
      hAccepted->Fill(x, y); // fill 2d histogram with input values
      hSin->Fill(x); // fill our sin dist histogram
    }
  }
  
  gStyle->SetOptStat(1111);
  gStyle->SetOptFit(1111);

  // create canvas for 2d histograms (help histograms)
  TCanvas* c1 = new TCanvas("c1", "help canvas", 600, 400);
  c1->cd();
  hThrown->SetMarkerStyle(20);
  hThrown->SetMarkerColor(kBlue);
  hThrown->SetMarkerSize(1.0);
  hThrown->Draw("P");

  hAccepted->SetMarkerStyle(29);
  hAccepted->SetMarkerColor(kRed);
  hAccepted->SetMarkerSize(1.5);
  hAccepted->Draw("P SAME");  

  // create canvas for hSin
  TCanvas* c2 = new TCanvas("c2", "sin canvas", 900, 600);
  hSin->SetMinimum(0);
  hSin->Draw();

  // create 1d function
  TF1* func = new TF1("func", "[0]*sin(x)", 0, TMath::Pi());
  func->SetParameter(0, 10);
  func->SetLineColor(kRed);
  hSin->Fit(func);

  c2->SaveAs("sinx_box.gif");
}

// Questions:
// What is the efficiency of the generation?
// does this agree with theory?

//________________________________________________________________________
void anagenerate(Int_t nEvents) 
{
  cout << "Generating " << nEvents << " events" << endl << endl;
  
  // create histogram objects
  TH1F* hSin = new TH1F("hSin", "Box generated sin(x) distribution", 
			100, 0, TMath::Pi());
  
  // create random generator object
  TRandom* randGen = new TRandom();
  
  // make a loop for the number of events
  for(Int_t n = 0; n < nEvents; n++) {
    
    if((n+1)%1000==0)
      cout << "event " << n+1 << endl;
    
    hSin->Fill(TMath::ACos(1-2*randGen->Rndm())); // fill our sin dist histogram
  }
  
  gStyle->SetOptStat(1111);
  gStyle->SetOptFit(1111);
  
  // create canvas for hSin
  TCanvas* c1 = new TCanvas("c2", "sin canvas", 900, 600);
  hSin->SetMinimum(0);
  hSin->Draw();
  
  // create 1d function
  TF1* func = new TF1("func", "[0]*sin(x)", 0, TMath::Pi());
  func->SetParameter(0, 10);
  func->SetLineColor(kRed);
  hSin->Fit(func);

  c1->SaveAs("sinx_analytical.gif");
}

//________________________________________________________________________
void rootfuncgenerate(Int_t nEvents) 
{
  cout << "Generating " << nEvents << " events" << endl << endl;

  // create histogram objects
  TH1F* hSin = new TH1F("hSin", "ROOT func generated sin(x) distribution", 
			100, 0, TMath::Pi());
  
  // make a loop for the number of events
  TF1* sinFunc = new TF1("sinFunc", "sin(x)", 0, TMath::Pi());

  for(Int_t n = 0; n < nEvents; n++) {
    
    if((n+1)%1000==0)
      cout << "event " << n+1 << endl;
    
    hSin->Fill(sinFunc->GetRandom()); // fill our sin dist histogram
  }
  
  gStyle->SetOptStat(1111);
  gStyle->SetOptFit(1111);

  // create canvas for hSin
  TCanvas* c1 = new TCanvas("c1", "sin canvas", 900, 600);
  hSin->SetMinimum(0);
  hSin->Draw();

  // create 1d function
  TF1* func = new TF1("func", "[0]*sin(x)", 0, TMath::Pi());
  func->SetParameter(0, 10);
  func->SetLineColor(kRed);
  hSin->Fit(func);

  c1->SaveAs("sinx_rootfunc.gif");
}