//-------------------------------------------------------------------------------------
// skeleton_cms_dimuon.cc
//
// This is a skeleton program to read/analyse CMS files via root or standalnoe c++
//
// * You can run this macro using
//     $ root skeleton_cms_dimuon.cc
//
// * or if you can compile this file as follows (you may require root-shell):
//     $ root-shell
//     $ g++ skeleton_cms_dimuon.cc -o skeleton_cms_dimuon `root-config --cflags --glibs`
//   and run:
//     $ ./skeleton_cms_dimuon
//
// Apr 2026
//-------------------------------------------------------------------------------------

#include "ROOT/RDataFrame.hxx"
#include "ROOT/RDFHelpers.hxx"
#include "ROOT/RVec.hxx"
#include "TCanvas.h"
#include "TH1D.h"
#include "TLatex.h"
#include "TStyle.h"
#include "TF1.h"
#include <iostream>
using namespace ROOT::VecOps;
void skeleton_cms_dimuon()
{
   // Enable multi-threading
   ROOT::EnableImplicitMT();

   // Create dataframe from NanoAOD files
   ROOT::RDataFrame df("Events", "/hepdata/CMS/DoubleMuon/Run2012BC_DoubleMuParked_Muons.root");
   // Get the number of events
   int event_count = static_cast<int>(*df.Count());

    // Output the number of events
   std::cout << "Number of events: " << event_count << std::endl;

   //auto df = df_original.Range(event_count);


   // For simplicity, select only events with exactly two muons and require opposite charge
   auto df_2mu = df.Filter("nMuon == 2", "Events with exactly two muons");
   auto df_os = df_2mu.Filter("Muon_charge[0] != Muon_charge[1]", "Muons with opposite charge");

   // Compute invariant mass of the dimuon system
   auto df_mass = df_os.Define("Dimuon_mass", InvariantMass<float>, {"Muon_pt", "Muon_eta", "Muon_phi", "Muon_mass"});

   // Make histogram of dimuon mass spectrum. Note how we can set title and axis labels in one go
   auto h = df_mass.Histo1D({"Dimuon_mass", "Dimuon mass;m_{#mu#mu} (GeV);N_{Events}", 30000, 0.25, 300}, "Dimuon_mass");

   // Request cut-flow report
   auto report = df.Report();

   // Produce plot
   gStyle->SetOptStat(0); gStyle->SetTextFont(42);
   auto c = new TCanvas("c", "", 800, 700);
   c->SetLogx(); c->SetLogy();

   h->GetXaxis()->SetTitleSize(0.04);
   h->GetYaxis()->SetTitleSize(0.04);
   h->DrawClone();

   // Define Breit-Wigner function: A / [(x^2 - M^2)^2 + (M * Gamma)^2]
   TF1 *bwFit = new TF1("bwFit", "[0] / ((x*x - [1]*[1])*(x*x - [1]*[1]) + ([1]*[2])*([1]*[2]))", 2.95, 3.25);
   bwFit->SetParameters(1e5, 3.1, 0.1);  // Initial guesses: norm, mass (j/psi), width

   // Perform the fit
   h->Fit(bwFit, "R");
   bwFit->Draw("same");

   TLatex label; label.SetNDC(true);
   label.DrawLatex(0.175, 0.740, "#eta");
   label.DrawLatex(0.205, 0.775, "#rho,#omega");
   label.DrawLatex(0.270, 0.740, "#phi");
   label.DrawLatex(0.400, 0.800, "J/#psi");
   label.DrawLatex(0.415, 0.670, "#psi'");
   label.DrawLatex(0.485, 0.700, "Y(1,2,3S)");
   label.DrawLatex(0.755, 0.680, "Z");
   label.SetTextSize(0.040); label.DrawLatex(0.100, 0.920, "#bf{CMS Open Data}");
   label.SetTextSize(0.030); label.DrawLatex(0.630, 0.920, "#sqrt{s} = 8 TeV, L_{int} = 11.6 fb^{-1}");

   c->SaveAs("dimuon_spectrum.pdf");
   c->SaveAs("dimuon_spectrum.png");
   c->SaveAs("dimuon_spectrum.root");

   // Print cut-flow report
   report->Print();

   // Print fit results
   std::cout << "Fit Results:    " << std::endl;
   std::cout << "Normalization:  " << bwFit->GetParameter(0) << std::endl;
   std::cout << "Jpsi Mass:      " << bwFit->GetParameter(1) << " GeV" << std::endl;
   std::cout << "Jpsi Width:     " << bwFit->GetParameter(2) << " GeV" << std::endl;
   std::cout << "Jpsi Life-Time: " << 6.582*1e-22/(1000*bwFit->GetParameter(2)) << " GeV" << std::endl;
}

int main()
{
  skeleton_cms_dimuon();
}