// Helper function
float DeltaR(float eta1,float phi1,float eta2,float phi2)
{
	float deltaPhi = TMath::Abs(phi1-phi2);
	float deltaEta = eta1-eta2;
	if(deltaPhi > TMath::Pi())
	deltaPhi = TMath::TwoPi() - deltaPhi;
		return TMath::Sqrt(deltaEta*deltaEta + deltaPhi*deltaPhi);
}




void runAnalysis(TString prefix="Wenu")
{
// Configure Job
int nEvents = 100000;  // Total # of events
int totalPassed = 27000;  // Total # of events that pass cuts
float wTransMassCut = 20.;
float metCut = 30.;
float minElecPt = 20.;
float maxElecEta = 1.1;
float minJetPt = 15.;
float maxJetEta = 2.;
float deltaRcut = 0.1;

TFile *infile = new TFile(prefix+".root");
TTree *hadronTree = infile->Get("HadronTree");
TTree *caloTree 		= infile->Get("CaloTree");
TFile *fout = new TFile(prefix+".HIST.root","RECREATE");

// Branch config=====================================================
// Basically just cut and pasted from tree->MakeClass() output
Int_t           PythiaInput_N;
vector<float>   *PythiaInput_pdgId;
vector<float>   *PythiaInput_eta;
vector<float>   *PythiaInput_phi;
vector<float>   *PythiaInput_e;
vector<float>   *PythiaInput_mass;
vector<float>   *PythiaInput_pt;
Int_t           AntiKt7_N;
vector<float>   *AntiKt7_eta;
vector<float>   *AntiKt7_phi;
vector<float>   *AntiKt7_e;
vector<float>   *AntiKt7_mass;
vector<float>   *AntiKt7_pt;
vector<int>     *AntiKt7_ind;
vector<float>   *AntiKt7_numC;
Int_t            CaloInput_N;
vector<float>   *CaloInput_eta;
vector<float>   *CaloInput_phi;
vector<float>   *CaloInput_e;
vector<float>   *CaloInput_mass;
vector<float>   *CaloInput_pt;
Int_t            AntiKt7_Calo_N;
vector<float>   *AntiKt7_Calo_eta;
vector<float>   *AntiKt7_Calo_phi;
vector<float>   *AntiKt7_Calo_e;
vector<float>   *AntiKt7_Calo_mass;
vector<float>   *AntiKt7_Calo_pt;
vector<int>     *AntiKt7_Calo_ind;
vector<float>   *AntiKt7_Calo_numC;
// List of branches
TBranch        *b_PythiaInput_N;   //!
TBranch        *b_PythiaInput_pdgId;   //!
TBranch        *b_PythiaInput_eta;   //!
TBranch        *b_PythiaInput_phi;   //!
TBranch        *b_PythiaInput_e;   //!
TBranch        *b_PythiaInput_mass;   //!
TBranch        *b_PythiaInput_pt;   //!
TBranch        *b_AntiKt7_N;   //!
TBranch        *b_AntiKt7_eta;   //!
TBranch        *b_AntiKt7_phi;   //!
TBranch        *b_AntiKt7_e;   //!
TBranch        *b_AntiKt7_mass;   //!
TBranch        *b_AntiKt7_pt;   //!
TBranch        *b_AntiKt7_ind;   //!
TBranch        *b_AntiKt7_numC;   //!
TBranch        *b_CaloInput_N;   //!
TBranch        *b_CaloInput_eta;   //!
TBranch        *b_CaloInput_phi;   //!
TBranch        *b_CaloInput_e;   //!
TBranch        *b_CaloInput_mass;   //!
TBranch        *b_CaloInput_pt;   //!
TBranch        *b_AntiKt7_Calo_N;   //!
TBranch        *b_AntiKt7_Calo_eta;   //!
TBranch        *b_AntiKt7_Calo_phi;   //!
TBranch        *b_AntiKt7_Calo_e;   //!
TBranch        *b_AntiKt7_Calo_mass;   //!
TBranch        *b_AntiKt7_Calo_pt;   //!
TBranch        *b_AntiKt7_Calo_ind;   //!
TBranch        *b_AntiKt7_Calo_numC;   //!
PythiaInput_pdgId = 0;
PythiaInput_eta = 0;
PythiaInput_phi = 0;
PythiaInput_e = 0;
PythiaInput_mass = 0;
PythiaInput_pt = 0;
AntiKt7_eta = 0;
AntiKt7_phi = 0;
AntiKt7_e = 0;
AntiKt7_mass = 0;
AntiKt7_pt = 0;
AntiKt7_ind = 0;
AntiKt7_numC = 0;
CaloInput_eta = 0;
CaloInput_phi = 0;
CaloInput_e = 0;
CaloInput_mass = 0;
CaloInput_pt = 0;
AntiKt7_Calo_eta = 0;
AntiKt7_Calo_phi = 0;
AntiKt7_Calo_e = 0;
AntiKt7_Calo_mass = 0;
AntiKt7_Calo_pt = 0;
AntiKt7_Calo_ind = 0;
AntiKt7_Calo_numC = 0;
// Set branch addresses and branch pointers
hadronTree->SetBranchAddress("PythiaInput_N", &PythiaInput_N, &b_PythiaInput_N);
hadronTree->SetBranchAddress("PythiaInput_pdgId", &PythiaInput_pdgId, &b_PythiaInput_pdgId);
hadronTree->SetBranchAddress("PythiaInput_eta", &PythiaInput_eta, &b_PythiaInput_eta);
hadronTree->SetBranchAddress("PythiaInput_phi", &PythiaInput_phi, &b_PythiaInput_phi);
hadronTree->SetBranchAddress("PythiaInput_e", &PythiaInput_e, &b_PythiaInput_e);
hadronTree->SetBranchAddress("PythiaInput_mass", &PythiaInput_mass, &b_PythiaInput_mass);
hadronTree->SetBranchAddress("PythiaInput_pt", &PythiaInput_pt, &b_PythiaInput_pt);
hadronTree->SetBranchAddress("AntiKt7_N", &AntiKt7_N, &b_AntiKt7_N);
hadronTree->SetBranchAddress("AntiKt7_eta", &AntiKt7_eta, &b_AntiKt7_eta);
hadronTree->SetBranchAddress("AntiKt7_phi", &AntiKt7_phi, &b_AntiKt7_phi);
hadronTree->SetBranchAddress("AntiKt7_e", &AntiKt7_e, &b_AntiKt7_e);
hadronTree->SetBranchAddress("AntiKt7_mass", &AntiKt7_mass, &b_AntiKt7_mass);
hadronTree->SetBranchAddress("AntiKt7_pt", &AntiKt7_pt, &b_AntiKt7_pt);
hadronTree->SetBranchAddress("AntiKt7_ind", &AntiKt7_ind, &b_AntiKt7_ind);
hadronTree->SetBranchAddress("AntiKt7_numC", &AntiKt7_numC, &b_AntiKt7_numC);

caloTree->SetBranchAddress("CaloInput_N",    &CaloInput_N,    &b_CaloInput_N);
caloTree->SetBranchAddress("CaloInput_eta",  &CaloInput_eta,  &b_CaloInput_eta);
caloTree->SetBranchAddress("CaloInput_phi",  &CaloInput_phi,  &b_CaloInput_phi);
caloTree->SetBranchAddress("CaloInput_e",    &CaloInput_e,    &b_CaloInput_e);
caloTree->SetBranchAddress("CaloInput_mass", &CaloInput_mass, &b_CaloInput_mass);
caloTree->SetBranchAddress("CaloInput_pt",   &CaloInput_pt,   &b_CaloInput_pt);
caloTree->SetBranchAddress("AntiKt7_Calo_N",   &AntiKt7_Calo_N,   &b_AntiKt7_Calo_N);
caloTree->SetBranchAddress("AntiKt7_Calo_eta", &AntiKt7_Calo_eta, &b_AntiKt7_Calo_eta);
caloTree->SetBranchAddress("AntiKt7_Calo_phi", &AntiKt7_Calo_phi, &b_AntiKt7_Calo_phi);
caloTree->SetBranchAddress("AntiKt7_Calo_e",   &AntiKt7_Calo_e,   &b_AntiKt7_Calo_e);
caloTree->SetBranchAddress("AntiKt7_Calo_mass",&AntiKt7_Calo_mass,&b_AntiKt7_Calo_mass);
caloTree->SetBranchAddress("AntiKt7_Calo_pt",  &AntiKt7_Calo_pt,  &b_AntiKt7_Calo_pt);
caloTree->SetBranchAddress("AntiKt7_Calo_ind", &AntiKt7_Calo_ind, &b_AntiKt7_Calo_ind);
caloTree->SetBranchAddress("AntiKt7_Calo_numC",&AntiKt7_Calo_numC,&b_AntiKt7_Calo_numC);
//===================================================================


// Histograms
TH1I *h_nJets = new TH1I("h_nJets","Jet multiplicity;N jets (p_{T}>15GeV and #eta<2.0)",10,0,10);
TH1F *h_jetPt = new TH1F("h_jetPt","Jet p_{T} distribution;Jet p_{T} [GeV]",100,0,300);
TH1I *h_nTrueJets = new TH1I("h_nTrueJets","Hadron-level Jet multiplicity;N jets (p_{T}>15GeV and #eta<2.0)",10,0,10);
TH1F *h_trueJetPt = new TH1F("h_trueJetPt","Hadron-level Jet p_{T} distribution;Jet p_{T} [GeV]",100,0,300);
TH1F *h_TrueWm 		= new TH1F("h_TrueWm",  "True W mass;Mass [GeV]",200,0,200);
TH1F *h_TrueWm_t	= new TH1F("h_TrueWm_t","True W transverse mass;Transverse Mass [GeV]",200,0,200);
TH1F *h_TrueWpt	= new TH1F("h_TrueWpt","True W pt;p_{T} [GeV]",200,0,200);
TH1F *h_TrueWeta	= new TH1F("h_TrueWeta","True W Pseudorapidity;#eta",161,-4,4);
TH1F *h_TrueWy	= new TH1F("h_TrueWy","True W Rapidity;Rapidity",81,-4,4);
TH1F *h_Wm 		= new TH1F("h_Wm",  "Reconstructed W mass;Mass [GeV]",200,0,200);
TH1F *h_Wm_t	= new TH1F("h_Wm_t","Reconstructed W transverse mass;Transverse Mass [GeV]",200,0,200);
TH1F *h_Wpt	= new TH1F("h_Wpt","Reconstructed W pt;p_{T} [GeV]",200,0,200);
TH1F *h_Weta	= new TH1F("h_Weta","Reconstructed W Pseudorapidity;#eta",161,-4,4);
TH1F *h_Wy	= new TH1F("h_Wy","Reconstructed W Rapidity;Rapidity",81,-4,4);
TH1F *h_elecPt= new TH1F("h_elecPt","Reconstructed Electron pt;p_{T} [GeV]",200,0,200);
TH1F *h_elecPhi= new TH1F("h_elecPhi","Reconstructed Electron phi;#phi",100,-3.2,3.2);
TH1F *h_elecEta= new TH1F("h_elecEta","Reconstructed Electron eta;#eta",161,-4,4);
TH1F *h_elecY= new TH1F("h_elecY","Reconstructed Electron Rapidity;Rapidity",81,-4,4);
TH1F *h_missEt= new TH1F("h_missEt","Missing Transverse Energy;E_{Tmiss} [GeV]",200,0,200);
TH1F *h_trueElecPt= new TH1F("h_trueElecPt","True Electron pt;p_{T} [GeV]",200,0,200);
TH1F *h_trueElecPhi= new TH1F("h_trueElecPhi","True Electron phi;#phi",100,-3.2,3.2);
TH1F *h_trueElecEta= new TH1F("h_trueElecEta","True Electron eta;#eta",161,-4,4);
TH1F *h_trueElecY= new TH1F("h_trueElecY","True Electron Rapidity;Rapidity",81,-4,4);
TH1F *h_trueNuPt= new  TH1F("h_trueNuPt","True Neutrino pt;p_{T} [GeV]",200,0,200);
TH1F *h_trueNuPhi= new TH1F("h_trueNuPhi","True Neutrino phi;#phi",100,-3.2,3.2);
TH1F *h_trueNuEta= new TH1F("h_trueNuEta","True Neutrino eta;#eta",161,-4,4);
TH1F *h_trueNuY=   new TH1F("h_trueNuY",  "True Neutrino Rapidity;Rapidity",161,-4,4);
TH1F *h_delPhi_nuMET = new TH1F("h_delPhi_nuMET","Phi separation between E_{Tmiss} and true neutrino;#Delta#phi",241,-6,6);
TH1F *h_deltaR_trueElecElec = new TH1F("h_deltaR_trueElecElec","Delta R between electron and true electron;#DeltaR",241,-6,6);
TH2F *h_trueElecEtatrueWeta_Wpt = new TH2F("h_trueElecEtatrueWeta_Wpt","Electron Eta - W Eta;W p_{T};#Delta #eta",200,0,200,161,-1,1);
TH1F *h_jetResponse = new TH1F("h_jetResponse","Jet Reponse;Calo p_{T}/True p_{T}",100,-10,10);
TH1F *h_metResponse = new TH1F("h_metResponse","E_{Tmiss} Response;E_{Tmiss}/neutrino E_{T}",100,-10,10);
TProfile *h_jetResponse_pt = new TProfile("h_jetResponse_pt","Jet Reponse vs Pt;True p_{T}[GeV];Calo/True",60,0,300);
TProfile *h_metResponse_et = new TProfile("h_metResponse_et","E_{Tmiss} Reponse vs neutrino Pt;neutrino E_{T}[GeV];MET/Neutrino E_{T}",200,0,200);
TH1F *h_sumJetPt = new TH1F("h_sumJetPt","Vector sum of Jet Pt;p_{T} [GeV]",100,0,300);
TH2F *h_sumJetPt_Wpt = new TH2F("h_sumJetPt_Wpt","Vector sum of Jet Pt vs W pt;#Sigma Jet p_{T} [GeV];W p_{T} [GeV]",100,0,300,100,0,300);
TH1F *h_sumTrueJetPt = new TH1F("h_sumTrueJetPt","Vector sum of Hadron-level Jet Pt;p_{T} [GeV]",100,0,300);
TH2F *h_sumTrueJetPt_TrueWpt = new TH2F("h_sumTrueJetPt_TrueWpt","Vector sum of Hadron-level Jet Pt vs W pt;#Sigma Jet p_{T} [GeV];W p_{T} [GeV]",100,0,300,100,0,300);
TH2F *h_Wmt_Wpt = new TH2F("h_Wmt_Wpt","W transverse mass vs W pt;W Transverse mass[GeV];W p_{T} [GeV]",100,0,300,100,0,300);

//  Find the total number of events to run
Int_t nEventsHadron = hadronTree->GetEntries();
Int_t nEventsCalo		 	= caloTree->GetEntries();
if(nEvents > nEventsHadron)
	nEvents = nEventsHadron;
if(nEvents > nEventsCalo)
	nEvents = nEventsCalo;

// Some data containers
TLorentzVector elecTrue = TLorentzVector();
TLorentzVector nuTrue 	= TLorentzVector();
TLorentzVector electron = TLorentzVector();
TLorentzVector neutrino = TLorentzVector();
TLorentzVector metVec 	= TLorentzVector();
TLorentzVector caloJet	= TLorentzVector();

// Useful flags and indexes 
bool elecFound = false;
bool nuFound = false;
int elecJetIndex = -1;
//for(Int_t e =0;e<nEvents;e++)
Int_t e=0;
int nPassed=0;
while(e<nEvents && nPassed<totalPassed)
{
	if(e%1000==0) std::cout << "Events processed: " << e << std::endl;
	
	// Reset some variables
	elecFound = false;
	nuFound = false;
	elecJetIndex = -1;
	
	// Load entries in both trees
	hadronTree->GetEntry(e);
	caloTree->GetEntry(e);
	
	// Find true electron and neutrino from pythia ouptut
	for(int p=0;p<PythiaInput_N;p++)
	{
		if(TMath::Abs(PythiaInput_pdgId->at(p))==11 && !elecFound)
		{
			//elecTrue.SetPtEtaPhiE(PythiaInput_pt->at(p),PythiaInput_eta->at(p),PythiaInput_phi->at(p),PythiaInput_e->at(p));
			elecTrue.SetPtEtaPhiM(PythiaInput_pt->at(p),PythiaInput_eta->at(p),PythiaInput_phi->at(p),0.000511);
			elecFound = true;
		}
		else if(TMath::Abs(PythiaInput_pdgId->at(p))==12 && !nuFound)
		{
			//nuTrue.SetPtEtaPhiE(PythiaInput_pt->at(p),PythiaInput_eta->at(p),PythiaInput_phi->at(p),PythiaInput_e->at(p));
			nuTrue.SetPtEtaPhiM(PythiaInput_pt->at(p),PythiaInput_eta->at(p),PythiaInput_phi->at(p),0.0);
			nuFound = true;
		}
		if (elecFound && nuFound) break;
	}
	
	// Compute Missing ET from Calorimeter "towers"
	double SumEx=0., SumEy=0.;
	for(int c=0;c<CaloInput_N;c++)
	{
		SumEx += CaloInput_pt->at(c)*TMath::Cos(CaloInput_phi->at(c));
		SumEy += CaloInput_pt->at(c)*TMath::Sin(CaloInput_phi->at(c));
	}
	//metVec.SetPxPyPzE(-1.*SumEx,-1.*SumEy,0.0,TMath::Sqrt(SumEx*SumEx+SumEy*SumEy));
	metVec.SetXYZM(-1.*SumEx,-1.*SumEy,0.0,0.0);

	
	// Loop over jets to find electron jet
	for(int cj=0;cj<AntiKt7_Calo_N;cj++)
	{
		if (elecFound)
		{
			caloJet.SetPtEtaPhiM(AntiKt7_Calo_pt->at(cj),AntiKt7_Calo_eta->at(cj),AntiKt7_Calo_phi->at(cj),0.000511);
			if(caloJet.DeltaR(elecTrue) < deltaRcut) 
			{
				electron = caloJet;
				elecJetIndex = cj;
			}
		}								
	}
	
	// Reconstruct  W
	// ---TrueW
	TLorentzVector trueW = elecTrue+nuTrue;
	// ---W
	TLorentzVector neutrino;
	neutrino.SetPtEtaPhiM(metVec.Pt(),electron.Eta(),metVec.Phi(),0.);
	TLorentzVector W = electron+neutrino;
	// Three definitions of Transverse mass of W
	//Double_t mW_t = TMath::Sqrt(2*electron.Et()*metVec.Et()*(1-TMath::Cos(electron.DeltaPhi(metVec))));
	//Double_t mW_t = W.Mt();
	Double_t mW_t = TMath::Sqrt(electron.Et()*electron.Et() + 
															metVec.Et()*metVec.Et() 
															- 2*electron.Et()*metVec.Et()*TMath::Cos(electron.DeltaPhi(metVec)));

	e++;
	// =================Make Event Selection==================
	if(
	// Electron cuts
	elecJetIndex == -1 || electron.Pt() < minElecPt || TMath::Abs(electron.Eta()) > maxElecEta 
	// Missing Et cut
	|| metVec.Et() < metCut 
	// Transverse M of W boson cut
	|| mW_t < wTransMassCut
	) continue;
	// =======================================================
	nPassed++;
	
  // Fill W Histograms
	h_TrueWm->Fill(trueW.M());
	h_TrueWm_t->Fill(trueW.Mt());
	h_TrueWpt->Fill(trueW.Pt());
	h_TrueWeta->Fill(trueW.Eta());
	h_TrueWy->Fill(trueW.Rapidity());
	
	h_Wm->Fill(W.M());
	h_Wm_t->Fill(mW_t);
	h_Wpt->Fill(W.Pt());
	h_Weta->Fill(W.Eta());
	h_Wy->Fill(W.Rapidity());

	// Fill lepton histograms
	h_trueElecPt->Fill(elecTrue.Pt());
	h_trueElecPhi->Fill(elecTrue.Phi());
	h_trueElecEta->Fill(elecTrue.Eta());
	h_trueElecY->Fill(elecTrue.Rapidity());
	
	h_trueNuPt->Fill(nuTrue.Pt());
	h_trueNuPhi->Fill(nuTrue.Phi());
	h_trueNuEta->Fill(nuTrue.Eta());
	h_trueNuY->Fill(nuTrue.Rapidity());
	
	h_elecPt->Fill(electron.Pt());
	h_elecPhi->Fill(electron.Phi());
	h_elecEta->Fill(electron.Eta());
	h_elecY->Fill(electron.Rapidity());
	
	h_missEt->Fill(metVec.Et());
	h_delPhi_nuMET->Fill(metVec.DeltaPhi(nuTrue));
	h_deltaR_trueElecElec->Fill(electron.DeltaR(elecTrue));
	h_trueElecEtatrueWeta_Wpt->Fill(trueW.Pt(),elecTrue.Eta()-trueW.Eta());
	h_metResponse->Fill(metVec.Et()/nuTrue.Et());
	h_metResponse_et->Fill(nuTrue.Et(),metVec.Et()/nuTrue.Et());
	
	// Fill TrueJet histograms
	double sumTrueJetPx=0.,sumTrueJetPy=0.;
	int nTrueJetsPassed = 0;
	for(int tj=0;tj<AntiKt7_N;tj++)
	{
		// Jet cuts
		if((AntiKt7_pt->at(tj) < minJetPt) || TMath::Abs(AntiKt7_eta->at(tj)) > maxJetEta) continue;
		sumTrueJetPx += AntiKt7_pt->at(tj)*TMath::Cos(AntiKt7_phi->at(tj));
		sumTrueJetPy += AntiKt7_pt->at(tj)*TMath::Sin(AntiKt7_phi->at(tj));

		h_trueJetPt->Fill(AntiKt7_pt->at(tj));
		nTrueJetsPassed++;
	}
	h_sumTrueJetPt->Fill(TMath::Sqrt(sumTrueJetPx*sumTrueJetPx+sumTrueJetPy*sumTrueJetPy));
	h_sumTrueJetPt_TrueWpt->Fill(TMath::Sqrt(sumTrueJetPx*sumTrueJetPx+sumTrueJetPy*sumTrueJetPy),trueW.Pt());
	h_nTrueJets->Fill(nTrueJetsPassed);

	// Fill Jet histograms
	double sumJetPx=0.,sumJetPy=0.;
	int nJetsPassed = 0;
	for(int cj=0;cj<AntiKt7_Calo_N;cj++)
	{
		// Don't use jet from electron
		if(elecJetIndex==cj) continue;
		// Jet cuts
		if((AntiKt7_Calo_pt->at(cj) < minJetPt) || TMath::Abs(AntiKt7_Calo_eta->at(cj)) > maxJetEta) continue;
		//  Loop over true jets to match with calo jets
		for(int tj=0;tj<AntiKt7_N;tj++)
		{
			float delR = DeltaR(AntiKt7_eta->at(tj),AntiKt7_phi->at(tj),AntiKt7_Calo_eta->at(cj),AntiKt7_Calo_phi->at(cj));
			if(delR < deltaRcut)
			{
				h_jetResponse->Fill(AntiKt7_Calo_pt->at(cj)/AntiKt7_pt->at(tj));
				h_jetResponse_pt->Fill( AntiKt7_Calo_pt->at(cj),AntiKt7_Calo_pt->at(cj)/AntiKt7_pt->at(tj));
			}
		}
	
		// Determine net jet pT
		sumJetPx += AntiKt7_Calo_pt->at(cj)*TMath::Cos(AntiKt7_Calo_phi->at(cj));
		sumJetPy += AntiKt7_Calo_pt->at(cj)*TMath::Sin(AntiKt7_Calo_phi->at(cj));
		
		h_jetPt->Fill(AntiKt7_Calo_pt->at(cj));
		nJetsPassed++;
	}
	h_sumJetPt->Fill(TMath::Sqrt(sumJetPx*sumJetPx+sumJetPy*sumJetPy));
	h_sumJetPt_Wpt->Fill(TMath::Sqrt(sumJetPx*sumJetPx+sumJetPy*sumJetPy),W.Pt());
	h_nJets->Fill(nJetsPassed);
	
	h_Wmt_Wpt->Fill(mW_t,W.Pt());
	
}// End event loop

fout->Write();
fout->Close();

}