///////////// SAMPLING //// media e varianza campionaria, loro errori per diverse pdf della popolazione ////////////
#include "TRandom.h"
#include <string>

void sampling(std::string pdfpop="Gauss")
{

  std::cout<<"pdf della popolazione = "<<pdfpop<<std::endl;
  
  // n cardinalita' dei campioni
  int Nins = 5;
  // N numero di campioni   
  int Nsample = 100000;

  double gausMed = 10.;
  double gausRadVar = 3.;

  double poissMed = 10.;
  double poissVar = 10.;
  double poissRadVar = sqrt(poissVar);

  double medpop =0.;
  double sigpop =0.;
  double varpop =0.;
  
  if (pdfpop=="Gauss")
    {
      medpop = gausMed;
      sigpop = gausRadVar;
      varpop = gausRadVar*gausRadVar;
    }
  else if (pdfpop=="Poisson")
    {
      
      medpop = poissMed;
      sigpop = poissRadVar;
      varpop = poissVar;
    }
  else 
    {
      std::cout<<"pdf "<<pdfpop<<" non disponibile"<<std::endl;
    }
  
  TH1D* hpop = new TH1D("hpop","Popolazione",200, medpop-7*sigpop, medpop+7*sigpop);
  TH1D* hmediaC = new TH1D("hmediaC","Media campionaria",
			   200, medpop-7*sigpop/sqrt(double(Nins)), medpop+7*sigpop/sqrt(double(Nins)));


  double dx=7*sigpop*sigpop/sqrt(double(Nins));
  double xmin = sigpop*sigpop - dx; 
  double xmax = sigpop*sigpop + dx; 
  TH1D* hvarC   = new TH1D("hvarC"  ,"Varianza campionaria (N-1 al denom)",
			   200, xmin, xmax );
  TH1D* hvarCN  = new TH1D("hvarCN"  ,"Varianza del campione (N al den)",
			   200, xmin, xmax );
  TH1D* hvarCT  = new TH1D("hvarCT"  ,"Varianza del campione con media pop nota",
			   200, xmin, xmax );

  
  TH1D* hsqrtvarC = new TH1D("hsqrtvarC"  ,"Radice quad Varianza campionaria (N-1 al denom)",
		       200, sigpop-7*sigpop/sqrt(double(Nins)), sigpop+7*sigpop/sqrt(double(Nins)));
  TH1D* hsqrtvarCN= new TH1D("hsqrtvarCN"  ,"Radice quad della Varianza del campione (N al denom)",
		       200, sigpop-7*sigpop/sqrt(double(Nins)), sigpop+7*sigpop/sqrt(double(Nins)));
  TH1D* hsqrtvarCT= new TH1D("hsqrtvarCT"  ,"Radice quad della Varianza -media nota- (N al denom)",
		       200, sigpop-7*sigpop/sqrt(double(Nins)), sigpop+7*sigpop/sqrt(double(Nins)));
    
  TRandom R(123456);

// loop principale per simulare gli N campioni 
    for (int k=0; k<Nsample; k++) {

      double mediaCamp  = 0.;
      double mediaCamp2 = 0.;
      double varCampT = 0.;
      double x =0.;
      double xarray[1000]; 
      // loop sugli elementi del campione 
      for (int j=0; j<Nins; j++) {

	if      (pdfpop=="Gauss") x = R.Gaus(gausMed,gausRadVar);
	else if (pdfpop=="Poisson") x = R.Poisson(poissMed);
	mediaCamp += x;
	mediaCamp2+= x*x;
	hpop->Fill(x);
	xarray[j]=x;

	varCampT += (x - gausMed)*(x - gausMed);
      }
      
      mediaCamp =  mediaCamp  / Nins;
      mediaCamp2 = mediaCamp2 / Nins;

      double varCamp1 = 0.;
      for (int j=0; j<Nins; j++) {

	x = xarray[j];
	varCamp1 += (x - mediaCamp)*(x - mediaCamp);
      }
      varCamp1 = varCamp1/(Nins-1);
      double varCamp1N = (Nins-1)*varCamp1/Nins;
      
      double varCampN  = mediaCamp2 - mediaCamp*mediaCamp;
      double varCamp   = Nins*varCampN/(Nins-1);

      //varCamp  = varCamp/(Nins-1);
      varCampT = varCampT/Nins;
      //varCampN = varCampN/Nins;
      //std::cout<<"sample "<<k<<"  varCamp , varCampN  "<<varCamp<<" "<< varCampN<<" varCampT "<<varCampT<<std::endl;
      //std::cout<<"sample "<<k<<"  varCamp1, varCamp1N "<<varCamp1<<" "<< varCamp1N<<std::endl;
      
      
      hmediaC->Fill(mediaCamp);
      hvarC->Fill(varCamp);
      hvarCN->Fill(varCampN);
      hvarCT->Fill(varCampT);
      hsqrtvarC->Fill(sqrt(varCamp));
      hsqrtvarCN->Fill(sqrt(varCampN));
      hsqrtvarCT->Fill(sqrt(varCampT));
    }


    TCanvas * c = new TCanvas("c","c",900,800);
    c->Divide(3,3);
    c->cd(1);
    hpop->Draw();
    c->cd(2);
    hmediaC->Draw();
    c->cd(4);
    hvarCN->Draw();
    c->cd(5);
    hvarC->Draw();
    c->cd(6);
    hvarCT->Draw();
    c->cd(7);
    hsqrtvarC->Draw();
    c->cd(7);
    hsqrtvarCN->Draw();
    c->cd(8);
    hsqrtvarC->Draw();
    c->cd(9);
    hsqrtvarCT->Draw();

    std::cout<<"========"<<std::endl;
    std::cout<<"Popolazione con media       "<<medpop<<std::endl;
    std::cout<<"                varianza    "<<varpop<<std::endl;
    std::cout<<"                devStandard "<<sigpop<<std::endl;
    std::cout<<"========"<<std::endl;
    std::cout<<"Valore d'aspett. media campionaria: "<< medpop << " da MC = "<<hmediaC->GetMean()<<std::endl;;
    std::cout<<"Errore sulla media campionaria: "<< sigpop/sqrt(Nins)<< " da MC = "<<hmediaC->GetRMS()<<std::endl;;
    
    std::cout<<"========"<<std::endl;
    std::cout<<"Varianza del campione             [N]   da MC: "<< hvarCN->GetMean() <<" +/- "<<hvarCN->GetRMS()<< std::endl;
    std::cout<<"Varianza campionaria              [N-1] da MC: "<< hvarC->GetMean() <<" +/- "<<hvarC->GetRMS()<< std::endl;
    std::cout<<"Varianza del camp. media-pop nota [N]   da MC: "<< hvarCT->GetMean() <<" +/- "<<hvarCT->GetRMS()<< std::endl;
    

    std::cout<<"========"<<std::endl;
    std::cout<<"sqrt-Varianza del campione             [N]   da MC: "<< hsqrtvarCN->GetMean() <<" +/- "<<hsqrtvarCN->GetRMS()<< std::endl;
    std::cout<<"sqrt-Varianza campionaria              [N-1] da MC: "<< hsqrtvarC->GetMean() <<" +/- "<<hsqrtvarC->GetRMS()<< std::endl;
    std::cout<<"sqrt-Varianza del camp. media-pop nota [N]   da MC: "<< hsqrtvarCT->GetMean() <<" +/- "<<hsqrtvarCT->GetRMS()<< std::endl;

    std::cout<<"========"<<std::endl;
    std::cout<<"Errore sulla      varianza campionaria da pdf della varianza campionaria : "<< sqrt(2./(Nins-1))*hvarC->GetMean()<< std::endl;
    std::cout<<"Errore sulla sqrt-varianza campionaria da pdf della varianza campionaria : "<< hsqrtvarC->GetMean()/sqrt(2.*(Nins-1))<< std::endl;
    std::cout<<"========"<<std::endl;
    
    /*
    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();
    */
}