#define ProcessaSimulazione_cxx
#include "ProcessaSimulazione.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>

Double_t efficiency(double x)
{
    double parEff[3];
    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% 

  
    //     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=parEff[0]/( 1.+
			 pow(81, (parEff[1]-x)/parEff[2] ) );

  return plat;
}



void ProcessaSimulazione::Loop()
{
//   In a ROOT session, you can do:
//      root> .L ProcessaSimulazione.C
//      root> ProcessaSimulazione t
//      root> t.GetEntry(12); // Fill t data members with entry number 12
//      root> t.Show();       // Show values of entry 12
//      root> t.Show(16);     // Read and show values of entry 16
//      root> t.Loop();       // Loop on all entries
//

//     This is the loop skeleton where:
//    jentry is the global entry number in the chain
//    ientry is the entry number in the current Tree
//  Note that the argument to GetEntry must be:
//    jentry for TChain::GetEntry
//    ientry for TTree::GetEntry and TBranch::GetEntry
//
//       To read only selected branches, Insert statements like:
// METHOD1:
//    fChain->SetBranchStatus("*",0);  // disable all branches
//    fChain->SetBranchStatus("branchname",1);  // activate branchname
// METHOD2: replace line
//    fChain->GetEntry(jentry);       //read all branches
//by  b_branchname->GetEntry(ientry); //read only this branch
   if (fChain == 0) return;

   // ss {
   /// inizialilizzazioni
   atTheStart();
   // }
   
   Long64_t nentries = fChain->GetEntriesFast();
   

   
   Long64_t nbytes = 0, nb = 0;
   for (Long64_t jentry=0; jentry<nentries;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      // nostro codice 
      ++nEventi;
      if (Cut(ientry) < 0) continue;
      // l'evento ha superato i tagli di selezione (implementati in Cut)
      ++nEventiPostTagli;

      // qui processiamo l'evento 
      if (nEventi%10000==0) std::cout<<"Event n. "<< nEventi<<" "<<" Eventi che passano i tagli finora "<<nEventiPostTagli<<";  Id, Evera, EMisurata "<<Id<<" "<<EVera<<" "<<EMisurata<<std::endl;
      ///


      // definisco h di riferimento per Evera
      bool EVeraOk  = EVera> Emin && EVera<Emax;
      bool EMisurataOk  = EMisurata> Emin && EMisurata<Emax;
      if ( EVeraOk )          fh_EVera->Fill(EVera);
      if ( EMisurataOk )      fh_EMisurata->Fill(EMisurata);

      if (EVeraOk && EMisurataOk)
	{
	  //
	  fh_PEC->Fill(EMisurata,EVera);
	}

   }

   /// finalizzazione della nostra logica
   atTheEnd();
   
}

// nuopstro codice da qui in poi 
void ProcessaSimulazione::atTheStart()
{
  // inizializziamo le variabili della classe che abbiamo aggiunto
  nEventi = 0;
  nEventiPostTagli = 0;
  fout = NULL;

  Emin = 16.;
  Emax = 25.;
  

  // creo istogrammi
  fh_hprova = new TH1F("hprova","hprova",200,0.,50.);

  // h EVera di riferimento
  fh_EVera = new TH1F("hEVera","EVera di riferimento", 200, Emin, Emax);
  // efficiency for EVera nel range delle misure considerate
  fh_EffInRange = new TH1F("hEffInRange","Efficienza per EVera nel range di misure 18-23 GeV", 200, Emin, Emax);
  // h EMisurata di partenza
  fh_EMisurata = new TH1F("hEMisurata","EMisurata", 200, Emin, Emax);
  // h EVeraStimata ottenuta dalla procedura di deconvoluzione
  fh_EVeraStimata = new TH1F("hEVeraStimata","Energia Stimata", 200, Emin, Emax);
  // h P0EVera ottenuta dalla procedura iterativa
  fh_P0EVera = new TH1F("hP0EVera","Prior(Energia)-iterativa", 200, Emin, Emax);
  // h EVeraStimata-new ottenuta dalla procedura di deconvoluzione
  fh_EVeraStimataN = new TH1F("hEVeraStimataN","Energia Stimata New", 200, Emin, Emax);
  // h P0EVera ottenuta dalla procedura iterativa
  fh_P0EVeraN = new TH1F("hP0EVeraN","Prior(Energia)-iterativa New", 200, Emin, Emax);

  // mappa per ProbNumero di eventi di EMisurata x ed EVera y 
  fh_PEC = new TH2F("hPEC","Evera vs EMisurata; E misurata (GeV); E vera (GeV)",200, Emin, Emax, 200, Emin, Emax); 
  // mappa per Prob di EMisurata condizionata ad EVera; x = EMisurata, y = EVera
  fh_PECNorm = new TH2F("hPECNorm","P(Emisurata | Evera); E misurata (GeV); E vera (GeV)",200, Emin, Emax, 200, Emin, Emax); 
  

}

