/////////////////////////
// accept_reject
/////////////////////////

#include "TStyle.h"
#include "TTimer.h"
#include "Riostream.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TROOT.h"
#include "TGraphErrors.h"
#include "TFrame.h"
#include "TAxis.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TRandom.h"
#include "TF1.h"
#include "TFormula.h"
#include "TFile.h"
#include "TTree.h"
#include "TMath.h"
#include "TMultiGraph.h"
#include "TPaveLabel.h"
#include "TLatex.h"
#include "TLegend.h"

Double_t constant (Double_t* x, Double_t* par) {

    return par[0];

}

Double_t fasym (Double_t* x, Double_t* par) {

    return par[0] * ( 1 + par[1]*TMath::Cos(x[0]) ) * 0.5/TMath::Pi();

}

int main(){

	using namespace std;

  // Parametri per la simulazione.
  const int Nrep  = 1000; // Numero di valori casuali da generare.
  double a = 0.5; // Parametro della funzione di distribuzione f(x) = (1+a*cos(x))/(2*pi).
  double max = 0.5*(1+TMath::Abs(a))/TMath::Pi(); // Massimo di f(x) nell'intervallo [0:2*pi].

  // Parametri per l'estetica dei grafici.
  int lf=133;  // fonts scales
  int tf=13;   // fonts label
  int cf=23;  // legend label
  int ls=20;   // size  scales
  int ts=25;   // size  label
  int cs=25;   // size  label

  double xOffset = 1.2;
  double yOffset = 1.0;
  double ms=1.2;
  int mStylepos=20, mStyleneg=22;

  TLatex *t2 = new TLatex();
  t2 = new TLatex();
  t2->SetNDC();
  t2->SetTextAlign(12);
  t2->SetTextFont(lf);
  t2->SetTextSize(ls);
  t2->SetTextColor(kBlack);

  // Assi del grafiuco di f(x).
  TH2F *hax = new TH2F("hax","",10,0,2*TMath::Pi(),10,0,1.2*max);
  hax->SetStats(kFALSE);
  hax->GetXaxis()->SetNdivisions(10);
  hax->GetYaxis()->SetNdivisions(10);
  hax->GetXaxis()->SetLabelFont(lf);
  hax->GetYaxis()->SetLabelFont(lf);
  hax->GetXaxis()->SetLabelSize(ls);
  hax->GetYaxis()->SetLabelSize(ls);
  hax->GetXaxis()->SetTitleFont(tf);
  hax->GetYaxis()->SetTitleFont(tf);
  hax->GetXaxis()->SetTitleSize(ts);
  hax->GetYaxis()->SetTitleSize(ts);
  hax->GetXaxis()->SetTitleOffset(1.0);
  hax->GetYaxis()->SetTitleOffset(1.0);
  hax->GetYaxis()->SetTitle( "f(x)" );
  hax->GetXaxis()->SetTitle( "x" );
  hax->SetLineColor(kBlack);
  hax->SetLineWidth(2);

  // Istogramma dei valori accettati per la variabile casuale x.
  TH1F * hist = new TH1F("hist","",16,0,2*TMath::Pi());
  hist->GetXaxis()->SetNdivisions(10);
  hist->GetYaxis()->SetNdivisions(10);
  hist->GetXaxis()->SetLabelFont(lf);
  hist->GetYaxis()->SetLabelFont(lf);
  hist->GetXaxis()->SetLabelSize(ls);
  hist->GetYaxis()->SetLabelSize(ls);
  hist->GetXaxis()->SetTitleFont(tf);
  hist->GetYaxis()->SetTitleFont(tf);
  hist->GetXaxis()->SetTitleSize(ts);
  hist->GetYaxis()->SetTitleSize(ts);
  hist->GetXaxis()->SetTitleOffset(1.0);
  hist->GetYaxis()->SetTitleOffset(1.0);
  hist->GetYaxis()->SetTitle( "counts" );
  hist->GetXaxis()->SetTitle( "x" );
  hist->SetFillColorAlpha(kGreen,0.8);
  hist->SetLineColor(kBlack);
  hist->SetLineWidth(2);


  // Funzione f(x).
  TF1* fx = new TF1("fx",fasym,0,2*TMath::Pi(),2);
  fx -> SetParameters(1, a);
  fx->SetLineColor(kBlack);
  fx->SetLineStyle(1);
  fx->SetLineWidth(2);

  // Funzione per il calcolo dei conteggi attesi.
  TF1* fexp = new TF1("fexp",fasym,0,2*TMath::Pi(),2);
  fexp -> SetParameters(1, a);
  fexp->SetLineColor(kBlack);
  fexp->SetLineStyle(1);
  fexp->SetLineWidth(2);

  // Funzione costante pari a f_max, per il disegno del rettangolo
  // che racchiude f(x).
  TF1* fc = new TF1("fc",constant,0,2*TMath::Pi(),1);
  fc->SetParameter(0,max);
  fc->SetLineColor(kBlue);
  fc->SetLineStyle(2);

  double ygen[2]={0.1,0.1};
  double xgen[2]={TMath::Pi(),1.2*TMath::Pi()};
  TGraph *gr = new TGraph(2,xgen,ygen);
  gr->SetMarkerStyle(20);
  gr->SetMarkerSize(4.0);
  gr->SetMarkerColor(kGreen);
  gr->RemovePoint(1);

  // Array con le coordinate dei punti accettati,
  // in particolare xtot_in contiene tutti i valori
  // accettati della variable casuale x.
  double xtot_in[Nrep]={0.0};
  double ytot_in[Nrep]={0.0};
  vector<double> xtot_out_vec;
  vector<double> ytot_out_vec;

  TCanvas *c1 = new TCanvas("c1","accept reject",1200,600);
  c1->cd();
  c1->Update();
  gPad->SetBottomMargin(0.2);
  gPad->SetLeftMargin(0.2);
  c1->Divide(2,0);
  c1->cd(1);
  hax->Draw();
  fx->Draw("same");
  fc->Draw("same");
  c1->cd(2);
  hist->Draw();
  int Ngen = 0;
  double rndx, rndfx;
  TRandom *rnd = new TRandom(); // Generatore di numeri casuali.
  double dx = hist->GetBinWidth(10);
  for(int irep = 0; irep < Nrep; irep++){
    while(1){
      rndx = rnd -> Uniform(2* TMath::Pi());  // Genero x con distribuzione uniforme in [0:2*pi].
      rndfx = rnd -> Uniform(max);  // Genero y con distribuzione uniforme in [0:f_max].
      Ngen += 1;  // Conto le coppie generate. 
      gSystem->ProcessEvents();
      gr->SetPoint(0,rndx,rndfx);
      if(rndfx < fx->Eval(rndx)) {  // Valuto la condizione y < f(x).
          hist->Fill(rndx); // Se la condizione e' vera, accetto x e lo inserisco nell'istogramma.
          fexp->SetParameters(hist->GetEntries()*dx,a);
          xtot_in[irep] = rndx; // Memorizzo tutti i valori di x accettati per un eventuale uso successivo.
          ytot_in[irep] = rndfx;
          c1->cd(1);
          gr->SetMarkerColor(kGreen);
          gr->Draw("p");
          c1->cd(2);
          hist->Draw("HIST,E");
          fexp->Draw("same");
          gSystem->ProcessEvents();
          c1->Update();
          break;
        } else {
          xtot_out_vec.push_back(rndx);
          ytot_out_vec.push_back(rndfx);
          c1->cd(1);
          gr->SetMarkerColor(kRed);
          gr->Draw("p");
          gSystem->ProcessEvents();
          c1->Update();
        }
    }
  }

  cout << "Efficienza del metodo = " << Nrep*1.0/Ngen << endl;

  // Le istruzioni seguenti servono per disegnari tutte le coppie di
  // punti (x,y) generate la cui x e' stata accettate, e tutte le coppie
  // rigettate.
  gr->SetPoint(0,1000,1000);
  TGraph *grIn = new TGraph(Nrep,xtot_in,ytot_in);
  grIn->SetMarkerStyle(20);
  grIn->SetMarkerSize(2.0);
  grIn->SetMarkerColor(kGreen);
  grIn->RemovePoint(1);

  TGraph *grOut = new TGraph();
  grOut->SetMarkerStyle(20);
  grOut->SetMarkerSize(2.0);
  grOut->SetMarkerColor(kRed);
  grOut->RemovePoint(1);

  for(int i = 0; i < xtot_out_vec.size(); i++){
    grOut->SetPoint(i,xtot_out_vec[i],ytot_out_vec[i]);
  }

  TLegend *leg = new TLegend(0.2,0.80,0.4,0.9);
  leg->AddEntry(grIn,"accepted","p");
  leg->SetTextFont(12);
  leg->SetTextSize(0.04);
  leg->SetLineColor(1);
  leg->SetLineStyle(1);
  leg->SetLineWidth(0);
  leg->SetFillColor(0);
  leg->SetFillStyle(0);

  TLegend *leg2 = new TLegend(0.4,0.80,0.6,0.9);
  leg2->AddEntry(grOut,"rejected","p");
  leg2->SetTextFont(12);
  leg2->SetTextSize(0.04);
  leg2->SetLineColor(1);
  leg2->SetLineStyle(1);
  leg2->SetLineWidth(0);
  leg2->SetFillColor(0);
  leg2->SetFillStyle(0);

  c1->cd(1);
  grIn->Draw("p");
  grOut->Draw("p");
  leg->Draw("same");
  leg2->Draw("same");

  return 0;
}