#define IonExc1_cxx
#include "IonExc1.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TGraph.h>
#include <TPad.h>

void IonExc1::Loop()
{
//   In a ROOT session, you can do:
//      root> .L IonExc1.C
//      root> IonExc1 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 = 4500;
   int nbins1 = 450;   
   gStyle->SetOptStat(1111);
   gStyle->SetPalette(57); //kRainBow
   gStyle->SetPalette(kRainBow);
   
   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);//*/
   
   // can
//   auto *ana= new TH2D("ana","ana", nbins, 0.0,0.0455, nbins, 0.0,1.); // edep
    auto *ana= new TH1D("ana","ana", nbins, 0.0,0.0455); //edep
////   // can1
//     auto *ana1= new TH1D("ana1","ana1", nbins, 0.0,0.0455); // edep
//   auto *ana1= new TH2D("ana1","ana1", nbins, 0.0,0.0455, nbins, 0.0,1.); // edep
//   // can2
//   auto *ana2 = new TGraph(); //dose - control 250000 = 10^-12
    auto *ana2= new TH2D("ana2","ana2", nbins, 0.0,0.0455, nbins, 0.0,1.); // edep
    //auto *ana2= new TH2D("ana2","ana2", nbins, 0.0,0.0455, nbins, 0.0,1.); // edep
   // can3
    auto *ana3= new TH2D("ana3","ana3", nbins, 0.0,0.0455, nbins, 0.0,1.); // edep
   //auto *ana3 = new TGraph(); //dose - control 250000 = 10^-12
   // can4
//   TH1D *ana4= new TH1D("ana4","ana4", nbins, 0.0,0.015); //edep
   //can1
   //TH2F *ana2= new TH2F("ana2","ana2", nbins, 0.0,0.0455, nbins, 0.0,23.); // edep
   //TH2F *ana3 = new TH2F("ana3","ana3", nbins,  0.0,23.0, nbins1, 0.0,10.0*pow(10,-9)); //dose - control 250000 = 10^-9
   // can2
   /*TH2F *ana4 = new TH2F("ana4","ana4",nbins1, 0.0,0.0455, nbins1, 0.0,1850.);
   TH2F *ana5 = new TH2F("ana5","ana5",nbins1, 0.0,0.0455, nbins1, 0.0,1850.);
   //can3
   TH2F *ana6= new TH2F("ana6","ana6", nbins, 0.0,0.0455, nbins, 0.0,1850.);
   TH2F *ana7 = new TH2F("ana7","ana7",nbins, 0.0,0.0455, nbins, 0.0,1850.); //  world//*/
   
      can->Divide(1,1);
//      can1->Divide(1,2);
      can2->Divide(1,3);
      can3->Divide(1,3);
//      can4->Divide(1,1);
   
//      can->cd(1);//*/
//      can1->cd(1);//*/
      can2->cd(1);//*/
      can3->cd(1);//*/
