#include "TRandom.h"



void FrequenzaElettroni(int NmaxEsp = 5000, int numeroParticelle=10000)
{

    float const pR1e = 0.1; // P(R1|e)
    float const pR2e = 0.8; // P(R2|e)
    float const pR1a = 0.9; // P(R1|a)
    float const pR2a = 0.05; // P(R2|a)

    float  frazioneEle = 0.2; //
    
    float  Nele_atteso  =  frazioneEle*numeroParticelle;
    float  eNele_atteso =  sqrt(frazioneEle*numeroParticelle*(1.-frazioneEle)); 
    float  Nalfa_atteso =  numeroParticelle-Nele_atteso;
    float  eNalfa_atteso = eNele_atteso;
    

    // se osservo solo R1:
    // P(nR1|e) = 0.9
    // P(R1 |e) = 0.1
    // P(nR1|a) = 0.1
    // P(R1 |a) = 0.9
    //P(e|nR1) = P(nR1|e)P(e)/[ P(nR1|e)P(e) + P(nR1|a)P(a)] = 0.9 * 0.2 /( 0.9*0.2 + 0.1*0.8 ) = 0.692
    //P(a|nR1) = P(nR1|a)P(a)/[ P(nR1|e)P(e) + P(nR1|a)P(a)] = 0.1 * 0.8 /( 0.9*0.2 + 0.1*0.8 ) = 0.308
    //P(e| R1) = P( R1|e)P(e)/[ P( R1|e)P(e) + P( R1|a)P(a)] = 0.1 * 0.2 /( 0.1*0.2 + 0.9*0.8 ) = 0.027
    //P(a| R1) = P( R1|a)P(a)/[ P( R1|e)P(e) + P( R1|a)P(a)] = 0.9 * 0.8 /( 0.1*0.2 + 0.9*0.8 ) = 0.973
    double pnR1e = 1.-pR1e;
    double pnR1a = 1.-pR1a;
    double pe_nR1 = pnR1e*frazioneEle      / ( pnR1e*frazioneEle + pnR1a*(1.-frazioneEle) );
    double pa_nR1 = pnR1a*(1.-frazioneEle) / ( pnR1e*frazioneEle + pnR1a*(1.-frazioneEle) );
    double pe_R1  = pR1e*frazioneEle      / ( pR1e*frazioneEle + pR1a*(1.-frazioneEle) );
    double pa_R1  = pR1a*(1.-frazioneEle) / ( pR1e*frazioneEle + pR1a*(1.-frazioneEle) );
    std::cout<<"xBayes:  P(e|nR1) = "<<pe_nR1 <<" P(a|nR1) "<<  pa_nR1 <<std::endl;
    std::cout<<"xBayes:  P(e| R1) = "<<pe_R1  <<" P(a| R1) "<<  pa_R1  <<std::endl;


    // seosservo R1 e insieme R2: 
    double pnR2e = 1.-pR2e;
    double pnR2a = 1.-pR2a;
    // qual e' la purezza del campione di elettroni selezionando eventi nR1R2  ?
    // qual e' la purezza del campione di elettroni selezionando eventi  R1R2  ? 
    // qual e' la purezza del campione di elettroni selezionando eventi nR1nR2 ? 
    // qual e' la purezza del campione di elettroni selezionando eventi  R1nR2 ? 
    // Bayes:
    //P(e|nR1R2)   P(e|nR1nR2)   P(e|R1R2)   P(e|R1nR2)   
    float pnR1R2e   = pnR1e*pR2e;
    float pR1R2e    = pR1e*pR2e;
    float pnR1nR2e  = pnR1e*pnR2e;
    float pR1nR2e   = pR1e*pnR2e;
    //P(a|nR1R2)   P(a|nR1nR2)   P(a|R1R2)   P(a|R1nR2)  
    float pnR1R2a   = pnR1a*pR2a;
    float pR1R2a    = pR1a*pR2a;
    float pnR1nR2a  = pnR1a*pnR2a;
    float pR1nR2a   = pR1a*pnR2a;
    
    float pe_nR1R2  = pnR1R2e*frazioneEle / ( pnR1R2e *frazioneEle + pnR1R2a *(1.-frazioneEle) );	
    float pe_R1R2   = pR1R2e*frazioneEle  / ( pR1R2e  *frazioneEle + pR1R2a  *(1.-frazioneEle) );
    float pe_nR1nR2 = pnR1nR2e*frazioneEle/ ( pnR1nR2e*frazioneEle + pnR1nR2a*(1.-frazioneEle) );
    float pe_R1nR2  = pR1nR2e*frazioneEle / ( pR1nR2e *frazioneEle + pR1nR2a *(1.-frazioneEle) );	
   
    // qual e' la purezza del campione di alfa selezionando eventi nR1R2  ? 
    // qual e' la purezza del campione di alfa selezionando eventi  R1R2  ? 
    // qual e' la purezza del campione di alfa selezionando eventi nR1nR2 ? 
    // qual e' la purezza del campione di alfa selezionando eventi  R1nR2 ? 
    // Bayes:
    float pa_nR1R2  = pnR1R2a*(1.-frazioneEle) / ( pnR1R2e *frazioneEle + pnR1R2a *(1.-frazioneEle) );	
    float pa_R1R2   = pR1R2a*(1.-frazioneEle)  / ( pR1R2e  *frazioneEle + pR1R2a  *(1.-frazioneEle) );
    float pa_nR1nR2 = pnR1nR2a*(1.-frazioneEle)/ ( pnR1nR2e*frazioneEle + pnR1nR2a*(1.-frazioneEle) );
    float pa_R1nR2  = pR1nR2a*(1.-frazioneEle) / ( pR1nR2e *frazioneEle + pR1nR2a *(1.-frazioneEle) );	
    std::cout<<"xBayes:  P(e|nR1 R2) = "<<pe_nR1R2 <<" P(a|nR1 R2) "<< pa_nR1R2  <<std::endl;
    std::cout<<"xBayes:  P(e| R1 R2) = "<<pe_R1R2  <<" P(a| R1 R2) "<< pa_R1R2   <<std::endl;
    std::cout<<"xBayes:  P(e|nR1nR2) = "<<pe_nR1nR2<<" P(a|nR1 R2) "<< pa_nR1nR2 <<std::endl;
    std::cout<<"xBayes:  P(e| R1nR2) = "<<pe_R1nR2 <<" P(a| R1nR2) "<< pa_R1nR2  <<std::endl;

    float pe_nR1R2_MC;
    float pe_R1R2_MC;
    float pe_nR1nR2_MC;
    float pe_R1nR2_MC;
    float pa_nR1R2_MC;
    float pa_R1R2_MC;
    float pa_nR1nR2_MC;
    float pa_R1nR2_MC;

    
    
    TH1D* hPe_nR1 = new TH1D("hPe_nR1","hPe_nR1",100,1.,-1.);
    TH1D* hPa_nR1 = new TH1D("hPa_nR1","hPa_nR1",100,1.,-1.);
    TH1D* hPe_R1  = new TH1D("hPe_R1","hPe_R1",100,1.,-1.);
    TH1D* hPa_R1  = new TH1D("hPa_R1","hPa_R1",100,1.,-1.);

    TH1F* hfe_fromR1  = new TH1F("hfe_fromR1","fe_fromR1",100,1.,-1.);
    TH1F* hfe_fromR2  = new TH1F("hfe_fromR2","fe_fromR2",100,1.,-1.);
    TH1F* hfe_fromR12 = new TH1F("hfe_fromR12","fe_fromR12",100,1.,-1.);
    TH2D* hfe_fromR1_fromR2 = new TH2D("hfe_fromR1_fromR2","fe_fromR2 vs fe_fromR1",100,1.,-1.,100,1.,-1.);

    TH1F* hPe_R1R2   = new TH1F("hPe_R1R2"    ,    "Pe_R1R2"    ,    100, 1.,-1.);
    TH1F* hPe_nR1R2  = new TH1F("hPe_nR1R2"   ,    "Pe_nR1R2"   ,    100, 1.,-1.);
    TH1F* hPe_R1nR2  = new TH1F("hPe_R1nR2"   ,    "Pe_R1nR2"   ,    100, 1.,-1.);
    TH1F* hPe_nR1nR2 = new TH1F("hPe_nR1nR2"  ,    "Pe_nR1nR2"  ,    100, 1.,-1.);

    TH1F* hPa_R1R2   = new TH1F("hPa_R1R2"    ,    "Pa_R1R2"    ,    100, 1.,-1.);
    TH1F* hPa_nR1R2  = new TH1F("hPa_nR1R2"   ,    "Pa_nR1R2"   ,    100, 1.,-1.);
    TH1F* hPa_R1nR2  = new TH1F("hPa_R1nR2"   ,    "Pa_R1nR2"   ,    100, 1.,-1.);
    TH1F* hPa_nR1nR2 = new TH1F("hPa_nR1nR2"  ,    "Pa_nR1nR2"  ,    100, 1.,-1.);

    int Nele = 0;       // Numero totale di elettroni 
    int Nalfa = 0;      // Numero totale di alfa 
    int NeleR1 = 0;     // Numero di e con segnale in R1 
    int NeleR2 = 0;     // Numero di e con segnale in R2
    int NalfaR1 = 0;    // Numero di alfa con segnale in R1
    int NalfaR2 = 0;    // Numero di alfa con segnale in R2

    TH1D* hNele    = new TH1D("hNele" ,  "Nele" ,  100,  Nele_atteso-5.*eNele_atteso, Nele_atteso+5.*eNele_atteso); // Numero totale di elettroni 
    TH1D* hNeleR1  = new TH1D("hNeleR1" ,"NeleR1" ,201,  99.5,300.5); // Numero di e con segnale in R1 
    TH1D* hNeleR2  = new TH1D("hNeleR2" ,"NeleR2" ,601,  1299.5,1900.5); // Numero di e con segnale in R2
    TH1D* hNalfa   = new TH1D("hNalfa",  "Nalfa",  601,  7699.5,8300.5); // Numero totale di alfa 
    TH1D* hNalfaR1 = new TH1D("hNalfaR1","NalfaR1",501,  6899.5,7400.5); // Numero di alfa con segnale in R1
    TH1D* hNalfaR2 = new TH1D("hNalfaR2","NalfaR2",201,  299.5,500.5); // Numero di alfa con segnale in R2

    int NR1 = 0;        // Numero di segnale in R1
    int NR2 = 0;        // Numero di segnale in R2
    int NR1R2 = 0;      // Numero di particelle che danno segnale in R1 e      in R2  
    int NR1nR2 = 0;     // Numero di particelle che danno segnale in R1 e non  in R2  
    int NnR1R2 = 0;     // Numero di particelle che danno segnale in R2 e non  in R1  
    int NnR1nR2 = 0;    // Numero di particelle che non danno segnale in R1 e  in R2
    TH1D* hNR1       = new TH1D("hNR1","NR1", 100, 100.,-100.);
    TH1D* hNR2       = new TH1D("hNR2","NR2", 100, 100.,-100.);
    TH1D* hNR1R2     = new TH1D("hNR1R2","NR1R2", 100, 100.,-100.);
    TH1D* hNR1nR2    = new TH1D("hNR1nR2","NR1nR2", 100, 100.,-100.);
    TH1D* hNnR1R2    = new TH1D("hNnR1R2","NnR1R2", 100, 100.,-100.);
    TH1D* hNnR1nR2   = new TH1D("hNnR1nR2","NnR1nR2", 100, 100.,-100.);

    int succR1=0;       // variabile temporanea 
    int succR2=0;       // variabile temporanea 
    
    float x;
    TRandom3 R(1556457);


    int Nesp = 1;
    while (Nesp <= NmaxEsp)
      {
	Nele = 0;       // Numero totale di elettroni 
	Nalfa = 0;      // Numero totale di alfa 
	NeleR1 = 0;     // Numero di e con segnale in R1 
	NeleR2 = 0;     // Numero di e con segnale in R2
	NalfaR1 = 0;    // Numero di alfa con segnale in R1
	NalfaR2 = 0;    // Numero di alfa con segnale in R2
	NR1 = 0;        // Numero di segnale in R1
	NR2 = 0;        // Numero di segnale in R2
	NR1R2 = 0;      // Numero di particelle che danno segnale in R1 e      in R2  
	NR1nR2 = 0;     // Numero di particelle che danno segnale in R1 e non  in R2  
	NnR1R2 = 0;     // Numero di particelle che danno segnale in R2 e non  in R1  
	NnR1nR2 = 0;    // Numero di particelle che non danno segnale in R1 e  in R2  
	succR1=0;       // variabile temporanea 
	succR2=0;       // variabile temporanea 


	int Ne_R1R2 = 0;
	int Ne_nR1R2 = 0;
	int Ne_nR1nR2 = 0;
	int Ne_R1nR2 = 0;

	int Na_R1R2 = 0;
	int Na_nR1R2 = 0;
	int Na_nR1nR2 = 0;
	int Na_R1nR2 = 0;
	bool is_ele=false;
	for (int i=0; i<numeroParticelle; ++i)
	  {
	    is_ele=false;
	    // simulo la particella i-esima 
	    x=R.Uniform();  
	    succR1 = 0;
	    succR2 = 0;
	    if (x<frazioneEle) // simulo la sua natura: e oppure alfa
	      {
		is_ele = true;
		// elettrone
		Nele++;
		x =  R.Uniform();
		if (x<pR1e)   // simulo la risposta del rivelatore 1 a elettrone 
		  {
		    // R1 per elettrone
		    NeleR1++;
		    NR1++;
		    succR1 = 1;
		  }
		x =  R.Uniform();
		if (x<pR2e)   // simulo la risposta del rivelatore 2 a elettrone con una nuova estrazione perche'  R1 e R2 sono indipendenti
		  {
		    // R2 per elettrone
		    NeleR2++;
		    NR2++;
		    succR2 = 1;
		  }
	      }
	    else
	      {
		// alfa
		Nalfa++;
		x =  R.Uniform(); 
		if (x<pR1a)   // simulo la risposta del rivelatore 1 a alfa
		  {
		    // R1 per alfa
		    NalfaR1++;
		    NR1++;
		    succR1 = 1;
		  }
		x =  R.Uniform();
		if (x<pR2a)   // simulo la risposta del rivelatore 2 a alfa con una nuova estrazione perche' la risposta di R1 e di R2 sono indipendenti
		  {
		    // R2 per alfa
		    NalfaR2++;
		    NR2++;
		    succR2 = 1;
		  }
	      }
	    
	    if (succR1+succR2 == 2)
	      {
		NR1R2++;
		if (is_ele) Ne_R1R2++;
		else        Na_R1R2++;
	      }
	    else if (succR1==1)
	      {
		NR1nR2++;
		if (is_ele) Ne_R1nR2++;
		else        Na_R1nR2++;
	      }
	    else if (succR2==1)
	      {
		NnR1R2++;
		//if (!is_ele) std::cout<<" Ev n "<<Nesp<<"  part n. "<<i<<" is_ele = "<<is_ele<<std::endl;
		if (is_ele) Ne_nR1R2++;
		else        Na_nR1R2++;
	      }
	    else
	      {
		NnR1nR2++;
		if (is_ele) Ne_nR1nR2++;
		else        Na_nR1nR2++;
	      }
	  }


	//if (Nesp%10000==0 || Nesp==NmaxEsp) {
	if (Nesp==1780) {
	  std::cout<<"\n\n ---------------------Esperimento n "<<Nesp<<" ---------------------------------"<<std::endl;
	  std::cout<<"Nparticelle          "<<numeroParticelle<<" \n";
	  std::cout<<"Frazione di e-       "<<frazioneEle<<" \n"; 
	  std::cout<<"Numero di e-  attesi "<<Nele_atteso<<" +/- "<<eNele_atteso <<" Numero di alfa  attesi "<< Nalfa_atteso<<" +/- "<<eNalfa_atteso<<" \n" <<std::endl;
	  std::cout<<"N ele / alfa = " << Nele<<" "<<Nalfa<<std::endl;
	  std::cout<<"N ele  con segnale in R1 " << NeleR1 <<"          N  ele con segnale in R2 " << NeleR2<<std::endl;
	  std::cout<<"N alfa con segnale in R1 " << NalfaR1<<"          N alfa con segnale in R2 " << NalfaR2<<std::endl<<std::endl;
	  std::cout<<"N segnali in  R1         "<<NR1      <<"          R2                       " << NR2<<std::endl;
	  std::cout<<"N segnali in  R1 &&  R2 "<<NR1R2<<std::endl;
	  std::cout<<"N segnali in  R1 && !R2 "<<NR1nR2<<" somma "<<NR1R2+NR1nR2<<" = NR1"<<std::endl;
	  std::cout<<"N segnali in  R2 && !R1 "<<NnR1R2<<" somma "<<NR1R2+NnR1R2<<" = NR2"<<std::endl;
	  std::cout<<"N segnali in !R1 && !R2 "<<NnR1nR2<<" somma "<<NR1R2+NR1nR2+NnR1R2+NnR1nR2<<std::endl;

	  std::cout<<"xMC:  P(e|nR1) = "<<double(Nele-NeleR1)/double(NnR1R2+NnR1nR2) <<" P(a|nR1) "<< double(Nalfa-NalfaR1)/double(NnR1R2+NnR1nR2) <<std::endl;
	  std::cout<<"xMC:  P(e| R1) = "<<double(NeleR1)/double(NR1) <<" P(a|R1) "<< double(NalfaR1)/double(NR1)<<std::endl;;
	  std::cout<<" ---------------------Esperimento n "<<Nesp<<" ---------------------------------"<<std::endl<<std::endl;


	}
	
	hNele   ->Fill(Nele);
	hNeleR1 ->Fill(NeleR1);
	hNeleR2 ->Fill(NeleR2);
	hNalfa  ->Fill(Nalfa);
	hNalfaR1->Fill(NalfaR1);
	hNalfaR2->Fill(NalfaR2);

	hNR1       ->Fill(NR1   );
	hNR2       ->Fill(NR2   );
	hNR1R2     ->Fill(NR1R2 );
	hNR1nR2    ->Fill(NR1nR2);
	hNnR1R2    ->Fill(NnR1R2);
	hNnR1nR2   ->Fill(NnR1nR2);

	hPe_nR1    ->Fill(double(Nele-NeleR1)/double(NnR1R2+NnR1nR2) );
	hPa_nR1	   ->Fill(double(Nalfa-NalfaR1)/double(NnR1R2+NnR1nR2) );
	hPe_R1 	   ->Fill(double(NeleR1)/double(NR1) );
	hPa_R1 	   ->Fill(double(NalfaR1)/double(NR1));


	//DA MC
	pe_nR1R2_MC  = double(Ne_nR1R2)/double(NnR1R2);
	pe_R1R2_MC   = double(Ne_R1R2) /double(NR1R2);
	pe_nR1nR2_MC = double(Ne_nR1nR2)/double(NnR1nR2);
	pe_R1nR2_MC  = double(Ne_R1nR2)/double(NR1nR2);
	//DA MC
	pa_nR1R2_MC  = double(Na_nR1R2)/double(NnR1R2);
	pa_R1R2_MC   = double(Na_R1R2) /double(NR1R2);
	pa_nR1nR2_MC = double(Na_nR1nR2)/double(NnR1nR2);
	pa_R1nR2_MC  = double(Na_R1nR2)/double(NR1nR2);

	hPe_R1R2    ->Fill(pe_R1R2_MC    );
	hPe_nR1R2   ->Fill(pe_nR1R2_MC   );
	hPe_R1nR2   ->Fill(pe_R1nR2_MC   );
	hPe_nR1nR2  ->Fill(pe_nR1nR2_MC  );
	
	hPa_R1R2    ->Fill(pa_R1R2_MC    );
	hPa_nR1R2   ->Fill(pa_nR1R2_MC   );
	hPa_R1nR2   ->Fill(pa_R1nR2_MC   );
	hPa_nR1nR2  ->Fill(pa_nR1nR2_MC  );
	
	//// stimo frazioneElettroni
	//// nota NR1 = N x fe x P(R1|e) - N x fe x P(R1|a) + N x P(R1|a) => N x fe (  P(R1|e) - P(R1|a) ) = NR1 - N P(R1|a)
	float fe_fromR1 = ( (double(NR1)-numeroParticelle*pR1a) /double(numeroParticelle) ) /( pR1e - pR1a );
	hfe_fromR1->Fill(fe_fromR1);

	//// stimo frazioneElettroni
	//// nota NR2 = N x fe x P(R2|e) - N x fe x P(R2|a) + N x P(R2|a) => N x fe (  P(R2|e) - P(R2|a) ) = NR2 - N P(R2|a) 
	float fe_fromR2 = ( (double(NR2)-numeroParticelle*pR2a) /double(numeroParticelle) ) /( pR2e - pR2a );
	hfe_fromR2->Fill(fe_fromR2);

	// media delle due stime (derivate da NR1 e da NR2)
	hfe_fromR12->Fill((fe_fromR1+fe_fromR2)/2.);
	// osservo le correlazioni 
	hfe_fromR1_fromR2->Fill(fe_fromR2,fe_fromR1);
	

	Nesp++;
      }


    TCanvas* cT = new TCanvas("cT","Truth",800,800);
    cT->Divide(2,3);
    cT->cd(1);
    hNele->Draw();
    cT->cd(2);
    hNalfa->Draw();
    cT->cd(3);
    hNeleR1->Draw();
    cT->cd(4);
    hNalfaR1->Draw();
    cT->cd(5);
    hNeleR2->Draw();
    cT->cd(6);
    hNalfaR2->Draw();
 
    TCanvas* cR = new TCanvas("cR","Detectors",800,800);
    cR->Divide(2,3);
    cR->cd(1);
    hNR1->Draw();
    cR->cd(2);
    hNR2->Draw();
    cR->cd(3);
    hNR1R2->Draw();
    cR->cd(4);
    hNR1nR2->Draw();
    cR->cd(5);
    hNnR1R2->Draw();
    cR->cd(6);
    hNnR1nR2->Draw();


    TCanvas* cB = new TCanvas("cB","Bayes",800,600);
    cB->Divide(2,2);
    cB->cd(1);
    hPe_nR1->Draw();
    cB->cd(2);
    hPa_nR1 ->Draw();
    cB->cd(3);
    hPe_R1 ->Draw();
    cB->cd(4);
    hPa_R1 ->Draw();

    TCanvas* cF = new TCanvas("cF","Frazione Elettroni",800,600);
    cF->Divide(2,2);
    cF->cd(1);
    hfe_fromR1->Fit("gaus");
    cF->cd(2);
    hfe_fromR2->Fit("gaus");
    cF->cd(3);
    hfe_fromR12->Fit("gaus");
    cF->cd(4);
    hfe_fromR1_fromR2->Draw("zcol");

    TCanvas* cPe = new TCanvas("cPe","Electron Purity in samples",800,600);
    cPe->Divide(2,2);
    cPe->cd(1);
    hPe_R1R2->Draw();
    cPe->cd(2);
    hPe_nR1R2->Draw();
    cPe->cd(3);
    hPe_R1nR2->Draw();
    cPe->cd(4);
    hPe_nR1nR2->Draw();
    TCanvas* cPa = new TCanvas("cPa","Alfa Purity in samples",800,600);
    cPa->Divide(2,2);
    cPa->cd(1);
    hPa_R1R2->Draw();
    cPa->cd(2);
    hPa_nR1R2->Draw();
    cPa->cd(3);
    hPa_R1nR2->Draw();
    cPa->cd(4);
    hPa_nR1nR2->Draw();

}