#include "TRandom3.h"
#include "TMath.h"
#include "TF1.h"
#include "TCanvas.h"
#include "TH1F.h"
#include <fstream>
#include "TNtuple.h"
#include "TH2F.h"
#include <vector>
#include <iostream>
#include <time.h>
#include "TStyle.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TBox.h"
#include "TFile.h"
#include <unistd.h>
#include <math.h>


//Funzione di distribuzione del background
double BackgroundDistribution(double tmin, double tmax)
{

	return TMath::Power(tmax - tmin, -1);
}


//Funzione di distribuzione dei tempi di decadimento dei muoni positivi
Double_t DecayTimeDistribution(double x, double tau){
	
	return TMath::Exp(-x/tau)/tau;
}

//Funzione di distribuzione dei tempi di decadimento dei muoni negativi
Double_t DecayTimeDistribution(double x, double tauPrimo, double tau){
	
	return TMath::Exp(-x/tauPrimo)/tau;
}


//Funzione di distribuzione totale (per il chi2)
Double_t TotalDistribution(double x, double alpha, double tau, double r, double tauPrimo, double norm, double tmin, double tmax)
{

	return (1 - alpha) * (DecayTimeDistribution(x, tau) + r * DecayTimeDistribution(x, tauPrimo, tau)) / norm + alpha * BackgroundDistribution(tmin, tmax);
}

