//      can4->cd(1);//*/

   Long64_t nentries = fChain->GetEntries();

   Long64_t nbytes = 0, nb = 0;
   for (Long64_t jentry=0; jentry<250;jentry++) { 
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      // if (Cut(ientry) < 0) cont1inue;
      
            cout<<"size = "<<nentries<<endl;
            cout<<" entry = "<<jentry<<endl;
                  cout<<endl;
                  
     if(Step_FinalMaterial->size()>0){mm++;
     
      Int_t cont1=0;//&&(Step_FinalMaterial->at(cont1)!="G4_WATER")&&(Step_FinalMaterial->size()==77043)
      while((cont1<Step_FinalMaterial->size()-1)){cont1++;}//*/
   
      
      //if((Step_AccumulatedEnergyDeposited->at(cont1)>0)&&(Step_FinalMaterial->at(cont1)!="G4_WATER")){mm1++;}
       //&&(Step_Particle->at(cont1)=="gamma")&&(Step_TrackID->at(cont1)==1)
       
       if((Step_FinalMaterial->at(cont1)=="G4_WATER")&&(Step_Particle->at(cont1)=="e-")){mm6++;
              
 
                for(Int_t l=0;l<Step_FinalMaterial->size();l++){mm7++;//&&(Step_Particle->at(l)=="e-") &&(Step_Particle->at(l)=="gamma") &&(Step_Particle->at(l)=="gamma") &&(Step_FinalMaterial->at(l)=="G4_WATER")
                        
                 if((Step_FinalProcess->at(l)=="e-_G4DNAIonisation")){Ion++;} 
                 if((Step_FinalProcess->at(l)=="e-_G4DNAExcitation")){Exc++;}
                 if((Step_FinalProcess->at(l)=="e-_G4DNAElastic")){Elast++;}   
                 if((Step_FinalProcess->at(l)=="e-_G4DNAVibExcitation")){VibExc++;}
                 if((Step_FinalProcess->at(l)=="e-_G4DNAAttachment")){Attac++;}
                 
                if(0.000010<(Step_InitialKineticEnergy->at(l)<0.000025)){pIon = (Ion)/(Exc+Ion+VibExc);}else{if(Step_InitialKineticEnergy->at(l)<=0.000010){pIon=(Ion+Attac)/(Exc+Ion+VibExc+Attac);}else{pIon = Ion/(Exc+Ion);}}
               if((0.000010<Step_InitialKineticEnergy->at(l)<0.000025)){pExc = (Exc+VibExc)/(Exc+Ion+VibExc);}else{if(Step_InitialKineticEnergy->at(l)<=0.000010){pExc=(Exc+VibExc+Attac)/(Exc+Ion+VibExc+Attac);}else{pExc = (Exc)/(Exc+Ion);}}
                 
                 //pDNAAttac = Attac/(Exc+Ion+Attac+VibExc);
                 //if(Step_InitialKineticEnergy->at(l)<0.01) {pVibExc = VibExc/(Exc+Ion+Attac+VibExc);
                 
//                cout<<"Energy = "<<Step_InitialKineticEnergy->at(cont1)<<" Ion = "<<Ion<<" Exc = "<<Exc<<" VIbExc = "<<VibExc<<" DNA Attac = "<<Attac<<endl;
                
                //ana1->Fill(Step_InitialKineticEnergy->at(l),pVIbExc);
                ana->Fill(Step_InitialKineticEnergy->at(l));
                ana2->Fill(Step_InitialKineticEnergy->at(l),pIon);
                ana3->Fill(Step_InitialKineticEnergy->at(l),pExc);
                 //ana3->Fill(Step_FinalKineticEnergy->at(cont1), pExc);
                 
//                        cout<<"Step size Final Kinetic = " << Step_FinalKineticEnergy->size()<<endl;
                        //cout<<"Step size Final Material = " << Step_FinalMaterial->size()<<endl;
                 }
            pIon1=pIon;
//            pIon = 0;
//            pIon2 = pIon2 + pIon1;
            pExc1=pExc;
//            pExc = 0;
//            pExc2 = pExc2 + pExc1;

//        DiffKin = Step_InitialKineticEnergy->at(cont1)-Step_FinalKineticEnergy->at(cont1);

//cout<<"Energy = "<<Step_InitialKineticEnergy->at(cont1)<<" Ion = "<<Ion<<" Exc = "<<Exc<<" VIbExc = "<<VibExc<<" DNA Attac = "<<Attac<<endl;

//                if((0.0000025<Step_InitialKineticEnergy->at(l)<0.000025)){pIon = (Ion)/(Exc+Ion+VibExc);}else{{pIon = Ion/(Exc+Ion);}}
//                if((0.0000025<Step_InitialKineticEnergy->at(l)<0.000025)){pExc = (Exc+VibExc)/(Exc+Ion+VibExc);}else{pExc = (Exc)/(Exc+Ion);}
               
//                if((0.0000025<Step_InitialKineticEnergy->at(cont1)<0.00001)) {pIon = (Ion)/(Exc+Ion+VibExc);}else{ {pIon = Ion/(Exc+Ion);}}
//                if((0.0000025<Step_InitialKineticEnergy->at(cont1)<0.00001)){pExc = (Exc+VibExc)/(Exc+Ion+VibExc);}else{pExc = (Exc)/(Exc+Ion);}

//                ana->Fill(Step_InitialKineticEnergy->at(cont1));
//                ana2->Fill(Step_InitialKineticEnergy->at(cont1),pIon1);
//                ana3->Fill(Step_InitialKineticEnergy->at(cont1),pExc1);
//            
//       ana->Fill(Step_FinalKineticEnergy->at(cont1));
//       ana1->Fill(Step_InitialKineticEnergy->at(cont1));
//       ana2->Fill(Step_FinalKineticEnergy->at(cont1),pIon1);
//       ana2->Fill(Step_FinalKineticEnergy->at(cont1),pIon);
//       sana3->Fill(Step_FinalKineticEnergy->at(cont1),pExc);
//       ana2->Fill(DiffKin);
//       ana3->Fill(Step_AccumulatedEnergyLost->at(cont1),Step_FinalKineticEnergy->at(cont1));
//       ana4->Fill(Step_AccumulatedEnergyLost->at(cont1));
       
        }//else{
//                if(Step_FinalLogicalVolume->at(cont1)=="OutOfWorld"){mm11++;}
//             }

      }
  }
  