void ProcessaSimulazione::atTheEnd()
{
  // stampiamo i contatori
  std::cout<<"**************************** Fine del loop sugli eventi *******************"<<std::endl;
  std::cout<<" N eventi letti               :  "<<nEventi<<std::endl;
  std::cout<<" N eventi cha passano i tagli :  "<<nEventiPostTagli<<std::endl;


  double NevMisInEVera_j = 0.;double NevInEVera_j = 0.;
  double eVbinc = 0.; // EVera = centro del bin y 
  TAxis* yax = (TAxis*)fh_PEC->GetYaxis();
  // Calcoliamo le probabilita' di Emisurata condizionata ad EVera 
  for (int j=1; j<=fh_PEC->GetNbinsY(); ++j)
    {
      // stabilire il contenuto totale dei bin i,j, per ogni i => N. eventi con EVera nel bin j
      // loop sui bin x e scalo per 1/Numero di eventi con EVera nel bin j
      // riempi un nuovo istogramma di probabilita' (a partira dalla distribuzione 2D di EVera e EMisurata)
      
      eVbinc = yax->GetBinCenter(j);
      NevMisInEVera_j = 0.;
      NevInEVera_j    = fh_EVera->GetBinContent(j);
      for (int im=1; im<=fh_PEC->GetNbinsX(); ++im)
	{
	  NevMisInEVera_j += fh_PEC->GetBinContent(im,j);
	}
      
      // calcolare l'efficienza totale per questa causa (j)
      double epsilon_j = 0.;
      epsilon_j = NevMisInEVera_j/NevInEVera_j;
      double epsilon_jInEMisRange = 0.;
      for (int im=1; im<=fh_PECNorm->GetNbinsX(); ++im)
	{
	  // calcolare la probabilita' che la causa j provochi l'effetto im,
	  //                               ossia che la EVera nel bin j sia misurata con una Emisurata nel bin im
	  fh_PECNorm->SetBinContent(im,j, ( fh_PEC->GetBinContent(im,j) )/NevInEVera_j);
	  // calcolare l'efficienza perche' questa causa (j) produca un effetto nel range considerato 
	  if (im>44 && im<158) epsilon_jInEMisRange +=  fh_PEC->GetBinContent(im,j);
	}
      epsilon_jInEMisRange = epsilon_jInEMisRange/NevInEVera_j;
      double effTh = efficiency(eVbinc);
      std::cout<<"Bin j = "<<j<<" EVera = "<<eVbinc<<" NevMisInEVera_j = "<<NevMisInEVera_j<<" NevInEVera_j = "<<NevInEVera_j<<" eff(16-25) = "<<epsilon_j<<" eff(18-23) = "<<epsilon_jInEMisRange<<" eff_teorica = "<<effTh<<std::endl;
      fh_EffInRange->SetBinContent(j,epsilon_jInEMisRange);


      /// inizializziamo l'istogramma della distribuzione vera stimata
      fh_P0EVera->SetBinContent(j,1./200.); // partiamo da ditribuzione uniforme 
      //fh_P0EVera->SetBinContent(j,fh_EMisurata->GetBinContent(j)/fh_EMisurata->GetEntries()); // partiamo da sitribuzione unifo
      //fh_P0EVera->SetBinContent(j,fh_EVera->GetBinContent(j)/fh_EVera->GetEntries()); // partiamo da sitribuzione uniforme 
    }
  

  
  
  // CAUSE:    200 bin di Evera
  // Effetti:  bin 45 - 157 compresi di Emisurata (Emisurata da 18 a 23 GeV)

  for (int myiter=0; myiter<1000; ++myiter)
    {

      std::cout<<" \n Unfolding::::: iteration n. "<< myiter<<std::endl;
      double pCiEj[200][200] = { 0. };
      double nCauseTot = 0.;
      for (int iC = 1; iC<=200; ++iC )
	{
	  int iCC = iC-1;
	  //std::cout<<" ic = "<<iC <<std::endl;
	  double nCause = 0.;
	  for (int jE=45; jE<158; ++jE)
	    {
 	      int jEE = jE-1;
	      pCiEj[iCC][jEE] = fh_PECNorm->GetBinContent(jE,iC) * fh_P0EVera->GetBinContent(iC);
	      if (fabs(iC-100)<5 && fabs(jE-100)<5) std::cout<<" iter, iC, jE "<<myiter<<" "<<iC<<" "<<jE<<" pCiEj "<<pCiEj[iCC][jEE] <<" "<<fh_P0EVera->GetBinContent(iC)<<std::endl;
	      // denominatore eq 3 
	      double pEj = 0.; 
	      for (int lcausa = 1; lcausa<=200; ++lcausa)
		{
		  pEj +=  fh_PECNorm->GetBinContent(jE,lcausa) * fh_P0EVera->GetBinContent(lcausa);
		}
	      if (fabs(iC-100)<5 && fabs(jE-100)<5) std::cout<<" pEj = "<<pEj<<std::endl;
	      //std::cout<<" iC, jE, pEj "<<iC<<" "<<jE<<" "<<pEj<<std::endl;
	      // Prob  della causa condizionata all'effetto 
	      if (pEj>0.0000001) pCiEj[iCC][jEE] = pCiEj[iCC][jEE] / pEj;
	       if (fabs(iC-100)<5 && fabs(jE-100)<5) std::cout<<" pCiEj = "<<pCiEj[iCC][jEE]<<std::endl;

	      // sommmiamo su tutti gli effetti per stimare il numero di cause: nCause (eq 4-)
	      nCause += fh_EMisurata->GetBinContent(jE) * pCiEj[iCC][jEE];

	    }
	  // correggiamo per l'effciienza (eq. 4)
	  if (nCause > 0.001 && fh_EffInRange->GetBinContent(iC)>0.00001)
	    {
	      nCause = nCause / double(fh_EffInRange->GetBinContent(iC));
	    }
	  else
	    {
	      std::cout<<" low Eff = "<<fh_EffInRange->GetBinContent(iC)<<" in bin "<<iC<<std::endl;
	      continue;
	    }
	  fh_EVeraStimataN->SetBinContent(iC,nCause);
	  if (fabs(iC-100)<5) std::cout<<" iC, nCause "<<iC<<" "<<nCause<<std::endl;

	  // NCauseTot
	  nCauseTot += nCause;
	  
	}
      if (nCauseTot<0.001)
	{
	  std::cout<<" nCauseTot too low "<<nCauseTot<<std::endl;
	}
      else
	{
	  // aggiorniamo la prob delle cause 
	  for (int iCn = 1; iCn<=200; ++iCn )
	    {
	      if (iCn%10==0) std::cout<<" iCn,  nCauseTot = "<<iCn<<" "<<nCauseTot<<" previous and current  prior "
		       <<fh_P0EVera->GetBinContent(iCn)<<" "<<fh_EVeraStimataN->GetBinContent(iCn)/nCauseTot<<std::endl;
	      fh_P0EVera->SetBinContent( iCn, fh_EVeraStimataN->GetBinContent(iCn)/nCauseTot );
	      // settiamo EVeraStimata a quella derivata in questa iterazione
	      fh_EVeraStimata->SetBinContent(iCn, fh_EVeraStimataN->GetBinContent(iCn));
	    }
	  std::cout<<" finito"<<std::endl; 
	}
      std::cout<<" finito1"<<std::endl; 
    }

  

  // definisco e apro il file 
  fout = new TFile("ProcSimul_output.root","RECREATE");

  
  // scriviamo gli istogrammi nel file 
  fh_hprova->Write();
  fh_EVera->Write();
  fh_EMisurata->Write();
  fh_PEC->Write();

  fh_PECNorm->Write();
  fh_EffInRange->Write();
  fh_P0EVera->Write();
  fh_EVeraStimata->Write();
  std::cout<<" finito2"<<std::endl; 

  // chiudiamo il file 
  fout->Close();  

  std::cout<<" out of Loop"<<std::endl;
}