#define elec_geo20_cxx
#include "elec_geo20.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMath.h>
#include <iostream>
#include <vector>

using namespace std;
using std::vector;


void elec_geo20::Loop(string num)
{

using namespace std;
using std::vector;
	
 if (!(gInterpreter->IsLoaded("vector")))
    gInterpreter->ProcessLine("#include <vector>");
  gSystem->Exec("rm -f AutoDict*vector*vector*float*");
  gInterpreter->GenerateDictionary("vector<vector<float> >", "vector");	
	
// Flags
bool writeOutput = true;
bool doReconstruction = false;	
bool doTriggerfilling = false;	

//TFile *newfile = new TFile("outpumt_electrons300mev_histo.root", "recreate");
if(writeOutput) {
TFile *newfile = new TFile(("output_electrons"+num+".recov08_histo.root").c_str(), "recreate");	
cout << "creating output file "+num+" mev" << endl;
}

const int ptbin = 2000;
const int etabin = 2000;
const int mubin = 43;

const int ptxmin = 0;
const int ptxmax = 300000;
const int etaxmin = -4;
const int etaxmax = 4;
const int muxmin = -1.5;
const int muxmax = 41.5;

TH1F *pt_mc = new TH1F("pt_mc", "pt_mc", 2000, 0, ptxmax);
TH1F *eta_mc = new TH1F("eta_mc", "eta_mc", 2000, etaxmin, etaxmax);
TH1F *phi_mc = new TH1F("phi_mc", "phi_mc", 2000, etaxmin, etaxmax);

TH1F *pt_mc_2lead = new TH1F("pt_mc_2lead", "pt_mc_2lead", 2000, 0, ptxmax);
TH1F *eta_mc_2lead = new TH1F("eta_mc_2lead", "eta_mc_2lead", 2000, etaxmin, etaxmax);
TH1F *phi_mc_2lead = new TH1F("phi_mc_2lead", "phi_mc_2lead", 2000, etaxmin, etaxmax);

TH1F *pt_mc_diff = new TH1F("pt_mc_diff", "pt_mc_diff", 2000, 0, ptxmax);

TH2F *ptdr_mc = new TH2F("ptdr_mc", "ptdr_mc", 2000, 0, ptxmax, 2000, 0, 0.5);
TH2F *ptdr_mc_cut = new TH2F("ptdr_mc_cut", "ptdr_mc_cut", 2000, 0, ptxmax, 2000, 0, 0.5);
TH2F *ptdr_reco = new TH2F("ptdr_reco", "ptdr_reco", 2000, 0, ptxmax, 2000, 0, 0.5);
TProfile *ptdrProf = new TProfile("ptdrProf", "ptdrProf", 2000, 0, 300000, 0, 0.5);


TH2F *recopt_eff = new TH2F("recopt_eff", "recopt_eff", 2000, 0, ptxmax, 2000, 0, 1.2);
TH2F *recoMpt_eff = new TH2F("recoMpt_eff", "recoMpt_eff", 2000, 0, ptxmax, 2000, 0, 1.2);
TProfile *recopt_prof = new TProfile("recopt_prof", "recopt_prof", 2000, 0, ptxmax, 0, 1.2 );
TProfile *recoMpt_prof = new TProfile("recoMpt_prof", "recoMpt_prof", 2000, 0, ptxmax, 0, 1.2 );


TH2F *recoeta_eff = new TH2F("recoeta_eff", "recoeta_eff", 2000, etaxmin, etaxmax, 2000, 0, 1.2);
TH2F *recoMeta_eff = new TH2F("recoMeta_eff", "recoMeta_eff", 2000, etaxmin, etaxmax, 2000, 0, 1.2);
TProfile *recoeta_prof = new TProfile("recoeta_prof", "recoeta_prof", 2000, etaxmin, etaxmax, 0, 1.2 );
TProfile *recoMeta_prof = new TProfile("recoMeta_prof", "recoMeta_prof", 2000, etaxmin, etaxmax, 0, 1.2 );


TH2F *recomu_eff = new TH2F("recomu_eff", "recomu_eff", mubin, muxmin, muxmax, 2000, 0, 1.2);
TH2F *recoMmu_eff = new TH2F("recoMmu_eff", "recoMmu_eff", mubin, muxmin, muxmax, 2000, 0, 1.2);
TProfile *recomu_prof = new TProfile("recomu_prof", "recomu_prof", mubin, muxmin, muxmax, 0, 1.2 );
TProfile *recoMmu_prof = new TProfile("recoMmu_prof", "recoMmu_prof", mubin, muxmin, muxmax, 0, 1.2 );

const int fcuts = 5;
const char* str_fhtcut[] = {"_0-003", "_003-005", "_005-01", "_01-03", "_03-1" };
TH1F *dr_mc[fcuts];
TH1F *dR_recMC[fcuts];
TH1F *dR_recMC_noDup[fcuts];
TH1F *dR_recMC_cl[fcuts];
TH1F *dR_TrkTruth1[fcuts];
TH1F *dR_TrkTruth1_ptcut[fcuts];
//~ TH1F *dR_TrkTruth1_etacut[fcuts];
//~ TH1F *dR_TrkCl1[fcuts];
//~ TH1F *dR_reco[fcuts];

for( int fcut=0;fcut<fcuts;fcut++){
	string fht_cut = str_fhtcut[fcut];
        dr_mc[fcut] = new TH1F(("dr_mc"+fht_cut).c_str(), ("dr_mc"+fht_cut).c_str(), 2000, 0, 0.5 );
	dR_recMC[fcut] = new TH1F(("dR_recMC"+fht_cut).c_str(), ("dR_recMC"+fht_cut).c_str(), 2000, 0, 0.5 );
	dR_recMC_noDup[fcut] = new TH1F(("dR_recMC_noDup"+fht_cut).c_str(), ("dR_recMC_noDup"+fht_cut).c_str(), 2000, 0, 0.5 );
	dR_recMC_cl[fcut] = new TH1F(("dR_recMC_cl"+fht_cut).c_str(), ("dR_recMC_cl"+fht_cut).c_str(), 2000, 0, 0.5 );
	dR_TrkTruth1[fcut] = new TH1F(("dR_TrkTruth1"+fht_cut).c_str(), ("dR_TrkTruth1"+fht_cut).c_str(), 4000, 0, 1 );
	dR_TrkTruth1_ptcut[fcut] = new TH1F(("dR_TrkTruth1_ptcut"+fht_cut).c_str(), ("dR_TrkTruth1_ptcut"+fht_cut).c_str(), 4000, 0, 1 );
	//~ dR_TrkTruth1_etacut[fcut] = new TH1F(("dR_TrkTruth1_etacut"+fht_cut).c_str(), ("dR_TrkTruth1_etacut"+fht_cut).c_str(), 4000, 0, 1 );
	//~ dR_TrkCl1[fcut] = new TH1F(("dR_TrkCl1"+fht_cut).c_str(), ("dR_TrkCl1"+fht_cut).c_str(), 2000, 0, 0.5 );
	//~ dR_reco[fcut] = new TH1F(("dR_reco"+fht_cut).c_str(), ("dR_reco"+fht_cut).c_str(), 2000, 0, 0.5 );
}

TH1F *pt_recMC = new TH1F("pt_recMC", "pt_recMC", 2000, 0, ptxmax);
TH1F *eta_recMC = new TH1F("eta_recMC", "eta_recMC", 2000, etaxmin, etaxmax);
TH1F *phi_recMC = new TH1F("phi_recMC", "phi_recMC", 2000, etaxmin, etaxmax);	

TH1F *pt_recMC_noDup = new TH1F("pt_recMC_noDup", "pt_recMC_noDup", 2000, 0, ptxmax);
TH1F *eta_recMC_noDup = new TH1F("eta_recMC_noDup", "eta_recMC_noDup", 2000, etaxmin, etaxmax);
TH1F *phi_recMC_noDup = new TH1F("phi_recMC_noDup", "phi_recMC_noDup", 2000, etaxmin, etaxmax);	

TH1F *pt_recMC_cl = new TH1F("pt_recMC_cl", "pt_recMC_cl", 2000, 0, ptxmax);
TH1F *eta_recMC_cl = new TH1F("eta_recMC_cl", "eta_recMC_cl", 2000, etaxmin, etaxmax);
TH1F *phi_recMC_cl = new TH1F("phi_recMC_cl", "phi_recMC_cl", 2000, etaxmin, etaxmax);	

TH1F *pt_reco_cl = new TH1F("pt_reco_cl", "pt_reco_cl", 2000, 0, ptxmax);
TH1F *eta_reco_cl = new TH1F("eta_reco_cl", "eta_reco_cl", 2000, etaxmin, etaxmax);



float El_mass = 0.511;
float El_m_Sq = El_mass*El_mass;


if (fChain == 0) return;
Long64_t nentries = fChain->GetEntriesFast();
cout << "nentries: " << nentries << endl;
Long64_t nbytes = 0, nb = 0;

for (Long64_t jentry=0; jentry<nentries;jentry++) {
	   //cout << "Start of event " << jentry << "*************" << endl;
	Long64_t ientry = LoadTree(jentry);
	if(ientry<0) break;

	nb += fChain->GetEntry(jentry);
	if(jentry % 1000 == 0) cout << jentry << endl;   
	//cout <<"actual mu: " << actualIntPerXing << " avg mu: " << averageIntPerXing << endl;
   	 
	float mcpt[2] = -99;
	float mceta[2] = -99;
 	float mcphi[2] = -99;   
	int j = 0;
	 
	if(mc_n < 2) continue;
	   
	for(Int_t i=0; i<mc_n; i++){
		int mcpdgid = (*mc_pdgId)[i];
	
		//cout << "pdgid: " << mcpdgid << endl;
		if(mcpdgid == 11 || mcpdgid == -11){
			mcpt[j] = (*mc_pt)[i];
			mceta[j] = (*mc_eta)[i];
			mcphi[j] = (*mc_phi)[i];						
		
			j++;
		} //End of if
	} //End of for
		

	float delR_truth = deltaR(mceta[0], mcphi[0], mceta[1], mcphi[1]);
	

	// Finding the leading and second leading pt lepton
	float mcptdiff = mcpt[0] - mcpt[1];
	float mcpt_lead = 0;
	float mcpt_2lead = 0;
	float mcpt_lead_inv = 0;
	float mceta_lead = -99;
	float mceta_2lead = -99;
	float mcphi_lead = -99;
	float mcphi_2lead = -99;
	
	if(mcptdiff > 0) { 
		mcpt_lead = mcpt[0];
		mceta_lead = mceta[0];
		mcphi_lead = mcphi[0];
		mcpt_2lead = mcpt[1];
		mceta_2lead = mceta[1];
		mcphi_2lead = mcphi[1];
	
		}
	else { 
		mcpt_lead = mcpt[1];
		mceta_lead = mceta[1];
		mcphi_lead = mcphi[1];
		mcpt_2lead = mcpt[0];
		mceta_2lead = mceta[0];
		mcphi_2lead = mcphi[0];
		}
		
	mcpt_lead_inv = 1/mcpt_lead;	 // Inverse of leading pt lepton
	
	float theta_lead = 2*atan(exp(-mceta_lead));	
	float theta_2lead = 2*atan(exp(-mceta_2lead));	
	
	ptdr_mc->Fill(mcpt_lead, delR_truth);
		
	if(mcpt_lead < 10000) continue;
	if(mcpt_2lead < 10000) continue;
	if(mceta_lead > 2.5 || mceta_2lead > 2.5) continue;	

	ptdr_mc_cut->Fill(mcpt_lead, delR_truth);
		
	pt_mc->Fill(mcpt_lead);
	eta_mc->Fill(mceta_lead);
	phi_mc->Fill(mcphi_lead);	
	
	pt_mc_2lead->Fill(mcpt_2lead);
	eta_mc_2lead->Fill(mceta_2lead);
	phi_mc_2lead->Fill(mcphi_2lead);
	pt_mc_diff->Fill(fabs(mcptdiff));	
		
		
	//Filling dR for different pT bins	
	float ptcut[] = {0, 0.003, 0.005, 0.01, 0.03, 1 };  // { inf-333GeV, 333-200GeV, 200-125GeV, 125-100GeV, 100-33GeV, 33-0GeV }
	float mcpt_inv_gev = mcpt_lead_inv*1000;
		
	//~ for(int kk=0; kk<fcuts; kk++) {
		//~ if(mcpt_inv_gev > ptcut[kk] && mcpt_inv_gev < ptcut[kk+1]) { dr_mc[kk]->Fill(delR_truth); }
	//~ }
		
	for(int el1=0; el1<el_n; el1++) {
		int nrefittrk =(*el_refittedTrack_n)[el1]; 
		cout << nrefittrk << endl;
	
		for(int tr1=0; tr1<nrefittrk; tr1++) {
			float trkqop = (*el_refittedTrack_qoverp)[el1][tr1];
			float trktheta = (*el_refittedTrack_theta)[el1][tr1];
			int refitauthor = (*el_refittedTrack_author)[el1][tr1];
			float trkthphi = (*el_refittedTrack_phi)[el1][tr1];
			
			float trkp = 1/trkqop;
			float trkpt = trkp*sin(trktheta);
			float trketa = -log(tan(trktheta/2));
			
			cout << "author: " << refitauthor << " pt: " << trkpt << " eta: " << trketa << endl;
			
			if(refitauthor == 4) cout << "GSF" << endl;
		}
	}

	

	
    //cout << "****** End of Event " << jentry << " *******\n" << endl;

} //End of jentry
   cout << "check2" << endl;
   
   if(writeOutput) {
   newfile->Write();
   delete newfile;
   }
   
   cout << "check3" << endl;
   

 } //End of Loop

 
 float deltaPhi(float phi1, float phi2) {
  float dphi = phi1 - phi2;
  while (dphi < -TMath::Pi())
    dphi += 2 * TMath::Pi();
  while (dphi > TMath::Pi())
    dphi -= 2 * TMath::Pi();
  return dphi;
}

float deltaVtx(float vx1, float vy1, float vz1, float vx2, float vy2, float vz2) {
  float dvx = vx1 - vx2;
  float dvy = vy1 - vy2;
  float dvz = vz1 - vz2;
  float dvr = TMath::Sqrt(dvx*dvx + dvy*dvy + dvz*dvz);
  return dvr;
}

float deltaR(float eta1, float phi1, float eta2, float phi2) {        
  float dphi = deltaPhi(phi1, phi2);
  float deta = fabs(eta1-eta2);
  float dr = TMath::Sqrt(dphi*dphi + deta*deta);
  return dr;
}