/************************************************************
* ROOT macro to generate artificial data file
*
* Note: function with the same name as the source code
*       can be executed in ROOT just by '.x upf_cv4_sim.cc'
*       or from terminal by 'root upf_cv4_sim.cc'
************************************************************/

void upf_cv4_sim(const char* fName = "sim.dat")
{
	const int N = 1000;
	const int nBins = 80;
	const double xMin = 0;
	const double xMax = 20;

	const double Pexp  =  0.75;	// fraction of events with exp. dist.
	                            // the rest has normal distribution
	const double tau   =  4;	// exp. dist.: f(x) = 1/tau exp(-x/tau)
	const double mu    = 10;
	const double sigma =  2;

	FILE* fOut = fopen(fName,"w");		// open file for writing
	TRandom3* rand = new TRandom3();	// initialize random number gen.
	// (https://root.cern.ch/doc/master/classTRandom.html)

	// OPTIONAL: compare the histogram of generated events with f(x):
	TH1D* h = new TH1D("hGen",";#it{x};entries",nBins,xMin,xMax);

	for (int i = 0; i < N; i++)
	{
		// fraction Pexp of events have exp. dist., the rest has normal dist.
		double r = rand->Rndm();	// random number uniformly dist. in (0,1)
		double x;
		if (r < Pexp)	// generate exponentially distributed number
			x = rand->Exp(tau);
		else	// generate normally distributed number
			x = rand->Gaus(mu,sigma);

		fprintf(fOut,"%lf\n",x);
		h->Fill(x);
	}
	fclose(fOut);	// close the output file

	// OPTIONAL: compare the histogram of generated events with f(x):
	gStyle->SetOptStat(0);	// do not draw the statistics box
	h->Draw();
	// functions (https://root.cern.ch/doc/master/classTF1.html)
	TF1* fExp  = new TF1("fExp","[0]/[1] * exp(-x/[1])",xMin,xMax);
	double norm = N * h->GetBinWidth(1);
	fExp->SetParameter(0,norm * Pexp);	// 1st par. is normalization
	fExp->SetParameter(1,tau);	// 1st parameter is tau
	fExp->Draw("same");

	TF1* fNorm = new TF1("fNorm","gaus",xMin,xMax);	// gaus = [0]*exp(-0.5*((x-[1])/[2])**2)
	fNorm->SetParameter(0,norm * (1-Pexp) / sqrt(2*TMath::Pi()) / sigma);
	fNorm->SetParameter(1,mu);
	fNorm->SetParameter(2,sigma);
	fNorm->SetLineColor(4);
	fNorm->Draw("same");

	TF1* fSum = new TF1("fSum","fExp + fNorm",xMin,xMax);
	fSum->SetLineColor(6);
	fSum->Draw("same");

	TLegend* leg = new TLegend(0.5,0.7,0.9,0.9);
	char buff[256];
	sprintf(buff,"%.0f generated events",h->GetEntries());
	leg->AddEntry(h,buff);
	sprintf(buff,"%.2f exponential distribution",Pexp);
	leg->AddEntry(fExp,buff);
	sprintf(buff,"%.2f normal distribution",1-Pexp);
	leg->AddEntry(fNorm,buff);
	leg->AddEntry(fSum,"sum of the two");
	leg->Draw();

	// Now you can try to draw the histogram of the data using 
	// the macro upf_cv4_hist.cc
}