%jsroot on

// Creates a TChain to be used by the Analysis.C class
#include <TChain.h>
#include <vector>
#include <TFile.h>
#include <iostream>
#include <string>
#include <stdio.h>

TString path = "https://atlas-opendata.web.cern.ch/atlas-opendata/samples/2020/GamGam/"

TChain* fChain = new TChain("mini");

fChain->AddFile(path+"Data/data_A.GamGam.root");
fChain->AddFile(path+"Data/data_B.GamGam.root");
fChain->AddFile(path+"Data/data_C.GamGam.root");
fChain->AddFile(path+"Data/data_D.GamGam.root");

UInt_t  Photon_n = -1;  //number of preselected photons

vector<float>   *Photon_pt;  //transverse momentum of the photon
vector<float>   *Photon_eta = 0;  //pseudorapidity of the photon
vector<float>   *Photon_phi = 0;  //azimuthal angle of the photon
vector<float>   *Photon_E   = 0;  //energy of the photon
vector<bool>    *Photon_isTightID = 0;

Bool_t          TrigP = 0;
vector<float>   *Photon_ptcone30 = 0;
vector<float>   *Photon_etcone20 = 0;

fChain->SetBranchAddress("photon_pt",        &Photon_pt);  
fChain->SetBranchAddress("photon_n",         &Photon_n);  

fChain->SetBranchAddress("photon_eta",       &Photon_eta);
fChain->SetBranchAddress("photon_phi",       &Photon_phi);
fChain->SetBranchAddress("photon_E",         &Photon_E);
fChain->SetBranchAddress("photon_isTightID", &Photon_isTightID);

fChain->SetBranchAddress("trigP",            &TrigP);
fChain->SetBranchAddress("photon_ptcone30",  &Photon_ptcone30);
fChain->SetBranchAddress("photon_etcone20",  &Photon_etcone20);

//Invariant mass histograms definition
TH1F *h_M_Hyy = new TH1F("h_M_Hyy","Diphoton invariant-mass ; M_{#gamma#gamma} [GeV] ; Events / bin", 30, 105, 160);

h_M_Hyy->SetMarkerSize(2.0);
h_M_Hyy->SetLineColor(kBlue);
h_M_Hyy->SetFillColor(kBlue-10);

int nentries, nbytes, i;
nentries = (Int_t)fChain->GetEntries();

std::cout << "* Total number of entries to analyse: " << nentries << std::endl;

// Preselection of good photons

for (i=0; i < nentries; i++)
{
    nbytes =  fChain->GetEntry(i);
    
    //# Cut: Diphoton trigger is satisfied
    if(TrigP)
    {
        int goodphoton_index[5]; // give the vector a size good enough to avoid segmentation faults due to events with more than 2 photons.
        int goodphoton_n = 0;    // the total number of good photon
        int photon_index = 0;
        
        for(unsigned int j=0; j<Photon_n; j++)
        {
            // photons are tight
            if(Photon_isTightID->at(j))
            {
                //# Cut: two photons with 25 GeV and excluding the transition region between the barrel and endcap calorimeters
                if(Photon_pt->at(j) >25000. && TMath::Abs(Photon_eta->at(j))<2.37 && 
                   (TMath::Abs(Photon_eta->at(j)) < 1.37 || TMath::Abs(Photon_eta->at(j)) > 1.52))
                {
                    goodphoton_n = goodphoton_n + 1;  // count
                    goodphoton_index[photon_index] = j;
                    photon_index++;
                }
            }
        }

        //# Cut: Exactly two photons   
        if(goodphoton_n==2)
        {
            int goodphoton1_index = goodphoton_index[0];
            int goodphoton2_index = goodphoton_index[1];
            // isolated photons
            if(((Photon_ptcone30->at(goodphoton1_index)/Photon_pt->at(goodphoton1_index)) < 0.065) && ((Photon_etcone20->at(goodphoton1_index) / Photon_pt->at(goodphoton1_index)) < 0.065 ))
            {
                if(((Photon_ptcone30->at(goodphoton2_index)/Photon_pt->at(goodphoton2_index)) < 0.065) && ((Photon_etcone20->at(goodphoton2_index) / Photon_pt->at(goodphoton2_index)) < 0.065 ))
                {
                    // TLorentzVector definitions
                    TLorentzVector Photon_1  = TLorentzVector();
                    TLorentzVector Photon_2  = TLorentzVector();
                    
                    Photon_1.SetPtEtaPhiE(Photon_pt->at(goodphoton1_index), Photon_eta->at(goodphoton1_index), Photon_phi->at(goodphoton1_index),Photon_E->at(goodphoton1_index));
                    Photon_2.SetPtEtaPhiE(Photon_pt->at(goodphoton2_index), Photon_eta->at(goodphoton2_index), Photon_phi->at(goodphoton2_index),Photon_E->at(goodphoton2_index));
                    
                    float dPhi_yy = fabs(Photon_phi->at(goodphoton1_index) - Photon_phi->at(goodphoton2_index));
                    dPhi_yy       = dPhi_yy < TMath::Pi() ? dPhi_yy : 2*TMath::Pi() - dPhi_yy;
                    float m_yy    = sqrt( 2 * Photon_1.Pt() * Photon_2.Pt() * 
                                      (cosh( Photon_1.Eta() - Photon_2.Eta()) - cos(dPhi_yy)));
                    
                    //Calculation of the Invariant Mass using TLorentz vectors
                    TLorentzVector Photon_12 = Photon_1 + Photon_2;
                    float mass_inv_GeV = Photon_12.M()/1000.;
                    
                    h_M_Hyy->Fill(mass_inv_GeV);
                }
            }
        }
    } // end TrigPhoton request
}

std::cout << "* Analysed a total of: " << nentries << " in this sample." << std::endl;

TCanvas *cz = new TCanvas("cz","cz",10,10,900,600);
TText tz; tz.SetTextFont(42); tz.SetTextAlign(21);
h_M_Hyy->Draw();
cz->Draw();

h_M_Hyy->Draw("E3");
cz->SetLogy();
cz->Draw();