#define EnerKin_cxx
#include "EnerKin.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>

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

   
   int nbins = 2000;
   gStyle->SetOptStat(1111);
   gStyle->SetPalette(57);
   
   TCanvas *can = new TCanvas("can","can",900,700); // control
   TCanvas *can1 = new TCanvas("can1","can1",900,700); //Ag
   /*TCanvas *can2 = new TCanvas("can2","can2",900,700); //I
   TCanvas *can3 = new TCanvas("can3","can3",900,700); //air
   TCanvas *can4 = new TCanvas("can4","can4",900,700); //Ti
   TCanvas *can5 = new TCanvas("can5","can5",900,700); //world */
     can->SetFrameBorderMode(1);
     can->SetLeftMargin(0.01);

     can1->SetFrameBorderMode(1);
     can1->SetLeftMargin(0.01);//*/
     
     /*can2->SetFrameBorderMode(1);
     can2->SetLeftMargin(0.01);
     
     can3->SetFrameBorderMode(1);
     can3->SetLeftMargin(0.01);
     
     can4->SetFrameBorderMode(1);
     can4->SetLeftMargin(0.01);

     can5->SetFrameBorderMode(1);
     can5->SetLeftMargin(0.01);//*/
   
   //bins more than 800 doesn't work to TH3F
   // can
   TH2F *ana= new TH2F("ana","ana", nbins, 0.0,0.0455, nbins, 0.0,23.0);
   TH2F *ana1 = new TH2F("ana1","ana1",nbins, 0.0,0.0455, nbins, 0.0,23.0);
   TH1F *an = new TH1F("an","an",nbins, 0.0,0.0455);
   //can1
   TH2F *ana2 = new TH2F("ana2","ana2",nbins, 0.0,0.0455, nbins, 0.0,1850.0);
   TH2F *ana3 = new TH2F("ana3","ana3",nbins, 0.0,0.0455, nbins, 0.0,1850.0);
   
   // can1
   /*TH2F *ana4= new TH2F("ana4","ana4", nbins, 0.0,0.0455, nbins, 0.0,0.5);
   TH2F *ana5 = new TH2F("ana5","ana5",nbins, 0.0,0.0455, nbins, 0.0,0.5);
   TProfile2D *ana6 = new TProfile2D("anaPos","anaPos",nbins, -1.0,1.0, nbins, -1.0,1.0);
   TProfile2D *ana7 = new TProfile2D("anaPos","anaPos",nbins, -1.0,1.0, nbins, -1.0,1.0);//*/
   
   // can2
   /*TH2F *ana8= new TH2F("ana8","ana8", nbins, 0.0,0.0455, nbins, 0.0,0.5);
   TH2F *ana9 = new TH2F("ana9","ana9",nbins, 0.0,0.0455, nbins, 0.0,0.5);
   TProfile2D *ana10 = new TProfile2D("anaPos","anaPos",nbins, -1.0,1.0, nbins, -1.0,1.0);
   TProfile2D *ana11 = new TProfile2D("anaPos","anaPos",nbins, -1.0,1.0, nbins, -1.0,1.0);//*/
   
   // can3
   /*TH2F *ana12= new TH2F("ana12","ana12", nbins, 0.0,0.0455, nbins, 0.0,0.5);
   TH2F *ana13 = new TH2F("ana13","ana13",nbins, 0.0,0.0455, nbins, 0.0,0.5);
   TProfile2D *ana14 = new TProfile2D("anaPos","anaPos",nbins, -1.0,1.0, nbins, -1.0,1.0);
   TProfile2D *ana15 = new TProfile2D("anaPos","anaPos",nbins, -1.0,1.0, nbins, -1.0,1.0);//*/
   
   // can4
   /*TH2F *ana16= new TH2F("ana16","ana16", nbins, 0.0,0.0455, nbins, 0.0,0.45);
   TH2F *ana17 = new TH2F("ana17","ana17",nbins, 0.0,0.0455, nbins, 0.0,0.45);
   TProfile2D *ana18 = new TProfile2D("anaPos","anaPos",nbins, -1.0,1.0, nbins, -1.0,1.0);
   TProfile2D *ana19 = new TProfile2D("anaPos","anaPos",nbins, -1.0,1.0, nbins, -1.0,1.0);//*/
   
   // can5
   /*TH2F *ana20= new TH2F("ana20","ana20", nbins, 0.0,0.0455, nbins, 0.0,1850.0);
   TH2F *ana21 = new TH2F("ana21","ana21",nbins, 0.0,0.0455, nbins, 0.0,1850.0);
   TProfile2D *ana22 = new TProfile2D("anaPos","anaPos",nbins, -1.0,1.0, nbins, -1.0,1.0);
   TProfile2D *ana23 = new TProfile2D("anaPos","anaPos",nbins, -1.0,1.0, nbins, -1.0,1.0);//*/

can->Divide(1,2);
can1->Divide(1,2);
/*can2->Divide(2,2);
can3->Divide(2,2);
can4->Divide(2,2);
can5->Divide(2,2);//*/
   
