// My examples with Box-Muller Transformation
// Read from WEB:  http://mathworld.wolfram.com/Box-MullerTransformation.html
// Generating Gaussian Random Numbers from uniform pseudo-random numbers


void mybm()
{
    int seed=12345; 
    const int neve=10000;
    float r1[neve], a=2.0, x[3], y[3] ;
    const float pi = 3.1415;
    
    gRandom->SetSeed(seed);  // root command

    TCanvas *c1 = new TCanvas("c1","multipads",900,700); // root command for plotting many figures
    gStyle->SetOptStat(0);
    c1->Divide(2,2);
    
  TH2F * huni = new TH2F("uni", "uni for Box-Muller Transformation", 100, 5., 5., 100, -5., 5.);  // root
  TH2F * hbm  = new TH2F("Box-Muller", "Box-Muller Transformation", 100, -5., 5., 100, -5., 5.);  // root command
  TH2F * hrot = new TH2F("Box-Muller", "Transfomed Variables",100, 2., 14.0, 100, -8., 5.);  // root command

    
    for(int ii=0; ii<neve; ii++){
        
        x[0] = gRandom->Uniform();  // root command
        y[0] = gRandom->Uniform();  // root command
        
        x[1] = pow(-2*log(x[0]), 0.5) * cos(2*pi*y[0]);
        y[1] = pow(-2*log(x[0]), 0.5) * sin(2*pi*y[0]);

        float z1 = 5.0 + 0.3*x[1];
        float z2 = 7.0 + 2.0*y[1];

        float w1 = (z1+z2)/sqrt(2.0);
        float w2 = (z1-z2)/sqrt(2.0);
        
        huni->Fill(x[0],y[0]);  // root command
        hbm->Fill(x[1],y[1]);  // root command
        hrot->Fill(w1,w2);  // root command
    }
    
    c1->cd(1);
    huni->Draw();   // root command
    
    c1->cd(2);
    hbm->Draw();   // root command
    
    c1->cd(3);
    hrot->Draw();   // root command

}