//
// somma di variabili uniformi
//



void MCDid1() {

  TCanvas *c1=new  TCanvas("c1", "distribution", 1);
     
     // create fill and draw
  gStyle -> SetOptStat(kFALSE) ;
  TH1F *h1 =  new TH1F("h1","gaussians", 100, -4., 4);
  Int_t i,k, nuni;
  Float_t unic=0.;
  nuni=100;
  for( i =1; i <=10000;i++) {
    unic=0.;
    for(k=1; k<=nuni; k++) {
      unic += gRandom->Uniform();
    }
    unic = (unic/nuni -0.5)*sqrt(12*nuni);
    h1->Fill(unic);
  }
  h1->Draw();
  c1->Update();

  // cumulative function
  TH1F *hint1 = new TH1F("hint1", "integral", 100, -4.,4.);
  Float_t sum = 0.;
  for( i =1; i <=100;i++) {
       sum += h1->GetBinContent(i);
       hint1->SetBinContent(i,sum);
  }

     // scale coordinates
  Float_t rightmax = 1.1*hint1->GetMaximum();
  Float_t scale = gPad->GetUymax()/rightmax;
  hint1->SetLineColor(kRed);
  hint1-> Scale(scale);
  hint1-> Draw("same");

  // right hand axis
  TGaxis *axis = new TGaxis(gPad->GetUxmax(), gPad->GetUymin(),
                            gPad->GetUxmax(),
                            gPad->GetUymax(),0,rightmax,510,"+L",0);
  axis->SetLineColor(kRed);
  axis->SetTextColor(kRed);
  axis->Draw();

}








void MCDid2() {

  //
  // distribuzione binomiale b(x; 10.,0.2)
  //

  TCanvas *c1=new  TCanvas("c1", "distribution", 1);
   

  // binomial parameters: p = probability, nuni = trials;
  //
  Float_t prob=0.50;
  Int_t nuni=30;
  
     // create fill and draw
  gStyle -> SetOptStat(kFALSE) ;
  TH1F *h1 =  new TH1F("h1","binomial", nuni+1, 0., (Double_t) nuni);
  Int_t i,k, nuni;
  Float_t unic=0.;

  for( i =1; i <=10000;i++) {
    unic=0.;
    for(k=1; k<=nuni; k++) {
      if(gRandom->Uniform() < prob) unic+=1.;
    }
    h1->Fill(unic);
  }
  h1->Draw();
  c1->Update();

  // cumulative function
  TH1F *hint1 = new TH1F("hint1", "integral", nuni+1, 0.,(Double_t) nuni);
  Float_t sum = 0.;
  for( i =1; i <=nuni+1;i++) {
       sum += h1->GetBinContent(i);
       hint1->SetBinContent(i,sum);
  }

     // scale coordinates
  Float_t rightmax = 1.1*hint1->GetMaximum();
  Float_t scale = gPad->GetUymax()/rightmax;
  hint1->SetLineColor(kRed);
  hint1-> Scale(scale);
  hint1-> Draw("same");

  // right hand axis
  TGaxis *axis = new TGaxis(gPad->GetUxmax(), gPad->GetUymin(),
                            gPad->GetUxmax(),
                            gPad->GetUymax(),0,rightmax,510,"+L",0);
  axis->SetLineColor(kRed);
  axis->SetTextColor(kRed);
  axis->Draw();

}




