#define Containment_cxx
#include "Containment.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>

#include <iostream>
#include <vector>
#include <cmath>

double calculateMean(const std::vector<double>& data) {
    double sum = 0.0;
    for (const double& value : data) {
        sum += value;
    }
    return sum / data.size();
}


double calculateStandardDeviation(const std::vector<double>& data) {
    if (data.empty()) {
        // Handle the case of an empty dataset.
        return 0.0; // You can choose an appropriate value or behavior here.
    }

    double mean = calculateMean(data);
    double sumOfSquaredDifferences = 0.0;

    for (const double& value : data) {
        double difference = value - mean;
        sumOfSquaredDifferences += difference * difference;
    }

    double variance = sumOfSquaredDifferences / data.size();
    return std::sqrt(variance);
}


void Containment::Loop()
{
//   In a ROOT session, you can do:
//      root> .L Containment.C
//      root> Containment 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;

   Long64_t nentries = fChain->GetEntriesFast();

   TString root_fileName = fileName;
   TFile *fout = new TFile(root_fileName, "RECREATE");
   fout->cd();
   int n_layer = 12;
   TH1D *Edev = new TH1D("Edev", "Edev", n_layer, 0.5, n_layer + 0.5);
   
   TH1D *Edep[n_layer];
   
   Long64_t nbytes = 0, nb = 0;
   for(int ll=0; ll<n_layer; ll++){
   double Egap_l = 0, Egap_l_err = 0;
   double Eabs_l = 0, Eabs_l_err = 0;
   double n = 0;
   vector<double> Eabs_temp;
   vector<double> Egap_temp;
   Eabs_temp.clear(); Egap_temp.clear();
   
   
   
   std::stringstream s1;
   s1 << "Edep_Layer" << ll;
   Edep[ll] = new TH1D(s1.str().c_str(), s1.str().c_str(), 80, 0., 2000.);
   
   for (Long64_t jentry=0; jentry<nentries;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      
      for(int evt=0; evt <nentries; evt++){ 
      Eabs_l = 0., Egap_l = 0.; 
     // if(FirstHadInteractionPosition>-15 && FirstHadInteractionPosition < -7 && (EventID==evt)){ // Here you need to modify the values on the z-axis according to your geometry
     if((EventID==evt)){
      	        n++;
      		
      		for(int nhit=0; nhit<Egap->size(); nhit++){
      			if(Gap_Layer->at(nhit) == ll){
      				
      				Egap_l += Egap->at(nhit);
      			}
      		} 
      	     	
      		for(int nhit=0; nhit<Eabs->size(); nhit++){
      			if(Abso_Layer->at(nhit) == ll){
      		
      				Eabs_l += Eabs->at(nhit);
      			}
      		} 
      		Eabs_temp.push_back(Eabs_l);
      		Egap_temp.push_back(Egap_l);
      		Edep[ll]->Fill(Eabs_l+Egap_l);
      		
      	} 
      	        
      	}      
      } 
      
      cout << " Eabs " << Eabs_l << " Egap " << Egap_l << endl;
      cout << "n " << Egap_temp.size() << " Layer " << ll +1 << " egap " << calculateMean(Egap_temp)  << " +- " << calculateStandardDeviation(Egap_temp) << " eabs " << calculateMean(Eabs_temp) << " +- " << calculateStandardDeviation(Eabs_temp) << endl;
      Edev->SetBinContent(ll+1, calculateMean(Eabs_temp)+calculateMean(Egap_temp));
      Edev->SetBinError(ll+1, calculateStandardDeviation(Eabs_temp) );
         
   } 
   
   fout->cd();
   //Write and close the file
   fout->Write();
   fout->Close();
}