void MonteCarloConsistence2(int nExperiments, int nMaxEvents)
{


	//EVENTS SIMULATION

	double temp = 0;
	double maxTime = 0;
	double minTime = 0;
	int binNumber = 25000;
	double alphaSimulato = 0.1;
	double tauSimulato = 2.2;		
	double tauPrimoSimulato = 0.88; 
	double rSimulato = 0.8;
	double normSimulato = TMath::Exp(-0.1 / tauSimulato) - TMath::Exp(-23.67 / tauSimulato) + rSimulato * tauPrimoSimulato / tauSimulato * (TMath::Exp(-0.1 / tauPrimoSimulato) - TMath::Exp(-23.67 / tauPrimoSimulato));
	
	//TF1 function which follows my distribution:
	//-----------------------------------------------------------------------------------------------------------------------------
	TF1 *f1 = new TF1("f1","((1-[4])*(TMath::Exp(-x/[0])/[0]+[2]*TMath::Exp(-x/[1])/[0])/[3])+[4]*TMath::Power(25, -1)",0.1,23.67);
	
	f1->SetParameter(0,tauSimulato);
		f1->SetParameter(1,tauPrimoSimulato);
			f1->SetParameter(2,rSimulato);
					f1->SetParameter(3,normSimulato);
						f1->SetParameter(4,alphaSimulato);
	//-----------------------------------------------------------------------------------------------------------------------------

	//ALPHA estimate
	//-----------------------------------------------------------------------------------------------------------------------------
	double eventiBkgrd = 0;
	double alpha = 0;
	//-----------------------------------------------------------------------------------------------------------------------------
	
	
	//CHI2 ESTIMATE VARIABLES
	//-----------------------------------------------------------------------------------------------------------------------------	
	double sensitivity = 0.001;
	double Mu = 0; //Valore di aspettazione degli eventi nel i-esimo bin
	double chiSquare1 = 0;
	double tau = 2.2; //Valore dalla letteratura
	double tauPrimo = 0;
	double r = 0;
	double norm = 0;
	double x = 0;
	double beginTauHisto=0.8;
	double endTauHisto = 1.;
	double beginRHisto=0.7;
	double endRHisto = 0.9;
	
	
	double minTauPrimo = 0;
	double minR = 0;
	std::vector<double> vecStimaTau;
	std::vector<double> vecStimaR;
	std::vector<double> vecEventi;
	//-----------------------------------------------------------------------------------------------------------------------------
	
	
	
	
	//Loop over event #
	auto start1 = chrono::steady_clock::now(); 
	
	for (int nEvents = 100000; nEvents <=nMaxEvents; nEvents += 5000)
	{
		
		cout << "Events = " << nEvents << endl;
		
		TH1F *hTauPrimo = new TH1F("hTauPrimo", "#tau_{-} estimate distribution", 200, beginTauHisto, endTauHisto);
		
		TH1F *hR = new TH1F("hR", "r estimate distribution", 200, beginRHisto, endRHisto);

		//Loop over experiments
		for (int p = 0; p < nExperiments; p++)
		{

			TH1F *h1 = new TH1F("h1", "tempi", binNumber, 0, 25);
			
			for (int k = 0; k < nEvents; k++)
			{
				
					temp =f1->GetRandom();    //Fill the histogram with the wanted distribution
					h1->Fill(temp);
				
			}
			
		maxTime =(double) h1->FindLastBinAbove()/1000;
		minTime =(double) (h1->FindFirstBinAbove()-1)/1000;		
			
			h1->Rebin(100);
			//cout << "vector max: " << maxTime << " || " << "bin max: " << maxTime << endl;
			//cout << "vector min: " << minTime << " || " << "bin min: " << minTime << endl;

			eventiBkgrd = 0;
			
			for (int k = 1; k <= binNumber/100; k++)
			{

				if (h1->GetXaxis()->GetBinCenter(k) > 15)
					eventiBkgrd += h1->GetBinContent(k); //Conteggi da 15 a maxTime
			}

			//aplha estimate:from 0 to maxTime	
			eventiBkgrd = eventiBkgrd * (maxTime - minTime) / (maxTime - 15);
			alpha = eventiBkgrd / h1->GetEntries();



	
	
	
					//CHI2 ESTIMATE METHOD
//========================================================================================================================================					
//========================================================================================================================================					

					
	double sensitivity1 = 0.001;
	double Mu1 = 0; //expectation value of events in the n-bin
	double chiSquare1 = 0;
	double tau1 = 2.2; //from PDG
	double tauPrimo1 = 0; // FIRST PARAMETER TO ESTIMATE
	double r1 = 0; // SECOND PARAMETER TO ESTIMATE
	double norm1 = 0; //normalization factor
	double x1 = 0;
	double beginTauHisto1 = 0.8;
	double beginRHisto1 = 0.7;
	int count=0;
	int index_chi=0;
	double chi[4];


	
	
	
			
			
			
			double MINN=500;
			int help=0;
			index_chi=0;
			chiSquare1 = 0;
			
			while (1) 
			{
				
						
								
					for (double i =-1 ; i<2 ; i+=2){ //cycle over tauPrimo1
						
						chiSquare1 = 0;
						
						tauPrimo1 = (double)(beginTauHisto1 + i*sensitivity1);
						r1 = (double)(beginRHisto1);
						
						
						
						norm1 = TMath::Exp(-minTime / tau1) - TMath::Exp(-maxTime / tau1) + r1 * tauPrimo1 / tau1 * (TMath::Exp(-minTime / tauPrimo1) - TMath::Exp(-maxTime / tauPrimo1)); 
						
						//cout << "tauprimo: " << tauPrimo1 << " r: " << r1 << " norm: " << norm1 << " minTime: " << minTime << " maxTime: " << maxTime << endl;
						
						for (int k = 2; k <=150; k++){ //chi2 calculation starting from the first non-zero bin (the distribution is generated from the second bin)
								
								x1 = h1->GetXaxis()->GetBinCenter(k);
								Mu1 = h1->GetEntries() * h1->GetBinWidth(k) * TotalDistribution(x1, alpha, tau, r1, tauPrimo1, norm1, minTime, maxTime);
								chiSquare1 += (h1->GetBinContent(k) - Mu1) * (h1->GetBinContent(k) - Mu1) / (h1->GetBinContent(k));
							
							}
							
						help=(i+1)/2;
						
						
						chiSquare1=round(chiSquare1*1000)/1000;
						
						
						
						chi[help]=chiSquare1; //first and second component of the array filled
						
												
						}
						
						
						for (double i =-1 ; i<2 ; i+=2){ //cycle over r1 
						chiSquare1 = 0;
						tauPrimo1 = (double)(beginTauHisto1);
						r1 = (double)(beginRHisto1 + i*sensitivity1);
						
						//cout << "r: " << r1 << " , ed tauPrimo: " << tauPrimo1 <<endl;
						
						norm1 = TMath::Exp(-minTime / tau1) - TMath::Exp(-maxTime / tau1) + r1 * tauPrimo1 / tau1 * (TMath::Exp(-minTime / tauPrimo1) - TMath::Exp(-maxTime / tauPrimo1)); //soluzione analitica di norm1
						
						//cout << "norm : " << norm1 <<endl;		
												
						for (int k = 2; k <=150; k++){ 
								 
								
								x1 = h1->GetXaxis()->GetBinCenter(k);
								
								//cout<<x1<< " " << h1->GetBinWidth(k) <<endl;
								
								Mu1 = h1->GetEntries() * h1->GetBinWidth(k) * TotalDistribution(x1, alpha, tau, r1, tauPrimo1, norm1, minTime, maxTime);
								chiSquare1 += (h1->GetBinContent(k) - Mu1) * (h1->GetBinContent(k) - Mu1) / (h1->GetBinContent(k));
							}
							
						help=2+(i+1)/2;
						
						chiSquare1=round(chiSquare1*1000)/1000;
												
						
						
						
						chi[help]=chiSquare1; //third and fourth component of the array filled
						
						
						
						}
						
						
						if(chi[0]==MINN || chi[1]==MINN || chi[2]==MINN|| chi[3]==MINN){
							
							//cout<<"got the min!!"<<endl;
							//cout<< "min: "<< MINN << endl;
							break;
							
						}
						
						//cout << "=============================================="<<endl;
						//cout << chi[0] <<" "<< chi[1] <<" "<< chi[2] <<" "<< chi[3] << endl;
						//cout << "=============================================="<<endl;
						
						index_chi=5;
					for(int j=0; j< 4; j++){
						
							if (chi[j]<MINN){
							
								MINN=chi[j];
								index_chi=j;
								//cout<< "min: "<< MINN << endl;
								//cout << "=============================================="<<endl;
							
							}
					}
					
					//cout << "tauprimo: " << tauPrimo1 << " r: " << r1 <<endl;
					//cout << "=============================================="<<endl;
					//cout<< "min: "<< MINN << endl;
					//cout << "=============================================="<<endl;
					
					//cout<< "Chi min: "<< MINN << endl;
					
					
						
						if (index_chi==0){
						beginTauHisto1 -= sensitivity1;
						beginRHisto1 = beginRHisto1;
						//cout << ".................................................................." << endl;
						//cout << "decreasing tauprimo" << endl;
						} else
						
						if (index_chi==1){
						beginTauHisto1 += sensitivity1;
						beginRHisto1 = beginRHisto1;
						//cout << ".................................................................." << endl;
						//cout << "increasing tauprimo"<< endl;
						} else
						
						if (index_chi==2){
						beginTauHisto1 = beginTauHisto1;
						beginRHisto1 -= sensitivity1;
						//cout << ".................................................................." << endl;
						//cout << "decreasing r"<< endl;
						} else
						
						if (index_chi==3){
						beginTauHisto1 = beginTauHisto1;
						beginRHisto1 += sensitivity1;
						//cout << ".................................................................." << endl;
						//cout << "increasing r" << endl;
						}
						
						
					
					count ++;
					if (count > 400){break;}
				
			}	

			
			minTauPrimo = tauPrimo1;
			minR = r1;
			
			//cout<< "===================================="<<endl;
			//cout<< "minTauPrimo= "<< minTauPrimo<<" ";
			//cout<< "minR= "<< minR<<endl;
			

			hTauPrimo->Fill(minTauPrimo);
			hR->Fill(minR);


			//timevector.clear();
			
			
			delete h1;
		}

		minTauPrimo = hTauPrimo->GetMean(1);
		minR = hR->GetMean(1);
		
		cout<< "===================================="<<endl;
		cout<< "average TauPrimo= "<< minTauPrimo<<" ";
		cout<< "average R= "<< minR<<endl;

		vecStimaTau.push_back(minTauPrimo);
		vecStimaR.push_back(minR);
		vecEventi.push_back(nEvents);

		delete hTauPrimo;
		delete hR;
		
	}

	double *stimaTau = vecStimaTau.data();
	double *stimaR = vecStimaR.data();
	double *xEventi = vecEventi.data();

	TGraphErrors *g1 = new TGraphErrors(vecEventi.size(), xEventi, stimaTau);
	TGraphErrors *g2 = new TGraphErrors(vecEventi.size(), xEventi, stimaR);

	TCanvas *c1 = new TCanvas("c1", "Monte Carlo Estimates", 1000, 1000);
	c1->Divide(1, 2);

	c1->cd(1);
	g1->Draw("AP");

	c1->cd(2);
	g2->Draw("AP");

	/*TFile *f = new TFile("MonteCarloConsistenza2.root", "RECREATE");
	f->cd();
	g1->Write("ConsistenzaTau");
	g2->Write("ConsistenzaR");
	f->Close();*/

	
	auto end1 = chrono::steady_clock::now();
	
	
	cout << ".................................................................." << endl;
	cout << "Elapsed time in milliseconds for METHOD 2 : " 
		<< chrono::duration_cast<chrono::milliseconds>(end1 - start1).count()
		<< " ms" << endl;
	cout << ".................................................................." << endl;
}