cout<<"event = "<<mm<<endl;
cout<<"event edep!=0 ="<<mm1<<endl;
cout<<"event edep!=0 e- ="<<mm12<<endl;

cout<<"event control ="<<mm6<<endl;
cout<<"STEP Total control="<<mm7<<endl;
cout<<"eventos OutOfWorld = "<<mm11<<endl;
cout<<"####------####"<<endl;
cout<<"####---control---####"<<endl;

cout<<"Total = " << mm7 << endl;
cout<<"Ionisation = " << Ion << endl;
cout<<"Perc. Ionisation = " << pIon << endl;
cout<<"VibExcitation = " << VibExc << endl;
cout<<"Perc. VibExcitation = " << pVibExc << endl;
cout<<"Excitation = " << Exc << endl;
cout<<"Perc. Excitation = " << pExc << endl;
cout<<"DNA Attachament = " << Attac << endl;
cout<<"Perc. DNA Attachament = " << pDNAAttac << endl;
cout<<"Elastic = " << Elast << endl;
cout<<"####---####"<<endl;

//#### draw histogram

//can control
can->cd(1);
//ana2 -> SetTitle("Final Kinetic Energy spectrum x Probability Ionisation;FinalKinEnergy [MeV];Probabilty Ionisation");
ana -> SetNameTitle("energy", "Inital Kinetic Energy spectrum - Electroes Totais; InitialKinEnergy [MeV]; Counts");
//ana1 -> SetNameTitle("energy", "Kinect Energy spectrum - Electroes Primários;  Energy[MeV]; Counts");
//ana1 -> SetNameTitle("energy", "Kinect Energy spectrum - Electroes Secundários;  Energy[MeV]; Counts");
//ana1 -> SetNameTitle("energy", "Kinect Energy spectrum - Fotoes;  Energy[MeV]; 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->SetLogx();
gPad->SetLogy();
//ana2 -> ProjectionY("")->Draw();
ana -> Draw("");

//can->cd(2);
////ana2 -> SetNameTitle("energy", "Final Kinetic Energy spectrum x Probability Ionisation - Electroes Totais; Energy Lost[MeV]; FinalKinEnergy[MeV]; 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->SetLogz();
//gPad->SetLogy();
//ana -> ProjectionY("ProjY")->Draw();

/*can->cd(3);
ana2 -> SetNameTitle("energy", "Edep x Distance in to center control per event;  Energy[MeV]; Sqrt(x^2+y^2+z^2) [mm]; Counts");
ana2-> SetFillColor(0); 
gPad->SetLeftMargin(0.1);
gPad->SetRightMargin(0.15);
ana2->GetZaxis()->SetTitleOffset(1.6); 
//anaPos->GetZaxis()->SetRangeUser(-0.6, 0.6);
//gPad->SetLogz();
//gPad->SetLogy();
ana2 -> Draw("colz");

can->cd(4);
//ana1 -> SetNameTitle("energy1", "Edep x distance in to center control per event;  Energy[MeV]; Sqrt(x^2+y^2+z^2) [mm]; Counts");
ana2-> SetFillColor(0); 
gPad->SetLeftMargin(0.1);
gPad->SetRightMargin(0.15);
ana2->GetZaxis()->SetTitleOffset(1.6); 
//anaPos->GetZaxis()->SetRangeUser(-0.6, 0.6); "energy", "PosX x PosY Ti;  PosX [mm];  PosY [mm]; Edep[MeV]"
gPad->SetLogy();
ana2 -> ProjectionX("PROJECTION IN X")->Draw();
//ana17 -> Draw("colz");//*/

