#include <cmath>        // std::abs

float               mean(float data[], int n);
float standard_deviation(float data[], int n);

void my_pi()
{
    int seed=34545, in, neve=0, sampling=100, jjmax=100;
    float x[2], xx[jjmax],yy[jjmax],zz[jjmax], mypi[sampling];
    const float pi = 3.1415;
    
    gRandom->SetSeed(seed);  // root command
    
    
    for (int jj=0; jj<jjmax; jj++) {
        

        neve = neve + 10;
        
        
 // =================== now do this calculation 100 times, and quote mean and sigma :
        
        for (int kk=0; kk<sampling; kk++) {
                    in = 0;

// ==
    for (int ii=0; ii<neve; ii++){
        
        x[0] = gRandom->Uniform();  // root command
        x[1] = gRandom->Uniform();
        
        if( sqrt(x[0]*x[0]+x[1]*x[1]) <= 1 ) in = in + 1 ;  // root command
    }
// ==
            mypi[kk] = 4.0 * (float(in)/float(neve));
            
        }
        
// ===================
        float my_pi = mean(mypi, sampling);
        float my_pi_sig = standard_deviation(mypi, sampling);
        
        xx[jj] = neve;
        yy[jj] = my_pi;  //  abs(pi - my_pi);
        zz[jj] = my_pi_sig;
        
        cout<<jj<<" trial = "<< neve <<" my mean pi = "<<my_pi<<" Standard Deviation = " << my_pi_sig <<endl;
        
                                 }
    // create graph
//    TGraph *gr = new TGraph(jjmax,xx,yy);
      TGraph *gr = new TGraphErrors(jjmax,xx,yy,0,zz);
    TCanvas *c1 = new TCanvas("c1","Graph Draw Options",200,10,600,400);
    // draw the graph with axis, contineous line, and put a * at each point
//    gr->Draw("AC*");
        gr->Draw("ALP");
    gr->SetLineWidth(1);
    gr->SetMarkerColor(4);
    gr->SetMarkerStyle(29);
    gr->SetTitle("MC calculation of Pi with error bar");
    gr->GetXaxis()->SetTitle("MC trial events");
    gr->GetYaxis()->SetTitle("estimated Pi value");
    
}





float standard_deviation(float data[], int n)
{
    float mean=0.0, sum_deviation=0.0;
    int i;
    for(i=0; i<n;++i)
    {
        mean+=data[i];
    }
    mean=mean/n;
    for(i=0; i<n;++i)
        sum_deviation+=(data[i]-mean)*(data[i]-mean);
    return sqrt(sum_deviation/n);
}

float mean(float data[], int n)
{
    float mean=0.0, sum_deviation=0.0;
    int i;
    for(i=0; i<n;++i)
    {
        mean+=data[i];
    }
    mean=mean/n;
    for(i=0; i<n;++i)
        sum_deviation+=(data[i]-mean)*(data[i]-mean);
    return mean;
}