void MCDid3() {

  // 
  // chi-square distribution
  //

  TCanvas *c1=new  TCanvas("c1", "distribution", 1);
  

  //degrees of freedom
  nuni=5;
   
     // create fill and draw
  gStyle -> SetOptStat(kFALSE) ;
  TH1F *h1 =  new TH1F("h1","chi square", 100, 0, 4.*nuni);
  Int_t i,k, nuni;
  Float_t unic=0.;
  

  for( i =1; i <=10000;i++) {
    unic=0.;
    for(k=1; k<=nuni; k++) {
      unic += pow(gRandom->Gaus(0,1),2.);
    }
    h1->Fill(unic);
  }
  h1->Draw();
  c1->Update();

  // cumulative function
  TH1F *hint1 = new TH1F("hint1", "integral", 100,0, 4.*nuni);
  Float_t sum = 0.;
  for( i =1; i <=100;i++) {
       sum += h1->GetBinContent(i);
       hint1->SetBinContent(i,sum);
  }

     // scale coordinates
  Float_t rightmax = 1.1*hint1->GetMaximum();
  Float_t scale = gPad->GetUymax()/rightmax;
  hint1->SetLineColor(kRed);
  hint1-> Scale(scale);
  hint1-> Draw("same");

  // right hand axis
  TGaxis *axis = new TGaxis(gPad->GetUxmax(), gPad->GetUymin(),
                            gPad->GetUxmax(),
                            gPad->GetUymax(),0,rightmax,510,"+L",0);
  axis->SetLineColor(kRed);
  axis->SetTextColor(kRed);
  axis->Draw();

}




 // 
 //        Maxwell/Boltzmann distributions
 //


Double_t maxw(Double_t *x, Double_t *par) {
  Double_t xi = x[0];
  return par[0]*sqrt(2./3.14159)*xi*xi * exp(-xi*xi/2.);
}

Double_t boltz(Double_t *x, Double_t *par) {
  Double_t xi = x[0];
  return par[0]* 2.*sqrt(1./3.14159)* sqrt(xi) * exp(-xi);
}

void MCDid4() {

  gStyle->SetOptStat(1000001);
  TCanvas *c1=new  TCanvas("c1", "Maxwell distribution", 1);
  TCanvas *c2=new  TCanvas("c2", "Boltzmann distribution", 1);
  
   
     // create fill 
  gStyle -> SetOptStat(kTRUE) ;
  TH1F *hm =  new TH1F("hm","maxwell", 100, 0, 6.);
  TH1F *hb =  new TH1F("hb","boltzmann", 100, 0, 12.);

  Int_t i,k, nuni=3;
  Float_t unic=0.;
  

  for( i =1; i <=100000;i++) {
    unic=0.;
    for(k=1; k<=nuni; k++) {
      unic += pow(gRandom->Gaus(0,1),2.);
    }
    hb->Fill(unic/2.);
    hm->Fill(sqrt(unic));
  }
  
  TF1 *fit = new TF1("maxw", maxw, 0,6,1);
  fit-> SetParameter(0,5000.);  
  fit-> SetLineColor(4);

  TF1 *fitb = new TF1("boltz", boltz, 0,12,1);
  fitb-> SetParameter(0,3000.);  
  fitb-> SetLineColor(4);
 
  c1->cd();
  //hm->Draw("error");
  hm-> Fit("maxw","V","error",0.,6.); 
  c1->Update();
  c2->cd();
  //hb->Draw();
  hb-> Fit("boltz","V","error",0.,12.); 
  c2->Update();

}



/*  ======================================================================
 
  TMinuit m(1);
 
  m.DefineParameter(0,"sigma",  1., 0.01, 0.,  1000.); 

  void fcn(Int_t &, Double_t *, Double_t &, Double_t *, Int_t);
  m.SetFCN(fcn);   
  m.Migrad();

  Double_t x0,x1,x2,x3,x4,x5,x6,x7,x8,x9;
  Double_t e0,e1,e2,e3,e4,e5,e6,e7,e8,e9;
  Int_t err0=m.GetParameter(0,x0,e0);
 
  Int_t nx = ghMC->GetNbinsX();
  Double_t Dx = 40./(Double_t) nx;
  Stat_t iMC=ghMC->Integral(); 
  
  Double_t xc,xa,yc,ff1,ff2,fmax;

  fmax = x0;
  for (Int_t i=0;i<1000000;i++)  {
    xc=-16 + 32.*gRandom->Rndm(); 
    xa= TMath::Abs(xc); 
    ff1 = pow(xa/x1,x2);
    ff2 = pow(xa/x3,x2);
    yc = x0*TMath::Exp(-ff1 );
    if(gRandom->Rndm() < yc/fmax) rej.Fill(xc);
  }

 Stat_t totrej=rej->Integral();

 rej.Scale(iMC/totrej);
 rej.Draw("same");


void fcn(Int_t& npar, Double_t *g, Double_t& fval, Double_t *x, Int_t iflag) {
  Double_t a=x[0];

  Double_t d1=x[1];
  Double_t e1=x[2];

  Double_t d2=x[3];
  Double_t e2=x[4];

  Int_t nx = ghMC->GetNbinsX();
 
  Double_t Dx = 40./(Double_t) nx;
  fval=0;
  Stat_t iMC=ghMC->Integral();       
 
  for (Int_t i=0;i<nx;i++)  {
    Double_t xi= -20 + i*Dx + 0.5*Dx;
    xi = TMath::Abs(xi);
    Double_t wMC=ghMC->GetBinContent(i);
    if(wMC>10.) {
      Double_t xx1 = pow(xi/d1,e1);
      Double_t xx2 = pow(xi/d2,e2);
      Double_t num = wMC - a*Dx*iMC*TMath::Exp(-xx1);
      Double_t den = wMC;
      if (den) fval += num*num/den;
    }
  }
}

  ==================================================================  */




 // 
 //        Monte Carlo Integration  
 //        of the Error function from xmin=-3.5 and xmax
 //


