// 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 arg1 = 14.0/22.0; // 2 over pi
  Double_t arg2 = par[1]*par[1]*par[2]*par[2]; //Gamma=par[1]  M=par[2]
  Double_t arg3 = ((x[0]*x[0]) - (par[2]*par[2]))*((x[0]*x[0]) - (par[2]*par[2]));
  Double_t arg4 = x[0]*x[0]*x[0]*x[0]*((par[1]*par[1])/(par[2]*par[2]));
  return par[0]*arg1*arg2/(arg3 + arg4);
}

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,1.0);
    func->SetParameter(2,5.0);
    func->SetParameter(1,95.0);
    h->Fit("fitFunc");
    h->DrawClone();

  

    return;
}