#include "TFile.h"
#include "TTree.h"
 
 double Edepentry=0;
 double Epentry;
 double Eg0entry;
 
 double simevents=1.e7;
 double DetEvent=0;
 double DetEventPD=0;
 double DetEventDD=0;
 double Eg0count=0;
 double Eg0acount=0;
 double Eg0bcount=0;
 double Eg0ccount=0;
 double Eg0dcount=0;
 double Eg0ecount=0;
 double Eg0fcount=0;
 double Eg0gcount=0;
 double Eg0aMin=1820;
 double Eg0aMax=1828;
 double Eg0bMin=4024;
 double Eg0bMax=4033;
  double Eg0cMin=5010;
 double Eg0cMax=5019;
  double Eg0dMin=5040;
 double Eg0dMax=5051;
 double Eg0eMin=5054;
 double Eg0eMax=5066;
 double Eg0fMin=6410;
 double Eg0fMax=6419;
 double Eg0gMin=6991;
 double Eg0gMax=7001;
 

void treefilter()
{
	ofstream results;
		results.open("Results.txt",  ios::out);
   TFile *fin = TFile::Open("Sim.root");
   if (fin == 0) {
      printf("Error: cannot open the file!\n");
      return;
   }
   fin->ls();
   TTree *tree = nullptr;
   fin->GetObject("Tree1", tree);
   if (!tree) {
      printf("Error: failed to get the tree from the file!\n");
      return;
   }
   tree->SetBranchStatus("*", 0);
   for (auto branchName : {"Edep", "Zint", "Egamma", "Egammacascade", "Pgx", "Pgy", "Pgz"}) {
      cout << "Start " << branchName << " copy" << endl;
      tree->SetBranchStatus(branchName, 1);
      cout << branchName << " Ntupla copied " << endl;
   }
    auto fout = TFile::Open("SimFiltered.root", "recreate");
   auto copytree = tree->CloneTree(0);
     const auto nentries = tree->GetEntries();
	 const auto Edepbranch = tree->GetBranch("Edep");
	 const auto Edep0branch = tree->GetBranch("Edep[0]");
	 const auto Edep1branch = tree->GetBranch("Edep[1]");
	 const auto Epbranch = tree->GetBranch("Ep");
	 const auto Eg0branch = tree->GetBranch("Eg[0]");
	 double Eg0total = Eg0branch->GetEntries();
     
       for (auto i : ROOT::TSeqI(nentries)) {
     Edepentry= Edepbranch->GetEntry(i);
     Epentry= Epbranch->GetEntry(i);
     Eg0entry=Eg0branch->GetEntry(i);
      if (Edepentry > 1000){
	         copytree->Fill();
  		}
  	  if (Edepentry > 0){
	         DetEvent++;
  		}
  		if (Edep0entry > 0){
	         DetEventPD++;
  		}
  			if (Edep1entry > 0){
	         DetEvent1D++;
  		}
  	  if (Epentry > 0){
	        Epcount++;
  		}
  		
		  if (Eg0entry > Eg0aMin && Eg0entry <Eg0aMax ){
	        Eg0acount++;
  		}
  		 else if (Eg0entry > Eg0bMin && Eg0entry <Eg0bMax ){
	        Eg0bcount++;
  		}
  		else if (Eg0entry > Eg0cMin && Eg0entry <Eg0cMax ){
	        Eg0ccount++;
  		}
  		else	if (Eg0entry > Eg0dMin && Eg0entry <Eg0dMax ){
	        Eg0dcount++;
  		}
  		else if (Eg0entry > Eg0eMin && Eg0entry <Eg0eMax ){
	        Eg0ecount++;
  		}
  		else if (Eg0entry > Eg0fMin && Eg0entry <Eg0fMax ){
	        Eg0fcount++;
  		}
  		else if (Eg0entry > Eg0gMin && Eg0entry <Eg0gMax ){
	        Eg0gcount++;
  		}
   }
      fout->Write();
      delete fout;
      
      results << "**************************************** "<< endl;
   			results << "Results "<< endl;
   	 		results << "Total simulated eventes= " << simevents <<endl;
   	 		results << "Fusion events= " << Epcount <<endl;
   	 		results << "Total Detected events  " << DetEvent<<endl;
   	 		results << "Detected events PD" << DetEventPD<<endl;
   	 		results << "Detected events DD" << DetEventDD<<endl;
   	 		results << "Total gamma= " << Eg0total <<endl;
   	 		results << "1824 keV gamma= " << Eg0acount <<endl;
   	 		results << "4028.7 keV gamma= " << Eg0bcount <<endl;
   	 		results << "5014.4 keV gamma= " << Eg0ccount <<endl;
   	 		results << "5045.8 keV gamma= " << Eg0dcount <<endl;
   	 		results << "5060.5 keV gamma= " << Eg0ecount <<endl;
   	 		results << "6414.3 keV gamma= " << Eg0fcount <<endl;
   	 		results << "6996.9 keV gamma= " << Eg0gcount <<endl;
}