// 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 minJetPt = 90;
float maxJetPt = 110;
float coneSize = 0.7;

TFile *infile = new TFile("Wplus"+prefix+".root");
TTree *spartyTree = infile->Get("SpartyJet_Tree");
TFile *fout = new TFile("Wplus"+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;
// 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;   //!
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;
// Set branch addresses and branch pointers
spartyTree->SetBranchAddress("PythiaInput_N", &PythiaInput_N, &b_PythiaInput_N);
spartyTree->SetBranchAddress("PythiaInput_pdgId", &PythiaInput_pdgId, &b_PythiaInput_pdgId);
spartyTree->SetBranchAddress("PythiaInput_eta", &PythiaInput_eta, &b_PythiaInput_eta);
spartyTree->SetBranchAddress("PythiaInput_phi", &PythiaInput_phi, &b_PythiaInput_phi);
spartyTree->SetBranchAddress("PythiaInput_e", &PythiaInput_e, &b_PythiaInput_e);
spartyTree->SetBranchAddress("PythiaInput_mass", &PythiaInput_mass, &b_PythiaInput_mass);
spartyTree->SetBranchAddress("PythiaInput_pt", &PythiaInput_pt, &b_PythiaInput_pt);
spartyTree->SetBranchAddress("AntiKt7_N", &AntiKt7_N, &b_AntiKt7_N);
spartyTree->SetBranchAddress("AntiKt7_eta", &AntiKt7_eta, &b_AntiKt7_eta);
spartyTree->SetBranchAddress("AntiKt7_phi", &AntiKt7_phi, &b_AntiKt7_phi);
spartyTree->SetBranchAddress("AntiKt7_e", &AntiKt7_e, &b_AntiKt7_e);
spartyTree->SetBranchAddress("AntiKt7_mass", &AntiKt7_mass, &b_AntiKt7_mass);
spartyTree->SetBranchAddress("AntiKt7_pt", &AntiKt7_pt, &b_AntiKt7_pt);
spartyTree->SetBranchAddress("AntiKt7_ind", &AntiKt7_ind, &b_AntiKt7_ind);
spartyTree->SetBranchAddress("AntiKt7_numC", &AntiKt7_numC, &b_AntiKt7_numC);
//===================================================================


// Histograms
const int nRbins = 100;
const int nEtaBins = 10;
float binStep = 0.5;
float binStart = 0.0;
float binEnd;
stringstream hname,htitle;
TProfile* h_PtProfile[nEtaBins];
TProfile* h_diffPtProfile[nEtaBins];
for(int b=0; b< nEtaBins;b++)
{
	binEnd = binStart+binStep;
	// differential profile
	hname.str(""); htitle.str("");
	hname << "h_diffPtProfile_eta" << binStart << "-" << binEnd;
	htitle <<"diff p_{T} profile for Jets of " << binStart << "< |#eta| < " << binEnd << ";r/R;d#Psi/dr";
	h_diffPtProfile[b] = new TProfile(hname.str().c_str(),htitle.str().c_str(),100,0.0,1.0);
	// profile
	hname.str(""); htitle.str("");
	hname << "h_PtProfile_eta" << binStart << "-" << binEnd;
	htitle <<"p_{T} profile for Jets of " << binStart << "< |#eta| < " << binEnd << ";r/R;#Psi";
	h_PtProfile[b] = new TProfile(hname.str().c_str(),htitle.str().c_str(),100,0.0,1.0);
	binStart = binEnd;
}
TH1F * h_jetPt = new TH1F("h_jetPt","Jet p_{T} distribution;p_{T} [GeV]",100,0,300);
TH1F * h_jetEta = new TH1F("h_jetEta","Jet #eta distribution;#eta",101,-5,5);
TH1F * h_jetMass = new TH1F("h_jetMass","Jet mass distribution;jet mass [GeV/c^{2}]",80,0,40);
TH1F * h_jetConstitDelR = new TH1F("h_jetConstitDelR","Jet - Constituent delta R;#Delta R",200,0,5);

Int_t nEvents = spartyTree->GetEntries();
float energySum[nRbins];
int		nConstit[nRbins];
for(int e=0;e<nEvents;e++)
{
	if(e%1000==0) std::cout << "Events processed: " << e << std::endl;
	
	spartyTree->GetEntry(e);
	
	// Loop over jets
	for(int j=0;j<AntiKt7_N;j++)
	{
		// Jet cuts
		if((AntiKt7_pt->at(j) < minJetPt) || (AntiKt7_pt->at(j) > maxJetPt)) continue;
		int etaBin = (int)(TMath::Abs(AntiKt7_eta->at(j))/binStep);
		if(etaBin>9) continue;
		
		// Clear temp store of energy
		//memset(energySum,0.,nRbins);
		//memset(nConstit,0,nRbins);
		for(int r=0;r<nRbins;r++)
		{
			energySum[r] = 0.0;
			nConstit[r] = 0;
		}


		//  Loop over input
		for(int p=0;p<PythiaInput_N;p++)
		{
			// Match constituent with Jet
			if(AntiKt7_ind->at(p) == j)  // then input p was clustered into jet j
			{
				float delR = DeltaR(PythiaInput_eta->at(p),PythiaInput_phi->at(p),AntiKt7_eta->at(j),AntiKt7_phi->at(j))/coneSize;
				h_jetConstitDelR->Fill(delR);
				int rBin = (int)(delR*nRbins);
				if(rBin<nRbins)
				{
					nConstit[rBin]++;
					energySum[rBin] += (PythiaInput_pt->at(p)/AntiKt7_pt->at(j));
				}
			}
		}
		
		// Fill Jet histograms
		h_jetPt->Fill(AntiKt7_pt->at(j));
		h_jetEta->Fill(AntiKt7_eta->at(j));
		h_jetMass->Fill(AntiKt7_mass->at(j));
		float enclosedE = 0.0;
		for(int r=0;r<nRbins;r++)
		{
			float delR = ((1.0*r)/nRbins);
			enclosedE += energySum[r];
			h_PtProfile[etaBin]->Fill(delR,enclosedE);
			//  Only fill differential profile if
			//  at least one particle in bin
			if(nConstit[r]>0)
			{
				h_diffPtProfile[etaBin]->Fill(delR,energySum[r]);
			}
		}
	} // End jet loop
}// End event loop

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

}