#define GeV2_NEG_PION_cxx
#include "GeV2_NEG_PION.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>

void GeV2_NEG_PION::Loop()
{
//   In a ROOT session, you can do:
//      root> .L GeV2_NEG_PION.C
//      root> GeV2_NEG_PION t
//      root> t.GetEntry(12); // Fill t data members with entry number 12
//      root> t.Show();       // Show values of entry 12
//      root> t.Show(16);     // Read and show values of entry 16
//      root> t.Loop();       // Loop on all entries
//

//     This is the loop skeleton where:
//    jentry is the global entry number in the chain
//    ientry is the entry number in the current Tree
//  Note that the argument to GetEntry must be:
//    jentry for TChain::GetEntry
//    ientry for TTree::GetEntry and TBranch::GetEntry
//
//       To read only selected branches, Insert statements like:
// METHOD1:
//    fChain->SetBranchStatus("*",0);  // disable all branches
//    fChain->SetBranchStatus("branchname",1);  // activate branchname
// METHOD2: replace line
//    fChain->GetEntry(jentry);       //read all branches
//by  b_branchname->GetEntry(ientry); //read only this branch
// Assuming fout is properly initialized before use, e.g.:
// TFile* fout = new TFile("output_file.root", "RECREATE");
// Assuming fout is properly initialized before use, e.g.:
// TFile* fout = new TFile("output_file.root", "RECREATE");
// TFile* fout = new TFile("output_file.root", "RECREATE");
if (!fChain) return;

// Set branch addresses
fChain->SetBranchAddress("Egap", &Egap);
fChain->SetBranchAddress("Eabs", &Eabs);
fChain->SetBranchAddress("EventID", &EventID);
fChain->SetBranchAddress("FirstHadInteractionPosition", &FirstHadInteractionPosition);
fChain->SetBranchAddress("Gap_Layer", &Gap_Layer);
fChain->SetBranchAddress("Abso_Layer", &Abso_Layer);

int n_radius_bins = 80;       
int n_layer = 12;             
double cylinder_radius = 0.8; 
double cylinder_length = 6.0; 
int n_layer = 12;  

TH1D* Edep = new TH1D("Edep_vs_Radius", "Energy Deposition vs Transverse Radius; Radius (cm); Energy (GeV)", 
                                n_radius_bins, 0., 120.); 


Long64_t nentries = fChain->GetEntriesFast();

for (Long64_t jentry = 0; jentry < nentries; jentry++) {
    Long64_t ientry = fChain->GetEntry(jentry);
    if (ientry < 0) break;
    
    
    if (FirstHadInteractionPosition > -15 && FirstHadInteractionPosition < -7) {
        
        
        for (int layer = 0; layer < n_layer; layer++) {
            double layer_energy_gap = 0.0;  
            double layer_energy_abs = 0.0;  

           
            for (size_t nhit = 0; nhit < Egap->size(); nhit++) {
                if (Gap_Layer->at(nhit) == layer) {
                    layer_energy_gap += Egap->at(nhit);
                }
            }

           
            for (size_t nhit = 0; nhit < Eabs->size(); nhit++) {
                if (Abso_Layer->at(nhit) == layer) {
                    layer_energy_abs += Eabs->at(nhit);
                }
            }
            
          
            double total_layer_energy = layer_energy_gap + layer_energy_abs;
            
            Edep->Fill(radius, total_layer_energy);
        }
    }
}

// Create the canvas and draw the histogram
TCanvas* c1 = new TCanvas("c1", "Energy Deposition vs Transverse Radius", 800, 600);
c1->SetGrid();
Edep->Draw();

}