#include "Riostream.h"
#include "TFile.h"
#include "TTree.h"
#include "TCanvas.h"
#include "TH1D.h"
#include "TF1.h"
#include "TLegend.h"
#include "TStyle.h"
#include "TMatrixDSym.h"
#include "TFitResult.h"
#include "TGraph.h"
#include "TMultiGraph.h"

using namespace std;

const double min_de = -0.15;
const double max_de =  0.15;
const int bins = 40;

double pdf_tot(double* x, double* par){
  
  //the function variable
  double B_de = x[0];
  
  //the array of parameters
  double N_sig = par[0];
  double mu    = par[1];
  double sigma = par[2];
  double N_misid = par[3];
  double mu_mis = par[4];
  double N_bkg = par[5];
  double slope = par[6];
  
  //signal pdf, normalised to 1
  double pdf_sig = exp(-(B_de - mu)*(B_de - mu)/2./sigma/sigma)/sqrt(2.*TMath::Pi())/sigma;
  //misid pdf, normalised to 1
  double pdf_misid = exp(-(B_de - mu_mis)*(B_de - mu_mis)/2./sigma/sigma)/sqrt(2.*TMath::Pi())/sigma;
  //backg pdf, normalised to 1
  double pdf_bkg = (1.+slope*B_de)/((max_de-min_de)+slope/2.*(max_de*max_de-min_de*min_de));
  
  double delta = (max_de-min_de)/bins;
  //total pdf
  double pdf_tot = (N_sig*pdf_sig + N_misid*pdf_misid + N_bkg*pdf_bkg)*delta;
  
  return pdf_tot;
}

void fitDeltaE_reparametrised(){
    
    //define an histogram to look at deltaE distribution
    TH1D* h_data = new TH1D("h_data",";#DeltaE [GeV]; Entries",bins,min_de,max_de);

    //open file and take the tree
    TFile* file = TFile::Open("data.root");
    TTree* tree = (TTree*) file->Get("treeData");
    
    int tot_entries = tree->GetEntries();
    cout << "Total entries in the tree: " << tot_entries << endl;
    
    //link the variables with tree banches
    double B_de;
    double bkg_killer;
    tree->SetBranchAddress("B_de",&B_de);
    tree->SetBranchAddress("bkg_killer",&bkg_killer);
    
    //loop over the entries
    for(int iEntry; iEntry<tot_entries; ++iEntry){
        
        tree->GetEntry(iEntry);
        
        //skip all candidates below the optimal cut point
        if(bkg_killer<0.92) continue;
        
        //fill the histograms
        h_data->Fill(B_de);
    }
    
    
    //Let's define the PDF for the fit, using TF1
    //https://root.cern.ch/doc/master/classTF1.html
    
    //The total function that describes our observed distribution
    TF1* pdf = new TF1("pdf",pdf_tot,min_de,max_de,7);
    
    //signal gauss, normalisation constant
    pdf->SetParName  (0,  "N_{sig}");
    pdf->SetParameter(0,  860);
    //signal gauss, mean fixed
    pdf->SetParName  (1,  "#mu_{sig}");
    pdf->FixParameter(1,   0.);
    //signal gauss, std dev fixed
    pdf->SetParName  (2,  "#sigma_{sig}");
    pdf->SetParameter(2,0.015);
    //mis-id gauss, normalisation constant
    pdf->SetParName  (3,  "N_{misid}");
    pdf->SetParameter(3,   220);
    //mis-id gauss, mean fixed
    pdf->SetParName  (4,  "#mu_{misid}");
    pdf->SetParameter(4,0.042);
    //background intercept and slope
    pdf->SetParName  (5,  "N_{bkg}");
    pdf->SetParameter(5,  1700);
    pdf->SetParName  (6,  "slope");
    //just some style for drawing...
    pdf->SetLineColor(kBlue-3);
    pdf->SetLineWidth(2);
    
    //and now fit
    cout << "\n First fit, fixing all possible parameters: \n\n";
    h_data->Fit("pdf","IR");
  

    //draw the result
    gStyle->SetOptStat(0);
    gStyle->SetOptFit(1111);
    TCanvas* c1 = new TCanvas("c1","c1",600,600);
    
    h_data->SetMinimum(0);
    h_data->SetMarkerColor(kBlack);
    h_data->SetMarkerStyle(8);
    h_data->SetMarkerSize(0.8);
    h_data->SetLineColor(kBlack);
    
    h_data->Draw("err");
    
  
    return;
}