Double_t func(Double_t xi) {
  return  sqrt(1./6.2832) * exp(- 0.5 *xi*xi);
}

Double_t funch(Double_t *x, Double_t *par) {
  Double_t xi = x[0];
  return  par[0]*sqrt(1./6.2832) * exp(- 0.5 *xi*xi);
}

void MCDid5(Double_t xmin, Double_t xmax) {

  // set parameters
  //
  Double_t fmax,area; 
  Int_t n=50;
  fmax = func(0.);
  
  area= fmax*(xmax-xmin);

  gStyle->SetOptStat(1000001);
  TCanvas *c1=new  TCanvas("c1", "Standard Gaussian", 1);
  TCanvas *c2=new  TCanvas("c2", "Sampled Gaussian", 1);

  Int_t succ=0, Ntot=500000;
  Double_t sum=0., sum2=0.;

  Double_t x[50], y[50];
  for(Int_t i=0; i<n;i++){
    x[i] = xmin + (xmax-xmin)*i/(Double_t) n;
    y[i] = func(x[i]);
  }
  TGraph  *g1=new TGraph(n,x,y);
  c1->cd();
  g1->SetLineColor(2);
  g1->Draw("Al");

     // create histogram 
  gStyle -> SetOptStat(kTRUE) ;
  TH1F *hs =  new TH1F("hs","sampled gaussian", 100, xmin, xmax);
 
  Int_t i,k;
  
  // rejection method
  //
 
  for( i =1; i <=Ntot;i++) {
    Double_t xc = xmin + (xmax-xmin)*gRandom->Rndm();
    if(fmax*gRandom->Rndm() < func(xc)) {
      hs->Fill(xc);
      succ++;
    }
    sum += func(xc);
    sum2 +=  func(xc)* func(xc);
  }
  Double_t integ,integ1, prob, error, error1,fNtot;
  fNtot = (Double_t) Ntot;
  prob=  (Double_t) succ/ fNtot;
  integ=  area*prob;
  error = sqrt( integ*(area-integ)/ fNtot ); 
 
  integ1 = (xmax-xmin)*sum/ fNtot ; 
  error1 = (pow((xmax - xmin),2)/(fNtot*(fNtot-1)))*(sum2 - sum*sum/fNtot) ;
  error1=sqrt(error1);
 

  TF1 * g2 = new TF1("funch",funch, xmin,xmax,1);
  g2-> SetLineColor(2); 

  // beware the integ value at the denominator
  g2-> SetParameter(0,(Double_t) succ * (xmax-xmin)/ (integ*100.) );  
 
  c2->cd();
  //hs->Draw();
  hs-> Fit("funch","V","error",xmin,xmax);
 
  c2->Update();

  cout << " ---------------------------------------------------------"<<endl;
  cout << " I  value of the integral between" <<xmin<<" and "
       <<xmax <<" = " << integ << " +- "<< error <<endl;
  cout << " II value of the integral between" <<xmin<<" and "
       <<xmax <<" = " << integ1 << " +- "<< error1 <<endl;
  cout << " ---------------------------------------------------------"<<endl;

}