//
// Correlazione tra variabili. L'esperienza della lente
//
#include "TRandom.h"
const double PI =3.141592653589793;

void lente(){
// posizione oggetti e distanze
    double P0k,Plk,Psk,x1k,x2k,fk;
// calcolo analitico senza tener conto della correlazione delle grandezze
    double P_0=0.0,P_l=30.0,P_s=72.0;
    double errP_0=0.05/sqrt(3);
    double errP_s=0.05/sqrt(3);
    double errP_l=1./3.;
    double errx1=sqrt(errP_l*errP_l+errP_0*errP_0);
    double errx2=sqrt(errP_l*errP_l+errP_s*errP_s);
    double x1=P_l-P_0;
    double x2=P_s-P_l;

    TH2F* hx12 = new TH2F ("hx12","Distanza immagine vs distanza oggetto",
			   100, x1-5*errx1, x1+5*errx1,
			   100, x2-5*errx2, x2+5*errx2);
   
    double fuoco=x1*x2/(x1+x2);
    double df_x1=(x2*(x1+x2)-x1*x2)/(x1+x2)/(x1+x2);
    double df_x2=(x1*(x1+x2)-x1*x2)/(x1+x2)/(x1+x2);
    double varF1=df_x1*df_x1*errx1*errx1+df_x2*df_x2*errx2*errx2;
    printf("\n Calcolo Analitico\n");
    printf("x1=%f +/- %f \n",x1,errx1);
    printf("x2=%f +/- %f \n",x2,errx2);
    printf("fuoco=%f \n",fuoco);
    printf("sigma=%f var(f)=Df_x1*var(x1)+Df_x2*var(x2))=%f\n",sqrt(varF1),varF1);

    TH1F* hf = new TH1F("hf", "Fuoco Misurato", 100, fuoco-5*sqrt(varF1), fuoco+5*sqrt(varF1) );

// numero di esperimenti
    int N=100000;
    TRandom R(123456);
    double mediaX1=0,mediaX2=0;
    double mediaX1_2=0,mediaX2_2=0;
    double mediaF=0,mediaF2=0;
    double cov=0;
// loop principale per simulare gli N esperimenti
    for (int k=0;k<N;k++) {
// simulazione delle fluttuazioni di ogni singola misura
      //P0k=R.Gaus(P_0,errP_0);
      P0k = P_0-errP_0+R.Uniform()*2.*errP_0;
      Plk=R.Gaus(P_l,errP_l);
      //Psk=R.Gaus(P_s,errP_s);
      Psk=P_s-errP_s+R.Uniform()*2.*errP_s;
// calcolo del fuoco dell'esperimento k-esimo
         x1k=Plk-P0k;
         x2k=Psk-Plk;
         fk=x1k*x2k/(x1k+x2k);

	 hx12->Fill(x1k,x2k);
	 hf->Fill(fk);
	 
 // per il calcolo di medie e varianze....
        mediaX1+=x1k;
        mediaX2+=x2k;
        mediaF+=fk;
        mediaX1_2+=x1k*x1k;
        mediaX2_2+=x2k*x2k;
        mediaF2+=fk*fk;
        cov+=(x1k-x1)*(x2k-x2);
    }
// medie e varianze delle due grandezze
    mediaX1=mediaX1/N;
    mediaX2=mediaX2/N;
    mediaX1_2=mediaX1_2/N;
    mediaX2_2=mediaX2_2/N;
    mediaF=mediaF/N;
    mediaF2=mediaF2/N;
    double varX1=(mediaX1_2-mediaX1*mediaX1);
    double varX2=(mediaX2_2-mediaX2*mediaX2);
    double varF=(mediaF2-mediaF*mediaF);
    double errX1=sqrt(varX1);
    double errX2=sqrt(varX2);
    double errF=sqrt(varF);
    cov=cov/N;
    printf("\n\n Verifica Simulazione\n");
    printf("media di X1=%f, sigma(X1)=%f e Var(x1)=%f\n",mediaX1,errX1, varX1);
    printf("media di X2=%f, sigma(X2)=%f e Var(x2)=%f\n",mediaX2,errX2, varX2);
    printf("covarianza %f\n",cov);

// calcolo varianza tenendo conto della correlazione
    double varF2=df_x1*df_x1*errx1*errx1+df_x2*df_x2*errx2*errx2+2.*df_x1*df_x2*cov;
    printf("\n\n Calcolo della varianza e sigma di F tenendo conto delle correlazione\n");
    printf("sigma=%f var(f)=Df_x1*var(x1)+Df_x2*var(x2)+2Df_x1*Df_x2*cov(X1,X2)=%f\n",sqrt(varF2),varF2);
    printf("\n\n Variabile aleatoria derivata. Risultati Simulazione");
    printf("\n media di f=%f, sigma(f)=%f e Var(f)=%f\n",mediaF,errF,varF);

    TCanvas* c1 = new TCanvas("c1","c1",1000,450);
    c1->Divide(2,1);
    c1->cd(1);
    hx12->Draw("zcol");
    c1->cd(2);
    hf->Draw();
}