#define GeantTreeAnalysis_cxx
#define n_quad 4
#define n_layer_per_quad 4
#define pos_layer_1 450
#define pos_layer_2 150
#define pos_layer_3 -150
#define pos_layer_4 -450

#include "Functions.h"
#include <cstdlib>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <TCanvas.h>
#include <TGaxis.h>
#include <TH2.h>
#include <TH3.h>
#include <TEfficiency.h>
#include <TGraph.h>
#include <TMultiGraph.h>
#include <TStyle.h>
#include <sstream>
#include "TMath.h"
#include <math.h>
#include <string>
#include <TFile.h>
#include <vector>
#include <TLegend.h>

using namespace std;


void GeantTreeAnalysis::Loop()
{
  if (fChain == 0)
    return;

	bool DELTA_RAY_CUTS=true;

  std::ostringstream tmp;
  Long64_t nentries = fChain->GetEntriesFast();
  ofstream results;
  results.open ("results_selection.txt");
  //Declaration before Loop
  TH1I* h0 = new TH1I("h0", "Cut flow", 10, 0.5, 10.5);
  TH1F* h1 = new TH1F("h1", "p_t in x-z plane", 100, 0, 1000);
  TH1F* h2 = new TH1F("h2", "Energy of Track", 100, 0, 1000);
  TH1F* h3 = new TH1F("h3", "Energydeposition along Track", 50, 0, 50);
  TH1F* h4 = new TH1F("h4", "angle between muon and electron x,z plane", 180, 0, 180);
  TH1F* Edep_layer[n_layer_per_quad * n_quad];
 // TH1F* Edep_layernocut[n_layer_per_quad * n_quad];
 // TH1F* Edep_layerallcut[n_layer_per_quad * n_quad];
  TLegend* leg = new TLegend(0.72,0.6,0.9,0.75);

  for(int i = 0; i < n_layer_per_quad * n_quad; i++){
	 char *histname = new char[30];
 	 sprintf(histname, "Edep_layer_%d",i);
 	 Edep_layer[i]=new TH1F(histname,"Energy deposited",25,0,5);
  }
  /*for(int i = 0; i < n_layer_per_quad * n_quad; i++){
	 char *histname = new char[30];
 	 sprintf(histname, "Edep_layer_%d_nocut",i);
 	 Edep_layernocut[i]=new TH1F(histname,"Energy deposited",50,0,50);
  }
  for(int i = 0; i < n_layer_per_quad * n_quad; i++){
	 char *histname = new char[30];
 	 sprintf(histname, "Edep_layer_%d_allcut",i);
 	 Edep_layerallcut[i]=new TH1F(histname,"Energy deposited",50,0,50);
  }*/
  //label bins in cut flow
h0->GetXaxis()->SetBinLabel(1, "All events");
h0->GetXaxis()->SetBinLabel(2, "Muon pos z");
h0->GetXaxis()->SetBinLabel(3, "Muon pos");
h0->GetXaxis()->SetBinLabel(4, "electron pos");
h0->GetXaxis()->SetBinLabel(5, "e Energy cut");
h0->GetXaxis()->SetBinLabel(6, "e tracklength cut");
h0->GetXaxis()->SetBinLabel(7, "electron pos 2");
h0->GetXaxis()->SetBinLabel(8, "straightness");
//Fill_TH(string HistoName, string sbin, float value_x, EventClass::HistoGraphType hgtype) --> HistoName->Fill(sbin.c_str(),value_x); Fill_TH("cutflow","All events", 1, EventClass::myTH1I);

  Long64_t jentry;
  Long64_t nbytes = 0, nb = 0;
  //Loop over all events to get a hitmap
  for ( jentry = 0; jentry < 20000; jentry++) { if(jentry%1000==0){std::cout<<" event "<<jentry<<endl;}
   Long64_t ientry = LoadTree(jentry);
   nb = fChain->GetEntry(jentry);
   nbytes += nb;
  h0->Fill("All events",1);

    if (Geant4Trajectories->at(0).GetPoints().at(Geant4Trajectories->at(0).GetPoints().size()-1).GetPosition().Z()>425){continue;}
    int passlayer = 0; h0->Fill("Muon pos z",1);
    for (int j = 0; j< Geant4Trajectories->at(0).GetPoints().size(); ++j ){
	if(Geant4Trajectories->at(0).GetPoints().at(j).GetPosition().Z()<425){passlayer= j; break;}
    }
    if (abs(Geant4Trajectories->at(0).GetPoints().at(passlayer).GetPosition().X())>1330.1||abs(Geant4Trajectories->at(0).GetPoints().at(passlayer).GetPosition().Y())>1330.1)      {continue;} //std::cout<<"Muon does not cross first active area in ev"<<jentry<<endl;continue;}
    h0->Fill("Muon pos",1);
    int numtrack = 0;
    for (int i = 0; i < Geant4Trajectories->size(); ++i){
	if(Geant4Trajectories->at(i).GetParticleName().substr(0,1) == "e"){
	  if(Geant4Trajectories->at(i).GetPoints().size()>1){
	     if (Geant4Trajectories->at(i).GetPoints().at(Geant4Trajectories->at(i).GetPoints().size()-1).GetPosition().Z()<425 & Geant4Trajectories->at(i).GetPoints().at(0).GetPosition().Z()>-425 & abs(Geant4Trajectories->at(0).GetPoints().at(0).GetPosition().X())<1330.1 & abs(Geant4Trajectories->at(0).GetPoints().at(0).GetPosition().Y())<1330.1 ){++numtrack;}
	  }
	}

    }
//std::cout<<numtrack<<endl;

if (numtrack == 0){continue;}//std::cout<<"No track found in active area in ev "<<jentry<<endl;} */
h0->Fill("electron pos",1);
    for (int k = 0; k < Geant4Trajectories->size(); ++k){
      double p_x = 0, p_z = 0, p_t = 0, depE = 0, theta = 0, twoDtot = 0, A = 0, B = 0 ;
      double v1[3], v2[3];
	if(Geant4Trajectories->at(k).GetParticleName().substr(0,1) != "e"){continue;}
	  if(Geant4Trajectories->at(k).GetPoints().size()<2){continue;}
		if (DELTA_RAY_CUTS){
			if (Geant4Trajectories->at(k).GetPoints().at(0).GetEnergy()<1) {continue;} h0->Fill("e Energy cut",1);
			if (Distance(Geant4Trajectories->at(k).GetPoints().at(0).GetPosition(),Geant4Trajectories->at(k).GetPoints().at(Geant4Trajectories->at(k).GetPoints().size()-1).GetPosition())<7){continue;} h0->Fill("e tracklength cut",1);
			if (Geant4Trajectories->at(k).GetPoints().at(Geant4Trajectories->at(k).GetPoints().size()-1).GetPosition().Z()>425 || Geant4Trajectories->at(k).GetPoints().at(0).GetPosition().Z()<-425){continue;} h0->Fill("electron pos 2",1);



			double realtracklength = 0;
				for (int m = 1; m < Geant4Trajectories->at(k).GetPoints().size(); ++m ){
					realtracklength = realtracklength + Distance(Geant4Trajectories->at(k).GetPoints().at(m-1).GetPosition(),Geant4Trajectories->at(k).GetPoints().at(m).GetPosition());
				}
				if (Distance(Geant4Trajectories->at(k).GetPoints().at(0).GetPosition(),Geant4Trajectories->at(k).GetPoints().at(Geant4Trajectories->at(k).GetPoints().size()-1).GetPosition())/realtracklength <0.98){continue;} h0->Fill("straightness",1);
		}

	    p_x = Geant4Trajectories->at(k).GetPoints().at(0).GetMomentum().X();
	    p_z = Geant4Trajectories->at(k).GetPoints().at(0).GetMomentum().Z();
	    p_t = sqrt(p_x*p_x +p_z*p_z);
	    //if (p_t!=0) {}
	  h1->Fill(p_t);

	  h2->Fill(Geant4Trajectories->at(k).GetPoints().at(0).GetEnergy());

	   //calculate angle: vetor 1 muon, vector 2 electron, starting point = vertex, end point

	v1[0]= Geant4Trajectories->at(k).GetPoints().at(0).GetPosition().X()-Geant4Trajectories->at(0).GetPoints().at(Geant4Trajectories->at(0).GetPoints().size()-1).GetPosition().X();
      //v1[1]= Geant4Trajectories->at(k).GetPoints().at(0).GetPosition().Y()-Geant4Trajectories->at(0).GetPoints().at(Geant4Trajectories->at(0).GetPoints().size()-1).GetPosition().Y();
	v1[2]= Geant4Trajectories->at(k).GetPoints().at(0).GetPosition().Z()-Geant4Trajectories->at(0).GetPoints().at(Geant4Trajectories->at(0).GetPoints().size()-1).GetPosition().Z();
	v2[0]= Geant4Trajectories->at(k).GetPoints().at(0).GetPosition().X()-Geant4Trajectories->at(k).GetPoints().at(Geant4Trajectories->at(k).GetPoints().size()-1).GetPosition().X();
      //v2[1]= Geant4Trajectories->at(k).GetPoints().at(0).GetPosition().Y()-Geant4Trajectories->at(k).GetPoints().at(Geant4Trajectories->at(k).GetPoints().size()-1).GetPosition().Y();
	v2[2]= Geant4Trajectories->at(k).GetPoints().at(0).GetPosition().Z()-Geant4Trajectories->at(k).GetPoints().at(Geant4Trajectories->at(k).GetPoints().size()-1).GetPosition().Z();

	twoDtot = ((v1[0]*v2[0])+(v1[2]*v2[2]));
	A = pow(v1[0],2)+pow(v1[2],2);
	A = sqrt(A);
	B = pow(v2[0],2)+pow(v2[2],2);
	B = sqrt(B);
	theta = acos(twoDtot/(A*B))*180.0 / M_PI;//min( acos(twoDtot/(A*B))*180.0 / M_PI, (180 - acos(twoDtot/(A*B))*180.0/M_PI));
	//std::cout<<theta<<" Degree"<<endl;
	h4->Fill(theta);

	  for (int l = 0; l < Geant4Hits->size(); ++l){
	   if (Geant4Hits->at(l).GetTrackID() == Geant4Trajectories->at(k).GetTrackID()){depE = depE + Geant4Hits->at(l).GetEnergyDep();}
	  }
	h3->Fill(depE);

	for (int layer = 1; layer < n_layer_per_quad * n_quad +1; ++layer){
	  for (int l = 0; l < Geant4Hits->size(); ++l){
	   if (Geant4Hits->at(l).GetTrackID() == Geant4Trajectories->at(k).GetTrackID()& WhichLayer(Geant4Hits->at(l))== layer)
			{depE = depE + Geant4Hits->at(l).GetEnergyDep();
		//if(depE>5& layer !=1) {std::cout<<"energy deposited"<<depE<<" in layer "<<layer<<" trackid"<<Geant4Hits->at(l).GetTrackID()<<"in ev "<<jentry<<" particle type "<<Geant4Hits->at(l).GetParticleName()<<"creation process "<<Geant4Hits->at(l).GetProcessName()<<endl;}
		string layer_name = Geant4Hits->at(l).GetSEName();
		int layer_number = 0;
		 if (layer_name.substr(0, 9) == "quadlayer"){
			layer_number = atoi(layer_name.substr(9, 2).c_str());
		 }
		//std::cout<<WhichQuad(Geant4Hits->at(l))<<" in quadlayer "<<WhichLayer(Geant4Hits->at(l))<< " in layer loop "<<layer<<" in ev "<<jentry<<" hit in "<<layer_name<<"Track ID "<<Geant4Hits->at(l).GetTrackID()<<" Hit z position "<<Geant4Hits->at(l).GetPoint().GetPosition().Z()<<endl;
	   }
	  }
	if (depE != 0){Edep_layer[layer-1]->Fill(depE);}
	  depE = 0;
	}


	    //h2->Fill(Geant4Trajectories->at(k).GetPoints().at(0).GetEnergy(),Distance(Geant4Trajectories->at(k).GetPoints().at(0).GetPosition(),Geant4Trajectories->at(k).GetPoints().at(Geant4Trajectories->at(k).GetPoints().size()-1).GetPosition()));
	   // h3->Fill(Geant4Trajectories->at(k).GetPoints().size(),Distance(Geant4Trajectories->at(k).GetPoints().at(0).GetPosition(),Geant4Trajectories->at(k).GetPoints().at(Geant4Trajectories->at(k).GetPoints().size()-1).GetPosition()));


	  //}
	  //std::cout<<p_t<<endl;

	//}
	//}}} //ifcheck
    }
  }
    ///end of brakets for eventdisplay loop
gStyle->SetOptStat();
leg->AddEntry(Edep_layer[1],"sTGC1","f");
leg->AddEntry(Edep_layer[5],"sTGC2","f");
leg->AddEntry(Edep_layer[9],"sTGC3","f");
leg->AddEntry(Edep_layer[13],"sTGC4","f");

  TCanvas* c1 = new TCanvas("c1", "p_t no cut", 800, 800);
  h1->Draw("E");
  h1->GetXaxis()->SetTitle("p_t (MeV)");
  h1->GetYaxis()->SetTitle("#tracks");
  h1->GetYaxis()->SetTitleOffset( 1.5 );
  h1->SetStats(1);
  c1->SetLogy();
  c1->SetLogx();
  if (DELTA_RAY_CUTS) c1->Print("p_t_cutlog.pdf");
	else c1->Print("p_t_nocutnolog.pdf");
  TCanvas* c2 = new TCanvas("c2", "Energy of Track", 800, 800);
  h2->Draw("E");
  h2->GetXaxis()->SetTitle("Energy (MeV)");
  h2->GetYaxis()->SetTitle("#tracks");
  h2->GetYaxis()->SetTitleOffset( 1.5 );
  h2->SetStats(1);
 c2->SetLeftMargin(1000);
c2->SetRightMargin(10);
  c2->SetLogz();
   c2->SetLogx();
   c2->SetLogy();
   if (DELTA_RAY_CUTS) c2->Print("Trackenergy_cutlog.pdf");
	 else c2->Print("Trackenergy_nocutlog.pdf");
  //
TCanvas* c3 = new TCanvas("c3", "Deposited_energy_nocut", 800, 800);
  h3->Draw("E");
  h3->GetXaxis()->SetTitle("Deposited_energy (MeV)");
  h3->GetYaxis()->SetTitle("#tracks");
  h3->GetYaxis()->SetTitleOffset( 1.5 );
  c3->SetLogy();
  h3->SetStats(1);
  if (DELTA_RAY_CUTS) c3->Print("depE_cutlog.pdf");
  else c3->Print("depE_nocutlog.pdf");
TCanvas* c4 = new TCanvas("c4", "Opening Angle 2D", 800, 800);
  h4->Draw("E");
  h4->GetXaxis()->SetTitle("Opening angle (deg)");
  h4->GetYaxis()->SetTitle("#tracks");
  h4->GetYaxis()->SetTitleOffset( 1.5 );
  h4->SetStats(1);
  //c4->SetLogy();
	if (DELTA_RAY_CUTS) c4->Print("theta_cut.pdf");
	else c4->Print("theta_nocut.pdf");
  results.close();
TCanvas* c0 = new TCanvas("c0", "Cutflow", 800, 800);
  h0->Draw("HISTTEXT");
  h0->GetXaxis()->SetTitle("Cuts");
  h0->GetYaxis()->SetTitle("#events");
  h0->GetYaxis()->SetTitleOffset( 1.5 );
  //c3->SetLogy();
  h0->SetStats(1);
  c0->Print("cutflowsel.pdf");
TCanvas* csTGC1 = new TCanvas("csTGC1", "Energydeposited M3SP", 800, 800);
  csTGC1->Divide(2,2);
  for (int i = 0 ; i < 4 ; ++i){
 csTGC1->cd(i+1);
 Edep_layer[i]->Draw("E");
 Edep_layer[i]->SetLineColor(1);
 Edep_layer[i+4]->Draw("SAME,E");
 Edep_layer[i+4]->SetLineColor(kRed);
 Edep_layer[i+8]->Draw("SAME,E");
 Edep_layer[i+12]->Draw("SAME,E");
 Edep_layer[i+12]->SetLineColor(kGreen+2);
 Edep_layer[i]->GetXaxis()->SetTitle("Deposited_energy (MeV)");
 Edep_layer[i]->GetYaxis()->SetTitle("#tracks");
 Edep_layer[i]->GetYaxis()->SetTitleOffset( 1.5 );
 Edep_layer[i]->SetStats(1);
 leg->Draw();}

 //csTGC1->SetLogy();
 csTGC1->Print("depE_sTGCsel.pdf");

TCanvas* csTGC2 = new TCanvas("csTGC2", "Energydeposited M3SC", 800, 800);
  csTGC2->Divide(2,2);
  for (int i = 4 ; i < 8 ; ++i){
 csTGC2->cd(i-3);
 Edep_layer[i]->Draw("E");
 Edep_layer[i]->GetXaxis()->SetTitle("Deposited_energy (MeV)");
 Edep_layer[i]->GetYaxis()->SetTitle("#tracks");
 Edep_layer[i]->GetYaxis()->SetTitleOffset( 1.5 );
 Edep_layer[i]->SetStats(1);}
 csTGC2->SetLogy();
 csTGC2->Print("depE_sTGC2.pdf");

TCanvas* csTGC3 = new TCanvas("csTGC3", "Energydeposited M2LP", 800, 800);
  csTGC3->Divide(2,2);
  for (int i = 8 ; i < 12 ; ++i){
 csTGC3->cd(i-7);
 Edep_layer[i]->Draw("E");
 Edep_layer[i]->GetXaxis()->SetTitle("Deposited_energy (MeV)");
 Edep_layer[i]->GetYaxis()->SetTitle("#tracks");
 Edep_layer[i]->GetYaxis()->SetTitleOffset( 1.5 );
 Edep_layer[i]->SetStats(1);}
 csTGC3->SetLogy();
 csTGC3->Print("depE_sTGC3.pdf");

TCanvas* csTGC4 = new TCanvas("csTGC4", "Energydeposited M2LC", 800, 800);
  csTGC4->Divide(2,2);
  for (int i = 12 ; i < 16 ; ++i){
 csTGC4->cd(i-11);
 Edep_layer[i]->Draw("E");
 Edep_layer[i]->GetXaxis()->SetTitle("Deposited_energy (MeV)");
 Edep_layer[i]->GetYaxis()->SetTitle("#tracks");
 Edep_layer[i]->GetYaxis()->SetTitleOffset( 1.5 );
 Edep_layer[i]->SetStats(1);}
 csTGC4->SetLogy();
 csTGC4->Print("depE_sTGC4.pdf");
}