TF1* expo;
vector<double> data;
TH1D* hData2;
TH1D* hPred2;
bool priorGaus=true;


//Function with unbinned logL
void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{
   expo->SetParameters(par[0],par[1]);
   f=0;
   double pred = expo->Integral(0,10);
   for(int i=0; i<data.size(); i++)
   {
      f+=-log(expo->Eval(data[i])/pred);
   }
   f+= pred-data.size()*log(pred);
   
}

//Function with binned logL
void fcn2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{
   expo->SetParameters(par[0],par[1]);
   f=0;
   for(int i=1; i<=hData2->GetNbinsX(); i++)
   {
      double dat = hData2->GetBinContent(i);
      double mc = expo->Integral(hData2->GetXaxis()->GetBinLowEdge(i),hData2->GetXaxis()->GetBinUpEdge(i));
      hPred2->SetBinContent(i,mc);
   
      if(dat>0)
	 f+=mc-dat+dat*log(dat/mc);
      else
	 f+=mc;
   }
}

//function to calculate log(likelihood*prior) for MCMC
double BayesLLH(double norm, double constant)
{
   expo->SetParameters(norm, constant);
   double pred = expo->Integral(0,10);
   double f=0;
   for(int i=1; i<=hData2->GetNbinsX(); i++)
   {
      double dat = hData2->GetBinContent(i);
      double mc = expo->Integral(hData2->GetXaxis()->GetBinLowEdge(i),hData2->GetXaxis()->GetBinUpEdge(i));
//      hPred2->SetBinContent(i,mc);
   
      if(dat>0)
	 f+=mc-dat+dat*log(dat/mc);
      else
	 f+=mc;
   }
   if(priorGaus)
      f+=0.5*pow((-0.5-constant)/0.1,2);
   return f;
}

