//
// verifica della relazione esistente tra
// una variabile aleatoria X e una variabile aleatoria y=f(x)
//
#include "TRandom.h"
const double PI =3.1415;
const double ThetaR=0.6283185;
// funzione per calcolare l'indice di rifrazione
double f(double x) {
    return sin(x)/sin(ThetaR);
}
void verifica2(){

  double ThetaR=0.6283185;
  // media della distribuzione degli angoli di incidenza
  double  mX=1.56;
// Sigma della distribuzione dei raggi
  double sX=0.08;

  TH1F* hn = new TH1F("hn","Misure n", 100,1.7-5.*0.01, 1.7+5.*0.001);
  TH1F* hrifang = new TH1F("hrif","Angolo di rifrazione", 100,mX-5.*sX,mX+5*sX);


// numero di esperimenti
    int N=100000;
    TRandom R(123456);
    double mediaX=0,mediaY=0;
    double mediaX2=0,mediaY2=0;
    double varYMC = 0.;
// loop principale per simulare gli N esperimenti
    for (int k=0;k<N;k++) {
// misura della grandezza indipendente ...
        double x=R.Gaus(mX,sX);
	hrifang->Fill(x);
// grandezza derivata 
        double y=f(x);
	hn->Fill(y);
// per il calcolo di medie e varianze....
        mediaX+=x;
        mediaY+=y;
        mediaX2+=x*x;
        mediaY2+=y*y;
	varYMC += ((y-1.69575)*(y-1.69575));
    }
// medie e varianze delle due grandezze
    mediaX=mediaX/N;
    mediaY=mediaY/N;
    mediaX2=mediaX2/N;
    mediaY2=mediaY2/N;
    varYMC = varYMC/N;
    double varX=(mediaX2-mediaX*mediaX);
    double varY=(mediaY2-mediaY*mediaY);
    double errX=sqrt(varX)/sqrt(N);  // perchè questa formula????
    double errY=sqrt(varY)/sqrt(N);  // perchè questa formula????
    double errYMC = sqrt(varYMC)/sqrt(N);
// calcolo teorico di media e errore sulla media della Y usando X
    double varYth=pow(cos(mediaX)/sin(ThetaR),2)*errX*errX;
    double Yth=f(mediaX);
    double correct=-0.5*(sin(mediaX)/sin(ThetaR))*varX;

    double vattesocorr = Yth+correct;
    TCanvas* c1 = new TCanvas("c1","c1",1000,450);
    c1->Divide(2,1);
    c1->cd(1);
    hrifang->Draw();
    c1->cd(2);
    hn->Draw();

// output
    printf("\n\n\n media di X e sigma %g %g\n",mediaX,sqrt(varX));
    printf("Y=PI*X^2\n");
    printf("varY=(2*PI*x)^2*varX\n\n");
    printf("Risultati attesi\n n=%g errore=%g \n\n\n",Yth, sqrt(varYth));
    printf("Risultati ottenuti con la simulazione\n indice di rifrazione ed errore %g %g = %g\n",mediaY,errY, errYMC);
    printf("Correction (for second order term)= %g\n",correct);
    printf("Difference Simulaton - Theory(1st order) = %g\n",Yth-mediaY);
    printf("Valore atteso + correzione = %g \n", vattesocorr);
}