//This macro produce manually the NLL and fit the parameters

//a dataset is created using the pdf mu*N_sig*S(x)+N_bkg*B(X),
//where the signal is a gaussian and the background is a uniform distribution

//From the fit on the dataset mu and mean(Higgs mass) are estrapolated

//the 2 NLL projections and a countor plot are plotted and saved

///////////////////////////////////
//TO RUN THE MACRO
//1. create directory model_2d
//2. run root -b contour_plots.C
////////////////////////////////

#include "RooGaussian.h"
#include "RooUniform.h"
#include "RooAbsPdf.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "TAxis.h"
#include "TH1.h"
#include "TRandom3.h"
#include "RooAddPdf.h"
#include "stdio.h"

using namespace RooFit ;
using namespace std;


void contour_plots(){

  const int dim = 5;

  Double_t alpha[dim]={1,10,20,50,100};
 
  TString title[dim]={"alpha=1","alpha=10","alpha=20","alpha=50","alpha=100"};
 

  TRandom3 a,b;
 
  Double_t mu_th=1; //theoretical value of mu
  Double_t Higgs_th=125;

 Double_t mu_fit[dim];
 Double_t mu_err[dim];
 Double_t mass_fit[dim];
 Double_t mass_err[dim]; 

  TCanvas *c = new TCanvas("c","NLL",1100,800);
  c->Divide(3,1);
  
  for(int i=0;i<dim;i++){ //loop alpha
    

      //Declaring  variables
      RooRealVar x("x","x",115,130) ;
      RooRealVar mean("mean","Higgs mass",Higgs_th,120,130);
      RooRealVar sigma("sigma","width of gaussian",2.5) ;
      RooRealVar mu("mu","signal strenght",mu_th,0,5);
      
      
      // Build gaussian p.d.f in terms of x,mean and sigma
      RooGaussian sig("sig","Signal shape",x,mean,sigma); 
      RooUniform bkg("bkg","Background shape",x);   
      
      //Setting Nsig, Nbkg,mu and n_s
      Double_t N_s=(alpha[i])*(a.Poisson(17.35));
      Double_t N_b=(alpha[i])*(b.Poisson(15.78));
      RooRealVar Nsig("Nsig","signal events",N_s);
      RooRealVar Nbkg("Nbkg","bkg events",N_b);
      RooFormulaVar n_s("n_s","@0*@1",RooArgList(Nsig,mu));
      
      
      //Creating model
      RooAddPdf model("model","model",RooArgList(sig,bkg),RooArgList(n_s,Nbkg));
      
      //Creating dataset
      //if number events not specified creates ev=sig+bkg
      RooDataSet *data = model.generate(x);
      
      //Building and drawing NLL
      RooAbsReal* nll = model.createNLL(*data) ;
      
      RooPlot* frame1 = mu.frame() ;
      RooPlot* frame2 = mean.frame() ;
      nll->plotOn(frame1,ShiftToZero()) ;
      nll->plotOn(frame2,ShiftToZero()) ;
      
      RooMinuit m(*nll);
      m.migrad(); //MIGRAD minimization
      
      mu_fit[i] = mu.getVal();
      mu_err[i] = mu.getError();
      mass_fit[i] = mean.getVal();
      mass_err[i] = mean.getError();

      RooPlot*   frame3 = m.contour(mu,mean,1,1.96);
      frame3->GetXaxis()->SetLimits((mu_fit[i]-(3*mu_err[i])),(mu_fit[i]+(3*mu_err[i]))); //setting axis limits
      frame3->SetMinimum(mass_fit[i]-(3*mass_err[i])); //setting axis limits
      frame3->SetMaximum(mass_fit[i]+(3*mass_err[i])); //setting axis limits
      frame3->SetTitle(title[i]);

      c->cd(1);frame1->Draw();
      c->cd(2);frame2->Draw();
      c->cd(3);frame3->Draw();
      c->Print("model_2d/NLL_"+title[i]+".eps");
  }//end loop alpha

  //printing fitted values  
  for(int i=0;i<dim;i++){
    cout<<""<<endl;
    cout<<"alpha: "<<alpha[i]<<endl;
    cout<<"                   mu fitted: "<<mu_fit[i]<<" +/- "<<mu_err[i]<<endl;     
    cout<<"           Higgs mass fitted: "<<mass_fit[i]<<" +/- "<<mass_err[i]<<endl;
  }
} //end void