#define GeV5_NEG_PION_cxx
#include "GeV5_NEG_PION.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>

void GeV5_NEG_PION::Loop()
{
//   In a ROOT session, you can do:
//      root> .L GeV5_NEG_PION.C
//      root> GeV5_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
 if (fChain == 0) return;
fChain->SetBranchAddress("Egap", &Egap);
fChain->SetBranchAddress("EAbs", &EAbs);
fChain->SetBranchAddress("EventID", &EventID);
fChain->SetBranchAddress("FirstHadInteractionPosition", &FirstHadInteractionPosition);
fChain->SetBranchAddress("Gap_X", &Gap_X);
fChain->SetBranchAddress("Gap_Y", &Gap_Y);
fChain->SetBranchAddress("Abs_X", &Abs_X);
fChain->SetBranchAddress("Abs_Y", &Abs_Y);
Long64_t nentries = fChain->GetEntriesFast();

int n_radius = 80;
TH1D *Edep = new TH1D("Edep", "Energy Deposition; Radius (cm); Energy (GeV)", n_radius, 0., 800.);

std::vector<double> radius(n_radius, 0.0); 
std::vector<double> totalEnergy(n_radius, 0.0);

int n = 0;


for (int i = 0; i < n_radius; i++) { 
 
    radius[i] = i * 10;


for (Long64_t jentry = 0; jentry < nentries; jentry++) {
    if (jentry > 100) break;
    Long64_t ientry = LoadTree(jentry);
    if (ientry < 0) break;
    fChain->GetEntry(jentry);
    
    for (int evt = 0; evt < nentries; evt++) { 
    
    
        std::vector<double> Eabs_r(n_radius, 0.0);
        std::vector<double> Egap_r(n_radius, 0.0);
     
    
   
        if (FirstHadInteractionPosition > -15 && FirstHadInteractionPosition < -7 && EventID == evt) {
            n++;
            
            for (size_t nhit = 0; nhit < Egap->size(); nhit++) {
                double radius_hit = sqrt(pow(Gap_X->at(nhit), 2) + pow(Gap_Y->at(nhit), 2));
                for (int i = 0; i < n_radius; i++) {
                    if (radius_hit < radius[i-1] && radius_hit > radius[i]) {
                        Egap_r[i] += Egap->at(nhit);
                    }
                }
            }
            
            for (size_t nhit = 0; nhit < EAbs->size(); nhit++) {
                double radius_hit = sqrt(pow(Abs_X->at(nhit), 2) + pow(Abs_Y->at(nhit), 2));
                for (int i = 0; i < n_radius; i++) {
                    if (radius_hit < radius[i-1] && radius_hit > radius[i]) {
                        Eabs_r[i] += EAbs->at(nhit);
                        std::cout << "DEBUG: Abs hit=" << radius_hit << " (bin " << i << "), E=" 
                                  << EAbs->at(nhit) << ", total Eabs_r[" << i << "]=" << Eabs_r[i] << std::endl;
                    }
                }
            }
        
            for (int i = 0; i < n_radius; i++) {
                double energy = (Eabs_r[i] + Egap_r[i]) / 2.0 * 1000.0;
                totalEnergy[i] += energy;
            }
       }
          for (int i = 0; i < n_radius; i++) {
    Edep->SetBinContent(i , totalEnergy[i]);
}
    }
 
}


}

TCanvas *c1 = new TCanvas("c1", "Energy Deposition", 800, 600);
Edep->Draw("HIST");

}