#include <iostream>
#include <cmath>
using namespace TMath;

// read the file you created, plot the Z boson mass peak and fit it to estimate the Z boson mass
Double_t fitFunc(Double_t* x, Double_t* par)
{
  Double_t gamma = Sqrt(Power(par[0],2)*(Power(par[0],2)+Power(par[1],2)));
  Double_t k1 = 2*Sqrt(2)*par[0]*par[1]*gamma;
  Double_t k2 = Pi()*Sqrt(Power(par[0],2)+gamma);
  Double_t k_fi = k1/k2;
  return par[2]*k_fi/(Power((x[0]*x[0]-Power(par[0],2)),2)+Power(par[0],2)*Power(par[1],2));
}

void ROOT_Intermediate_part2(){

    // read input file
    TFile *fInput = new TFile("zjet_filtered.root", "read");

    // read tree
    TTree *t = (TTree*)fInput->Get("Tdata_output");

    // declare variables to read tree
    double invar_mass = 0.;
    t->SetBranchAddress("m", &invar_mass);

    // TH1 object to see Z boson's mass spectrum
    TH1D *h = new TH1D("h", "Z boson's mass spectrum; Invariant mass (GeV); Counts", 100, 0., 150.);

    // go through tree
    int NEvents = t->GetEntries();
    for(int i = 0; i <= NEvents; i++){
        t->GetEntry(i);
        h->Fill(invar_mass);
    }

    // verify spectrum
    new TCanvas();
    h->DrawClone("hist");

    // fit gaussian
    new TCanvas();
    h->SetTitle("Gaussian fit");
    h->Fit("gaus");
    h->DrawClone();

    // fit Breit-Wigner function
    new TCanvas();
    h->SetTitle("Breit-Wigner fit");
    TF1 *func = new TF1("fitFunc",fitFunc,50,140,3);
    func->SetParameter(0,95.0);
    func->SetParameter(1,1);
    func->SetParameter(2,663050);
    h->Fit("fitFunc");
    h->DrawClone();

  

    return;
}