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

#include <TFile.h>
#include <TCanvas.h>
#include <TChain.h>
#include <TH1F.h>
#include <TCanvas.h>
#include <TLorentzVector.h>

#include <vector>
#include <iostream>
using namespace std;

#include "analysis.C"
//#include "mini.C"

#define GeV 1e-3

void skeleton_atlas()
{
    // Enable multi-threading
    ROOT::EnableImplicitMT();

    // Open root file to save histogram
    TFile *fout = new TFile("mgg.root","recreate");

    // Create chain
    //TChain *chain = new TChain("mini");
    //chain->Add("GamGam/Data/*.root");
    TChain *chain = new TChain("analysis");
    chain->Add("ODEO_*.root");

    // Photon branches
    std::vector<float> *ph_pt  = nullptr;
    std::vector<float> *ph_eta = nullptr;
    std::vector<float> *ph_phi = nullptr;
    std::vector<float> *ph_E   = nullptr;

    chain->SetBranchAddress("photon_pt",  &ph_pt);
    chain->SetBranchAddress("photon_eta", &ph_eta);
    chain->SetBranchAddress("photon_phi", &ph_phi);
    chain->SetBranchAddress("photon_e",   &ph_E);

    // Histogram
    TH1F *h_mgg = new TH1F(
        "h_mgg",
        "Diphoton invariant mass; m_{#gamma#gamma} [GeV]; Events",
        240, 100, 160
    );

    analysis *ana = new analysis(chain);
    //mini *ana = new mini(chain);


    // Event loop
    Long64_t nentries = chain->GetEntries();
    cout  << "Nevents: " << nentries << endl;

    for (Long64_t ev = 0; ev < nentries; ev++) {
        ana->GetEntry(ev);

        if(ev%50000==0) cout << ev << " / " << nentries <<" [" << ceil(100.0*ev / nentries) << "%]" << "\r"<<flush;

        //if (ana->photon_pt->size() != 2) continue;

        if( abs(ana->photon_eta->at(0))>1.37 && abs(ana->photon_eta->at(0))<1.52) continue;
        if( abs(ana->photon_eta->at(1))>1.37 && abs(ana->photon_eta->at(1))<1.52) continue;
        if( ana->photon_isTightID->at(0)==false) continue;
        if( ana->photon_isTightID->at(1)==false) continue;
        if( ana->photon_ptcone20->at(0)/ana->photon_pt->at(0) > 0.055) continue;
        if( ana->photon_ptcone20->at(1)/ana->photon_pt->at(1) > 0.055) continue;
        if( ana->photon_e->at(0) < 50) continue;
        if( ana->photon_e->at(1) < 30) continue;

        TLorentzVector g1, g2;

        g1.SetPtEtaPhiE(
            ana->photon_pt->at(0),
            ana->photon_eta->at(0),
            ana->photon_phi->at(0),
            ana->photon_e->at(0)
        );

        g2.SetPtEtaPhiE(
            ana->photon_pt->at(1),
            ana->photon_eta->at(1),
            ana->photon_phi->at(1),
            ana->photon_e->at(1)
        );

        TLorentzVector gg = g1 + g2;
        double m = gg.M();

        if(ana->photon_pt->at(0)/m < 0.35) continue;
        if(ana->photon_pt->at(1)/m < 0.35) continue;

        h_mgg->Fill( m );
    }

    fout->Write("h_mgg");

    TCanvas *c = new TCanvas("c", "Diphoton mass", 800, 600);
    h_mgg->Draw("error");
    c->SaveAs("mgg_chain.png");
    std::cout << "\nMerged data plotted" << std::endl;
}

int main(){
  skeleton_atlas();
  return 0;
}