void exponentialDecay()
{
   gRandom->SetSeed(6);

   //set up model
   expo = new TF1("expo","[0]*[1]/(exp([1]*10)-1)*exp([1]*x)",0,10);
   expo->SetParameters(20,-0.5);
   expo->SetNpx(10000);

   TF1* fixedexpo = new TF1("expo","[0]*[1]/(exp([1]*10)-1)*exp([1]*x)",0,10);
   fixedexpo->SetParameters(20,-0.5);
   fixedexpo->SetNpx(10000);
   fixedexpo->SetLineColor(kOrange);


   double Npred = expo->Integral(0,10);
   int Ndata = gRandom->Poisson(Npred);

   TH1D* hData1 = new TH1D("hData1",";x",50000,0,10);
   double binedges[6] = {0,0.5,1.5,3,6,10};
   hData2 = new TH1D("hData2",";x",5,binedges);
   TH1D* hPred = new TH1D("hPred",";x",5,binedges);

   //calculate logL before minimization (unbinned)
   double unbinnedLLH=0;
   for(int i=0; i<Ndata; i++)
   {
      double datapt = expo->GetRandom();
      data.push_back(datapt);
      hData1->Fill(datapt);
      hData2->Fill(datapt);
      
   }

   //calculate logL before minimization (binned)
   double binnedLLH=0;
   for(int i=1; i<=hData2->GetNbinsX(); i++)
   {
      double dat = hData2->GetBinContent(i);
      double mc = expo->Integral(binedges[i-1],binedges[i]);
      hPred->SetBinContent(i,mc);
      if(dat>0)
	 binnedLLH+=mc-dat+dat*log(dat/mc);
      else
	 binnedLLH+=mc;
      
   }
   hPred2 = (TH1D*)hPred->Clone("hPred2");

   //plots and minimization unbinned case
   TCanvas* c0 = new TCanvas("c0","",0,0,1200,600);
   c0->Divide(2,1);
   c0->cd(1);
   hData1->SetMaximum(11);
   hData1->Draw();
   fixedexpo->Draw("SAME");
   expo->Draw("SAME");
   std::cout << hData1->GetEntries() << "," << unbinnedLLH <<  std::endl;

   //do the magical incantations for minuit
   TMinuit *gMinuit = new TMinuit(1);
   gMinuit->SetFCN(fcn);
   Double_t arglist[10];
   Int_t ierflg = 0;

   arglist[0] = 1;
   gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
   gMinuit->DefineParameter(0,"normalization",Npred,0.001,0,0);
   gMinuit->DefineParameter(1,"constant",1,0.001,0,0);

   gMinuit->Migrad();
   
   double outpar[2], err[2];
   gMinuit->GetParameter(0,outpar[0],err[0]);
   gMinuit->GetParameter(1,outpar[1],err[1]);

   gMinuit->DefineParameter(1,"constant",-1E-12,0.001,0,0);
   gMinuit->FixParameter(1);
   gMinuit->Migrad();

   //make scan over llh in the constant parameter
   double scanpar[2];
   const int scanpts=21;
   double parval[scanpts], llhval[scanpts];
   
   double* grad;
   double f;
   int flag;
   for(int i=0; i<scanpts; i++)
   {
      parval[i]=outpar[1]-2*err[1]+i*4*err[1]/double(scanpts);
      scanpar[0]=outpar[0];
      scanpar[1]=parval[i];
      gMinuit->Eval(2,grad,f,scanpar,flag);
      llhval[i]=f;
   }

   
   c0->cd(2);			       
   TGraph* gScan = new TGraph(scanpts,parval,llhval);
   gScan->Draw("AL");
   gScan->SetLineWidth(3);
   gScan->GetXaxis()->SetTitle("constant");
   gScan->GetYaxis()->SetTitle("-ln(L)");


   //make plots and do minimization for binned case
   TCanvas* c1 =  new TCanvas("c1","",1200,0,600,600);
   hData2->Draw();
   hPred->SetLineColor(kOrange);
   hPred->Draw("SAME");
   std::cout << binnedLLH << std::endl;

   TMinuit *gMinuit2 = new TMinuit(1);
   gMinuit2->SetFCN(fcn2);

   arglist[0] = 1;
   gMinuit2->mnexcm("SET ERR", arglist ,1,ierflg);
   gMinuit2->DefineParameter(0,"normalization",Npred,0.001,0,0);
   gMinuit2->DefineParameter(1,"constant",1,0.001,0,0);

   gMinuit2->Migrad();

   hPred2->SetLineColor(kBlue);
   hPred2->Draw("SAME");


   //MCMC!
   TH2D* hMCMC = new TH2D("hMCMC",";normalization;constant",100,0,50,100,-1,0);

   double norm=20;
   double con = -0.5;
   double LLH = BayesLLH(norm,con);
   double normp, conp, LLHp;

   //Do Metropolis-Hastings algorithm
   for(int i=0; i<500000; i++)
   {
      normp=gRandom->Gaus(norm,0.2);
      conp=gRandom->Gaus(con,0.05);
      LLHp=BayesLLH(normp,conp);

      double accProb = TMath::Min(1.,TMath::Exp(LLH-LLHp));
      double fRandom = gRandom->Rndm();
      if ( fRandom <= accProb )
      {
	 norm=normp;
	 con=conp;
	 LLH=LLHp;
      }

      if(i>1000)
	 hMCMC->Fill(norm,con);
      
   }

   //Make plots from MCMC
   TCanvas* c2 = new TCanvas("c2","",0,600,1800,600);
   c2->Divide(3,1);
   c2->cd(1);
   hMCMC->Scale(1/hMCMC->Integral("width"));
   hMCMC->Draw("colz");
   
   c2->cd(2);
   TH1D* px = hMCMC->ProjectionX();
   px->Scale(1/px->Integral("width"));
   TH1D* pxc = (TH1D*)px->Clone("pxc");
   pxc->Reset();
   while(px->Integral("width")>(1-0.68))
   {
      int maxbin = px->GetMaximumBin();
      pxc->SetBinContent(maxbin,px->GetBinContent(maxbin));
      px->SetBinContent(maxbin,0);
   }
   pxc->SetFillColor(kGray);
   pxc->Draw("hist");
   px->Draw("hist SAME");
   
   
   c2->cd(3);
   TH1D* py = hMCMC->ProjectionY();
   py->Scale(1/py->Integral("width"));
   TH1D* pyc = (TH1D*)py->Clone("pyc");
   pyc->Reset();
   while(py->Integral("width")>(1-0.68))
   {
      int maxbin = py->GetMaximumBin();
      pyc->SetBinContent(maxbin,py->GetBinContent(maxbin));
      py->SetBinContent(maxbin,0);
   }
   pyc->SetFillColor(kGray);
   pyc->Draw("hist");
   py->Draw("hist SAME");

   expo->SetParameters(outpar[0],outpar[1]);
   
}