////////////////////////////////***********************************////////////////////////////////////
//can1->cd(1);
////ana2 -> SetTitle("Final Kinetic Energy spectrum x Probability Ionisation;FinalKinEnergy [MeV];Probabilty Ionisation");
//ana1 -> SetNameTitle("energy1a", "Final Kinetic Energy spectrum x Probability Vibration Excitation - Electroes Totais; FinalKinEnergy[MeV]; Probabilty Ionisation;  Counts");
////ana1 -> SetNameTitle("energy", "Kinect Energy spectrum - Electroes Primários;  Energy[MeV]; Counts");
////ana1 -> SetNameTitle("energy", "Kinect Energy spectrum - Electroes Secundários;  Energy[MeV]; Counts");
////ana1 -> SetNameTitle("energy", "Kinect Energy spectrum - Fotoes;  Energy[MeV]; 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();
////ana2 -> ProjectionY("")->Draw();
//ana1 -> Draw("colz");

//can1->cd(2);
////ana2 -> SetNameTitle("energy", "Final Kinetic Energy spectrum x Probability Ionisation - Electroes Totais; Energy Lost[MeV]; FinalKinEnergy[MeV]; 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 -> ProjectionY("ProjY1a")->Draw();

can2->cd(1);
//ana2 -> SetTitle("Final Kinetic Energy spectrum x Probability Ionisation;FinalKinEnergy [MeV];Probabilty Ionisation");
ana2 -> SetNameTitle("energy1", "Initial Kinetic Energy spectrum x Probability Ionisation - Electroes Totais; InitialKinEnergy [MeV]; Probabilty Ionisation;  Counts");
//ana1 -> SetNameTitle("energy", "Kinect Energy spectrum - Electroes Primários;  Energy[MeV]; Counts");
//ana1 -> SetNameTitle("energy", "Kinect Energy spectrum - Electroes Secundários;  Energy[MeV]; Counts");
//ana1 -> SetNameTitle("energy", "Kinect Energy spectrum - Fotoes;  Energy[MeV]; Counts");
ana2 -> SetFillColor(0); 
gPad -> SetLeftMargin(0.1);
gPad -> SetRightMargin(0.15);
ana2 -> GetZaxis()->SetTitleOffset(1.6); 
//anaPos->GetZaxis()->SetRangeUser(-0.6, 0.6);
gPad->SetLogx();
//gPad->SetLogy();
//ana2 -> ProjectionY("")->Draw();
ana2 -> Draw("colz");

can2->cd(2);
//ana2 -> SetNameTitle("energy", "Final Kinetic Energy spectrum x Probability Ionisation - Electroes Totais; Energy Lost[MeV]; FinalKinEnergy[MeV]; Counts");
ana2 -> SetFillColor(0); 
gPad -> SetLeftMargin(0.1);
gPad -> SetRightMargin(0.15);
ana2 -> GetZaxis()->SetTitleOffset(1.6); 
//anaPos->GetZaxis()->SetRangeUser(-0.6, 0.6);
//gPad->SetLogz();
gPad->SetLogy();
ana2 -> ProjectionY("ProjY2")->Draw();

can2->cd(3);
//ana2 -> SetNameTitle("energy", "Final Kinetic Energy spectrum x Probability Ionisation - Electroes Totais; Energy Lost[MeV]; FinalKinEnergy[MeV]; Counts");
ana2 -> SetFillColor(0); 
gPad -> SetLeftMargin(0.1);
gPad -> SetRightMargin(0.15);
ana2 -> GetZaxis()->SetTitleOffset(1.6); 
//anaPos->GetZaxis()->SetRangeUser(-0.6, 0.6);
//gPad->SetLogz();
gPad->SetLogy();
ana2 -> ProjectionX("ProjX2")->Draw();


