///=> immaginiamo un processo fisico che produca fotoni monocromatici con un'energia distrubuita gaussianamente con media 60 GeV e sigma 10 GeV
// => i fotoni "di segnale" sono immersi in un fondo di fotoni con spettro esponenzialmente decrescente con media 100;
// => l'abbondanza relativa di fotoni di segnale rispetto a fotoni di fondo e' 10% 
// => un apparato sperimentale puo'rivelare fotoni a partire da energie di 15 GeV (a 15 GeV eff = 50%, l'eff sale da 10 a 90% in un intervallo di 10 GeV e raggiunge il valore asintotico di 93%)
// => la risoluzione in energia dell'apparato e' s/E = 0.1 / sqrt(E[GeV])
// 
///generiamo uno spettro di energia di fotoni reali (binnato)
// produciamo uno spettro di energia di fotoni misurati (binnato)
// determiniamo la funzione di risposta del nostro apparato:  qual e' la probabilita' che un fotone reale nel bin i-esimo, sia rivelato e misurato con energia appartenente al bin j-esimo ?
// e' possibile ottenere lo spettro reale a partire da quello misurato ?




Double_t efficiency(double *x, double *par)
{
    //     par(0) -> asymptotic efficiency
    //     par(1) -> x @ 50% of asymptotic efficiency 
    //     par(2) -> DeltaX per eff. dal 10 al 90%
    //Double_t plat=par[0]/( 1.+ 81**( (par[1]-x[0])/par[2] ) );
  Double_t plat=par[0]/( 1.+
			 pow(81, (par[1]-x[0])/par[2] ) );

  return plat;
}



// Simulatore di tree
// due classi di oggetti id=0 e 1

void EsempioSimulazione() {


    // Rapporto tra segnale e fondo
    double Srate = 0.02;
    double Srate1 = 0.05;
    // EVera media / sigma
    double EMediaS = 21.;
    double EsigmaS = 0.20;
    // EVera media / sigma
    double EMediaS1 = 19.9;
    double EsigmaS1 = 0.05;
    // Energia media fondo
    double Etau = 5.;

    double parEff[3];
    double xEff[1];
    parEff[0] = 0.93; // efficienza asintotica
    parEff[1] = 15.;  // Energia per la quale l'eff vale 50% 
    parEff[2] = 5.;  // Intervallo di Energia (attorno a parEff[1]) nel quale l'eff. va da 10 a 90% 
  
// Numero di eventi da simulare
    int N=1000000;
// Apro file
    TFile * f=new TFile("EsempioSimulazione.root","RECREATE");
// Creo il TTree
    TTree * T= new TTree("T","Simulazione");
// Le variabili del TTree
    int Id;  //  per fotone di segnale e 0 per fotone di fondo 
    double EVera, EMisurata, ThVero, ThMisurato;
    T->Branch("Id",        &Id,         "Id/I");
    T->Branch("EVera",     &EVera,      "EVera/D");
    T->Branch("EMisurata", &EMisurata,  "EMisurata/D");
    T->Branch("ThVero",    &ThVero,     "ThVero/D");
    T->Branch("ThMisurato",&ThMisurato, "ThMisurato/D");

// Innesco Random
    TRandom3 R(122342342);
// Loop principale
    for (int i=0; i<N; i++)  {

      // generiamo il processo reale
      // scelgo a caso se di classe 0 o 1
        double r=R.Uniform();
	if (r<Srate)
	  {
	    // fotone di segnale
	    Id = 1;

	    // generare l'energia secondo lo spettro del segnale
	    EVera = 0.;
	    r=R.Gaus(EMediaS,EsigmaS);
	    EVera = r;
	    ThVero = 0.;
	    
	  }
	else if (r<Srate1)
	  {
	    // fotone di segnale
	    Id = 2;

	    // generare l'energia secondo lo spettro del segnale
	    EVera = 0.;
	    r=R.Gaus(EMediaS1,EsigmaS1);
	    EVera = r;
	    ThVero = 0.;
	    
	  }
	else
	  {
	    // fotone di fondo
	    Id = 0;

	    // generare l'energia secondo lo spettro del fondo
	    EVera = 0.;
	    r=R.Exp(Etau);
	    EVera = r+14.;
	    ThVero = 0.;
	  }

	bool  isDetected = false;
	// stabilire se il fotone e' rivelato
	// il mio apparato e' efficiente o meno per il fotone prodotto ?

	xEff[0] = EVera;
	double eff = efficiency(xEff, parEff);
        r=R.Uniform();
	if (r<eff) isDetected = true;
	
	EMisurata =  -1.;
	ThMisurato = -1.;
	if (isDetected)
	 {
	   // estrarre da una gaussiana con media = EVera, sigma = 0.1 * sqrt(EVera)
	   double sigma = 0.075*sqrt(EVera); // risoluzione del mio calorimetro e.m. che misura l'energia del fotone
	   r=R.Gaus(EVera,sigma);
	   EMisurata = r;
	   ThMisurato = 0.;
	 }
	
	
// carico evento generato nel TTree
        T->Fill();
    }
// scrivo TTree nel file e chiudo il file
    T->Write();
    f->Close();

}