//
// Program: buffon.C
//
// Description: Simulate the Buffon's needle problem with a Monte Carlo
//
// Author: A. Garfagnini
//
// Date: 30 November 2013
//
#include "TApplication.h"     // SA version
#include "TROOT.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TH2F.h"
#include "TGraph.h"
#include "TLine.h"
#include "TLegend.h"

#include "boost/random.hpp"


#include <iostream>
#include <cmath>
#include <fstream>
#include <unistd.h>
#include <libgen.h>
#include <stdlib.h>

using namespace std;

void usage(const char * pname)
{
    cerr << "Usage: " << pname << " [-o file_name]" << endl;
}

void buffon()
{
  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(0);
  gStyle->SetOptTitle(kFALSE);
  gStyle->SetLabelOffset(0.01,"x");
  gStyle->SetLabelOffset(0.005,"y");
  gStyle->SetTitleOffset(1.2,"y");
  gStyle->SetTitleOffset(1.2,"x");
}


int main(int argc, char * argv[])
{
   int opt;

    // Handle command line options
    while((opt = getopt(argc, argv, "h?")) != -1) {

        switch (opt) {

        case 'h':
        case '?':
        default:
            usage( basename( argv[0] ) );
            return 1;
        }
    }

    // Handle user data input
    cout << "Random generator seed: ";
    int seed;
    cin >> seed;

    boost::uniform_real<> uni_u(0.0, 1.0);
    boost::uniform_real<> uni_v(0.0, M_PI);
    boost::mt19937 rng;
    rng.seed(seed);
    boost::variate_generator < boost::mt19937 &,
                               boost::uniform_real<> > uni_rng_u(rng, uni_u);
    boost::variate_generator < boost::mt19937 &,
                               boost::uniform_real<> > uni_rng_v(rng, uni_v);


    // Startup ROOT
    TApplication myapp("myapp", &argc, argv);

    gROOT->Reset();
    buffon();

    TCanvas * c3 = new TCanvas("c3","test");
    gPad->SetTopMargin(0.01);
    gPad->SetBottomMargin(0.12);
    gPad->SetRightMargin(0.01);
    gPad->SetLeftMargin(0.12);

    TH2F * hs = new TH2F("hs", "", 50, 0., 49999., 50, 0.52, 0.79);
    hs->Draw();
    hs->GetXaxis()->SetTitle("ntry");
    hs->GetXaxis()->SetLabelSize(0.04);
    hs->GetYaxis()->SetTitle("result");

    double result[100] = {0.0};
    double ntry[100] = {0.0};
    int n_tot = 50000;
    int n_hit = 0;
    int i = 0;
    while (i<n_tot) {

        i++;

        double u, v;

        u = uni_rng_u();
        v = uni_rng_v();

        // cout << " u = " << u << " v = " << v << endl;

        if ( u < sin(v) ) {
            n_hit++;
        }

        if (i%500 == 0) {
            int j = i/500 - 1;
            ntry[j] = i;
            result[j] = double(n_hit)/ntry[j];

            //cout << j << " : " << ntry[j]
            //     << " : " << result[j] << endl;
        }

    }

    TGraph * g = new TGraph(100, ntry, result);
    g->SetMarkerColor(1);
    g->SetLineColor(1);
    g->SetMarkerStyle(25);
    g->SetMarkerSize(0.9);

    g->Draw("plsame");

    TLine * m_2_pi = new TLine(0., M_2_PI, 50100., M_2_PI);
    m_2_pi->SetLineColor(2);
    m_2_pi->SetLineStyle(1);
    m_2_pi->Draw();


    c3->Update();

    myapp.Run();

    return 0;
}