//
// Correlazione tra variabili. L'esperienza della lente
//
#include "TRandom.h"
const double PI =3.141592653589793;
#include "TMath.h"
// usare pi -> TMath::Pi()

void lente(){

  // usiamo il cm come unita' di lunghezza  
  
// 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_0max=0.05;//   /// errore massimo di lettura su metro graduato al mm = 0.1 cm  
    double errP_lmax=0.05;//   /// errore massimo di lettura su metro graduato al mm = 0.1 cm  
    double errP_smax=1.;//   

    double errP_0=0.05/sqrt(3.);//   /// errore statistco di lettura su metro graduato al mm = 0.1 cm  
    double errP_l=0.05/sqrt(3.);//   
    double errP_s=errP_smax/3.;      //// errore max = 1 cm /// errore stat. da distr. Gaussiana
    
    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;  /// distanza oggetto vera 
    double x2=P_s-P_l;  /// diatanza immagine vera  

    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);   /// distanza focale vera
    
    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; // var su f assummendo no correlazione

    
    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) );
    TH1F* hpo = new TH1F("hpo", "Posizione Oggetto", 100, P_0-5.*errP_0, P_0+5*errP_0 );
    TH1F* hps = new TH1F("hps", "Posizione Schermo", 100, P_s-5.*errP_s, P_s+5*errP_s );
    TH1F* hpl = new TH1F("hpl", "Posizione Lente",   100, P_l-5.*errP_l, P_l+5*errP_l );

// numero di esperimenti
    int N=10000;
    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

      /// produrre distribuzioni (istogrammi)
      P0k = P_0 -errP_0max + R.Uniform()*2.*errP_0max;
      hpo->Fill(P0k);

      Psk=R.Gaus(P_s,errP_s);
      hps->Fill(Psk);
      
      Plk = P_l -errP_lmax + R.Uniform()*2.*errP_lmax;
      hpl->Fill(Plk);
      
// 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();

    TCanvas* c2 = new TCanvas("c2","Posizioni",1000,450);
    c2->Divide(2,2);
    c2->cd(1);
    hpo->Draw();
    c2->cd(2);
    hps->Draw();
    c2->cd(3);
    hpl->Draw();
}