can->cd(1);//*/
can1->cd(1);//*/
//can2->cd(1);//*/
//can3->cd(1);//*/
//can4->cd(1);//*/
//can5->cd(1);//*/

   Long64_t nentries = fChain->GetEntries(); //para 25000 event => 202287 steps, deu 209435 1997243

   for (Long64_t jentry=0; jentry<nentries;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;  
      nb = fChain->GetEntry(jentry);   nbytes += nb;
       // if (Cut(ientry) < 0) continue;
      
      cout<<"size total = " <<nentries<<endl;
      cout<<"entry = "<<jentry<<endl; 
      cout<<endl;
      
      if(Step_FinalMaterial->size()>0){ mm++; // garantir que o size não é zero 
      
      Int_t cont=0; Int_t count=0;
      while((cont<Step_FinalMaterial->size()-1)||(cont<Step_AccumulatedEnergyDeposited->size()-1)){cont++;}
      
      //ana1->Fill(Step_FinalKineticEnergy->at(cont),TMath::Sqrt(pow(Step_FinalPosX->at(cont),2)+pow(Step_FinalPosY->at(cont),2)+pow(Step_FinalPosZ->at(cont),2))); ||((TMath::Sqrt(pow(Step_FinalPosX->at(cont),2)+pow(Step_FinalPosY->at(cont),2)+pow(Step_FinalPosZ->at(cont),2))>0.40)&&((TMath::Sqrt(pow(Step_FinalPosX->at(cont),2)+pow(Step_FinalPosY->at(cont),2)+pow(Step_FinalPosZ->at(cont),2)))<22.))
      
      if((Step_FinalMaterial->at(cont)=="G4_WATER")){mm2++;
            
            ana->Fill(Step_FinalKineticEnergy->at(cont),TMath::Sqrt(pow(Step_FinalPosX->at(cont),2)+pow(Step_FinalPosY->at(cont),2)+pow(Step_FinalPosZ->at(cont),2)));
           for(Int_t l=0;l<Step_FinalMaterial->size();l++){
            ekin=ekin+Step_FinalKineticEnergy->at(l);
            ana->Fill(Step_FinalKineticEnergy->at(l),TMath::Sqrt(pow(Step_FinalPosX->at(l),2)+pow(Step_FinalPosY->at(l),2)+pow(Step_FinalPosZ->at(l),2)));
            }
            sumKin=sumKin+Step_FinalKineticEnergy->at(cont);
            ana1->Fill(Step_FinalKineticEnergy->at(cont),TMath::Sqrt(pow(Step_FinalPosX->at(cont),2)+pow(Step_FinalPosY->at(cont),2)+pow(Step_FinalPosZ->at(cont),2)));
           
            
         }
      }
   }
    cout<<endl;
   
   cout<< "number first if = "<<mm<<endl;
   cout<< "number of edep dif of zero = "<<mm1<<endl;
   cout<< "event control = "<<mm2<<endl;
      cout<< "step control = "<<mm8<<endl;
         cout<< "step control edep != 0 = "<<mm9<<endl;
   
   /*cout<< "event world = "<<mm7<<endl;
         cout<< "step world = "<<mm18<<endl;
         cout<< "step world edep != 0 = "<<mm19<<endl;
   cout<< "event OOW = "<<mm20<<endl;
         cout<< "step OOW = "<<mm21<<endl;
         cout<< "step OOW edep != 0 = "<<mm22<<endl;//*/
         
//cout<< "creator process in control phot= "<<mm23<<endl;    

//can
can->cd(1);
ana -> SetNameTitle("energy", "Ekin x distance in to center Control STEP;  Energy[MeV];  Sqrt(x^2+y^2+z^2) [mm]; Counts");
ana -> SetFillColor(0); 
gPad->SetLeftMargin(0.1);
gPad->SetRightMargin(0.15);
ana->GetZaxis()->SetTitleOffset(1.6); 
//anaPos->GetZaxis()->SetRangeUser(-0.6, 0.6); 
//gPad->SetLogy();
ana -> Draw("colz");

can->cd(2);
//ana -> SetNameTitle("energy", "Edep x distance in to center Control projection;  Energy[MeV]; Sqrt(x^2+y^2+z^2) [mm]; Counts");
ana-> SetFillColor(0); 
gPad->SetLeftMargin(0.1);
gPad->SetRightMargin(0.15);
ana->GetZaxis()->SetTitleOffset(1.6); 
//anaPos->GetZaxis()->SetRangeUser(-0.6, 0.6); "energy", "PosX x PosY Ti;  PosX [mm];  PosY [mm]; Edep[MeV]"
gPad->SetLogy();
ana -> ProjectionX("pjcX")->Draw();
//ana17 -> Draw("colz");//*/

can1->cd(1);
ana1-> SetNameTitle("energy", "Ekin x distance in to center Control Event;  Energy[MeV]; Sqrt(x^2+y^2+z^2) [mm]; Counts");
ana1 -> SetFillColor(0); 
gPad->SetLeftMargin(0.1);
gPad->SetRightMargin(0.15);
ana1->GetZaxis()->SetTitleOffset(1.6); 
//anaPos->GetZaxis()->SetRangeUser(-0.6, 0.6); 
//gPad->SetLogy();
ana1 -> Draw("colz");

can1->cd(2);
//ana4 -> SetNameTitle("energy", "Edep x distance in to center Ag;  Energy[MeV]; Sqrt(x^2+y^2+z^2) [mm]; Counts");
ana1-> SetFillColor(0); 
gPad->SetLeftMargin(0.1);
gPad->SetRightMargin(0.15);
ana1->GetZaxis()->SetTitleOffset(1.6); 
//anaPos->GetZaxis()->SetRangeUser(-0.6, 0.6);
//gPad->SetLogz();
gPad->SetLogy();
ana1 -> ProjectionX("")->Draw();
//ana4 -> Draw("colz");


// out file
TFile *p = new TFile("EnerKin.root", "RECREATE"); // Sam 0.05mm

ana->Write();

}