#if !defined(__CINT__) || defined(__MAKECINT__)
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#include <signal.h>
#include <unistd.h>
#include <stdio.h>
#include <iostream>
#include <algorithm>
#include <fstream>
#include <TROOT.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TH1.h>
#include <TF1.h>
#include <TGraph.h>
#include <TGraphErrors.h>
#include <TRandom.h>
#include <TRandom1.h>
#include <TSystem.h>
#include <TApplication.h>
#include <TCanvas.h>
#include <TCanvasImp.h>
#include <TH1.h>
#include <TFile.h>
#include <TInterpreter.h>
#include <TBenchmark.h>
#include <TNtuple.h>
#include <TMatrixD.h>
#endif

/*
  esempio: maglia con V=5 V, R1=20 Ohm, R2 = 30 Ohm
           -> previsto: V1=2V , V2=3V

    +----R1----+
    |          |   
    V         R2   
    |          |
    +----------+


*/

void macro8(UInt_t seed=1234) {

  TRandom1 *rand=new TRandom1(seed);

  // valori misurati
  double V=5.003;
  double R1=20.5;
  double R2=29.7;

  double V1=2.042;
  double V2=2.960;

  // dati strumentali
  double DeltaDgtR=0.2;   // 2 digits = 0.2 ohm
  double DeltakR=0.01;    // scala: 1%
  
  double DeltaDgtV=0.001; // 1 digit = 1mV
  double DeltakV=0.005;   // 0.5%


  // qualche previsione calcolata
  double RS=R1+R2;
  double spercR1=R1*DeltakR*0.58;
  double spercR2=R2*DeltakR*0.58;
  double sdgtR1=DeltaDgtR*0.58;
  double sdgtR2=DeltaDgtR*0.58;
  double sR1=TMath::Sqrt(spercR1*spercR1+sdgtR1*sdgtR1);
  double sR2=TMath::Sqrt(spercR2*spercR2+sdgtR2*sdgtR2);
  double sRS=TMath::Sqrt(sdgtR1*sdgtR1+sdgtR1*sdgtR1+(RS*DeltakR*0.58)*(RS*DeltakR*0.58));
  printf("R1 = %f +- %f\n",R1,sR1);
  printf("R2 = %f +- %f\n",R2,sR2);
  printf("RS = %f +- %f\n",RS,sRS);

  TFile *fout=TFile::Open("rete.root","RECREATE");
  TNtuple *nt=new TNtuple("nt","simul. rete","strum:V:R1:R2:RS:V1:V2:dV1:dV");

  TH1F *hR1=new TH1F("hR1","R1",100,19,21);
  TH1F *hR2=new TH1F("hR2","R2",100,29,31);
  TH1F *hRS=new TH1F("hRS","RS",100,48,52);

  TH1F *hV=new TH1F("hV","V",100,4.95,5.05);
  TH1F *hV1=new TH1F("hV1","V1",100,1.95,2.05);
  TH1F *hV2=new TH1F("hV2","V2",100,2.95,3.05);
  TH1F *hdV1=new TH1F("hdV1","deltaV1",100,-0.05,0.05);
  TH1F *hdV=new TH1F("hdV","deltaV",100,-0.005,0.005);

  int nStrum=1000;
  int nMis=1000;

  // loop sugli strumenti  
  for (int iStrum=0; iStrum<nStrum; iStrum++) {
    if (!(iStrum%100)) printf("strumento n. %d ...\n",iStrum);

    double kV=rand->Uniform(1-DeltakV,1+DeltakV);
    double kR=rand->Uniform(1-DeltakR,1+DeltakR);

    // loop sulle misure con un dato strumento
    for (int iMis=0; iMis<nMis; iMis++) {

      double v = V*kV + rand->Uniform(-DeltaDgtV,+DeltaDgtV);      
      double r1 = R1*kR + rand->Uniform(-DeltaDgtR,+DeltaDgtR);
      double r2 = R2*kR + rand->Uniform(-DeltaDgtR,+DeltaDgtR);

      double v1 = V1*kV + rand->Uniform(-DeltaDgtV,+DeltaDgtV);      
      double v2 = V2*kV + rand->Uniform(-DeltaDgtV,+DeltaDgtV);      

      double rs=r1+r2;
      double v1r=v*r1/(r1+r2); 
      
      // confronto V1 previsto da V,R1,R2 e V1 misurato
      double dv1=v1-v1r;
      
      // confronto col V1+V2 misurato
      double dv=v-(v1+v2);


      hR1->Fill(r1);
      hR2->Fill(r2);
      hRS->Fill(rs);
      hV->Fill(v);
      hV1->Fill(v1);
      hV2->Fill(v2);
      hdV1->Fill(dv1);
      hdV->Fill(dv);

      nt->Fill(iStrum,v,r1,r2,rs,v1,v2,dv1,dv);
    
    } // fine loop sulle misure
  
  } // end loop sugli strumenti


  double dV1mean=hdV1->GetMean();
  double dV1rms=hdV1->GetRMS();
  printf("\nDeltaV1 atteso-misurato:  %f  +-  %f\n",dV1mean,dV1rms);
  printf("Compatibilita': %f\n",TMath::Abs(dV1mean)/dV1rms);

  double dVmean=hdV->GetMean();
  double dVrms=hdV->GetRMS();
  printf("\nDeltaV atteso-misurato:  %f  +-  %f\n",dVmean,dVrms);
  printf("Compatibilita': %f\n",TMath::Abs(dVmean)/dVrms);

  fout->cd();
  fout->Write();
  fout->Close();
}