can3->cd(1);
//ana3 -> SetTitle("Final Kinetic Energy spectrum x Probability Excitation - Electroes Totais; FinalKinEnergy [MeV]; Probabilty Excitation");
ana3 -> SetNameTitle("energy2", "Initial Kinetic Energy spectrum x Probability Excitation - Electroes Totais; InitialKinEnergyss [MeV]; Probabilty Excitation; Counts");
ana3-> SetFillColor(0); 
gPad->SetLeftMargin(0.1);
gPad->SetRightMargin(0.15);
//ana3->GetZaxis()->SetTitleOffset(1.6); 
//gPad->SetLogy();
//anaPos->GetZaxis()->SetRangeUser(-0.6, 0.6);
gPad->SetLogx();
ana3 -> Draw("colz");

can3->cd(2);
//ana3 -> SetNameTitle("energy", "Final Kinetic Energy spectrum x Probability Excitation - Electroes Totais; Energy Lost[MeV]; FinalKinEnergy[MeV]; Counts");
ana3-> SetFillColor(0); 
gPad->SetLeftMargin(0.1);
gPad->SetRightMargin(0.15);
ana3->GetZaxis()->SetTitleOffset(1.6); 
//anaPos->GetZaxis()->SetRangeUser(-0.6, 0.6);
//gPad->SetLogz();
gPad->SetLogy();
ana3 -> ProjectionY("ProjY3")->Draw();

can3->cd(3);
//ana3 -> SetNameTitle("energy", "Final Kinetic Energy spectrum x Probability Excitation - Electroes Totais; Energy Lost[MeV]; FinalKinEnergy[MeV]; Counts");
ana3-> SetFillColor(0); 
gPad->SetLeftMargin(0.1);
gPad->SetRightMargin(0.15);
ana3->GetZaxis()->SetTitleOffset(1.6); 
//anaPos->GetZaxis()->SetRangeUser(-0.6, 0.6);
//gPad->SetLogz();
gPad->SetLogy();
ana3 -> ProjectionX("ProjX3")->Draw();

//can4->cd(1);
//ana4 -> SetNameTitle("energy", "Energy Lost spectrum - Electroes Totais;  Energy[MeV]; Counts");
////ana1 -> SetNameTitle("energy", "Kinect Energy spectrum - Electroes Primários;  Energy[MeV]; Counts");
////ana1 -> SetNameTitle("energy", "Kinect Energy spectrum - Electroes Secundários;  Energy[MeV]; Counts");
////ana1 -> SetNameTitle("energy", "Kinect Energy spectrum - Fotoes;  Energy[MeV]; Counts");
//ana4-> SetFillColor(0); 
//gPad->SetLeftMargin(0.1);
//gPad->SetRightMargin(0.15);
//ana4->GetZaxis()->SetTitleOffset(1.6); 
////anaPos->GetZaxis()->SetRangeUser(-0.6, 0.6);
////gPad->SetLogz();
//gPad->SetLogy();
////ana2 -> ProjectionY("")->Draw();
//ana4 -> Draw();

// out file
TFile *p = new TFile("SPEC_TotElecDNA3A.root", "RECREATE"); // Sam 0.05mm
//TFile *p = new TFile("SPEC_1arioElectroesA.root", "RECREATE"); // Sam 0.05mm
//TFile *p = new TFile("SPEC_2arioElectroes.root", "RECREATE"); // Sam 0.05mm
//TFile *p = new TFile("SPEC_Fotoes.root", "RECREATE"); // Sam 0.05mm

ana->Write();
//ana1->Write();
//ana->SaveAs("SPEC_TotElecDNA.txt");
//ana1->SaveAs("SPEC_TotalElectroes1.txt");
//ana->SaveAs("SPEC_1arioElectroes1.txt");
//ana->SaveAs("SPEC_2arioElectroes.txt");
//ana->SaveAs("SPEC_Fotoes.txt");
//ana1->Write();
ana2->Write();
ana3->Write();
//ana4->Write();//*/

p->Write();
//p->Close();

}