#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
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();
 //   std::cout << "DEBUG: Total entries in chain = " << nentries << std::endl; ///////////////////////////// PASS
    int n_radius = 80;

    TH1D *Edep = new TH1D("Edep", "Energy Deposition; Radius (cm); Energy (GeV)", 5, 0., 800.);
 //   std::cout << "DEBUG: Created histogram with " << Edep->GetNbinsX() << " bins from " 
   //           << Edep->GetXaxis()->GetXmin() << " to " << Edep->GetXaxis()->GetXmax() << std::endl;     ///////////////////////////// PASS

    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;
     //  std::cout << "DEBUG: radius[" << i << "] = " << radius[i] << std::endl;  ///////////////////////////// PASS
    

    for (Long64_t jentry = 0; jentry < nentries; jentry++) {
        if (jentry >= 10) break;
        Long64_t ientry = LoadTree(jentry);
        if (ientry < 0) break;
        fChain->GetEntry(jentry);
    //    std::cout << "DEBUG: Processing entry " << jentry << ", EventID = " << EventID 
      //           << ", FirstHadInteractionPosition = " << FirstHadInteractionPosition << std::endl;  ///////////////////////////// PASS

        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++;
            //    std::cout << "DEBUG: Event " << evt << " passed cuts (n=" << n << ")" << std::endl;  ///////////////////////////// PASS

              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));

    // Debug: print every radius hit
   // std::cout << "DEBUG: Checking Gap hit - r=" << radius_hit << std::endl;

    if (radius_hit-1 >= radius[i] && radius_hit-1 < radius[i + 1]) {
        Egap_r[i] += Egap->at(nhit);

        // Debug: print when a hit is actually recorded
        std::cout << "DEBUG: Gap hit recorded - r=" << radius_hit 
                  << " (bin " << i << "), E=" << Egap->at(nhit) 
                  << ", total Egap_r[" << i << "]=" << Egap_r[i] << std::endl;
   }
}

                

                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));
                    if (radius_hit > radius[i] && radius_hit < radius[i + 1]) {
                        Eabs_r[i] += EAbs->at(nhit);
              //         std::cout << "DEBUG: Abs hit - r=" << 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;
     //           std::cout << "DEBUG: Bin " << i << " - Eabs_r=" << Eabs_r[i] 
                 //         << ", Egap_r=" << Egap_r[i] << ", energy=" << energy 
                   //       << ", totalEnergy=" << totalEnergy[i] << std::endl;
            }
        }
    }
  }      
    Edep->SetBinContent(i + 1, totalEnergy[i]);
//    std::cout << "DEBUG: Setting bin " << i+1 << " (r=" << radius[i] << "-" << radius[i]+10 
  //            << " cm) to " << totalEnergy[i] << std::endl;
}

gStyle->SetOptStat(1111111);

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