void MediaPesata(int N, double *x, double *sigmax, double &xmed, double &sigma) {
  //void MediaPesata(int N, double *x, double *sigmax) {

  double sx2=0;
  double ssm1=0;

  for (int i=0; i<N; i++) {
    sx2 += x[i]/sigmax[i]/sigmax[i];
    ssm1 += 1./sigmax[i]/sigmax[i];
  }
  xmed = sx2/ssm1;
  sigma = TMath::Sqrt(1/ssm1);
  
}



void macro6()
{
  double I = 30.00; // mA
  double V = 3000; // mA

  double sigmarI = 0.006;
  double sigmarV = 0.6;

  double sigmakI = 0.005;
  double sigmakV = 0.005;

  int nMis = 10000;
  int nStrum = 1;

  double R[10000];
  double sigmaR[10000];

  TH1F *hR = new TH1F("hR","",200,95,105);

  for (int j=0; j<nStrum; j++) {

    double kI = gRandom->Gaus(1,sigmakI);
    double kV = gRandom->Gaus(1,sigmakV);

    for (int n=0; n<nMis; n++) {
      
      double Ii = gRandom->Gaus(I,sigmarI)*kI;
      double Vi = gRandom->Gaus(V,sigmarV)*kV;
      double sigmaI = TMath::Sqrt(sigmarI*sigmarI + sigmakI*Ii*sigmakI*Ii);
      double sigmaV = TMath::Sqrt(sigmarV*sigmarV + sigmakV*Vi*sigmakV*Vi);


      int iMis = n+j*nMis;

      R[iMis] = Vi/Ii;
      
      sigmaR[iMis] = R[iMis]*TMath::Sqrt(sigmaI*sigmaI/Ii/Ii + sigmaV*sigmaV/Vi/Vi);
      
      hR->Fill(R[iMis]);
    }

  }

  int Nmisure = nMis*nStrum;

  double Rpes, sigmaRpes;
  
  MediaPesata(Nmisure, R, sigmaR, Rpes, sigmaRpes);
  //MediaPesata(Nmisure, R, sigmaR);


  printf("\nmedia pesata: %f +- %f\n",Rpes,sigmaRpes);

  hR->Draw();

  double Rmean = hR->GetMean();
  double sigmaRmean = hR->GetRMS();

  printf("\nR = %f +- %f \n",Rmean,sigmaRmean);


}