#include "SBAnalyzer.h"

void SBAnalyzer::BVertexProbability()
{
	// MC (Signal) //
	std::vector<TH1D*> MCBVtxPHistograms;
	const int NBVtxP = 10;
	const double BVtxPStep = 0.01;
	const double BVtxPStart = 0.00;

	for (int iV=0; iV<NBVtxP; iV++) {
		double BVP = BVtxPStart+(iV*BVtxPStep);
		std::string MCBVtxPHistName = "MCBVtxP"+to_string(BVP);
		std::string MCBVtxPHistTitle = "B-Vertex Probability > "+to_string(BVP);
		TH1D* MCBVtxPHistogram = new TH1D(MCBVtxPHistName.c_str(), MCBVtxPHistTitle.c_str(), 100, 0.0, 1.0);
		MCBVtxPHistograms.push_back(MCBVtxPHistogram);
	}

	int AllMCEntries = SInputTree->GetEntries();
	std::cout << "Number of Entries (MC Sample) = " << AllMCEntries << endl;
	for(int iEn=0; iEn<AllMCEntries; iEn++) { 		// <---<< Main loop over Input Tree entries
		if(iEn%100==0)
			std::cout << "\r" << (double)iEn/(double)AllMCEntries*100 << "\% Processed" << flush;

		ClearAllMCVariables();
		SInputTree->GetEntry(iEn);

		for (int Index=0; Index<(int)SIn_EventNumber->size(); Index++) { 		// <---<< Loop over combinations in an entry

			if (SIn_TMuonPt->at(Index)>TMuonPtCut) { 		// <---<< Triggering Muon Pt

				if (abs(SIn_TMuonEta->at(Index))<TMuonEtaCut) { 		// <---<< Triggering Muon Eta (The trigger has L1 seed with eta<1.5)

					if ((abs(SIn_TMuondxy->at(Index))/SIn_TMuondxyError->at(Index))>TMuonTIPSCut) { 		// <---<< Triggering Muon TIPS

						if (SIn_TMuonIsMediumMuon->at(Index)) { 		// <---<< Triggering Muon ID

							if (SIn_SLeptonPt->at(Index)>SLeptonPtCut) { 		// <---<< Non-Triggering Lepton Pt

								if (abs(SIn_SLeptonEta->at(Index))<SLeptonEtaCut) { 		// <---<< Non-Triggering Lepton Eta

									if ((abs(SIn_SLeptondxy->at(Index))/SIn_SLeptondxyError->at(Index))>SLeptonTIPSCut) { 		// <---<< Non-Triggering Lepton TIPS

										if (SIn_SLeptonNoIsoPFMVAElectronWP90V2->at(Index)) { 		// <---<< Non-Triggering Lepton ID
										//if (SIn_SLeptonIsMediumMuon->at(Index)) {

											for (int jV=0; jV<NBVtxP; jV++) {
												double BVP = BVtxPStart+(jV*BVtxPStep);
												if (SIn_BVtxProbability->at(Index)>=BVP) {
													MCBVtxPHistograms[jV]->Fill(SIn_BVtxProbability->at(Index));
												}
											}
										} 		// <---<< End of "Non-Triggering Lepton ID"
									} 		// <---<< End of "Non-Triggering Lepton TIPS"
								} 		// <---<< End of "Non-Triggering Lepton Eta"
							} 		// <---<< End of "Non-Triggering Lepton Pt"
						} 		// <---<< End of "Triggering Muon ID"
					} 		// <---<< End of "Triggering Muon TIPS"
				} 		// <---<< End of "Triggering Muon Eta"
			} 		// <---<< End of "Triggering Muon Pt"
		} 		// <---<< End of "Combination Loop"
	} 		// <---<< End of "Main Loop"

	// Data (Background) //
	std::vector<TH1D*> DBVtxPHistograms;
	for (int iV=0; iV<NBVtxP; iV++) {
		double BVP = BVtxPStart+(iV*BVtxPStep);
		std::string DBVtxPHistName = "DBVtxP"+to_string(BVP);
		std::string DBVtxPHistTitle = "B-Vertex Probability > "+to_string(BVP);
		TH1D* DBVtxPHistogram = new TH1D(DBVtxPHistName.c_str(), DBVtxPHistTitle.c_str(), 100, 0.0, 1.0);
		DBVtxPHistograms.push_back(DBVtxPHistogram);
	}

	int AllDEntries = BInputTree->GetEntries();
	std::cout << "\nNumber of Entries (Data Sample) = " << AllDEntries << endl;
	for(int iEn=0; iEn<AllDEntries; iEn++) { 		// <---<< Main loop over Input Tree entries
		if(iEn%100==0)
			std::cout << "\r" << (double)iEn/(double)AllDEntries*100 << "\% Processed" << flush;

		ClearAllDataVariables();
		BInputTree->GetEntry(iEn);

		for (int Index=0; Index<(int)BIn_EventNumber->size(); Index++) { 		// <---<< Loop over combinations in an entry

			if (BIn_TMuonPt->at(Index)>TMuonPtCut) { 		// <---<< Triggering Muon Pt

				if (abs(BIn_TMuonEta->at(Index))<TMuonEtaCut) { 		// <---<< Triggering Muon Eta (The trigger has L1 seed with eta<1.5)

					if ((abs(BIn_TMuondxy->at(Index))/BIn_TMuondxyError->at(Index))>TMuonTIPSCut) { 		// <---<< Triggering Muon TIPS

						if (BIn_TMuonIsMediumMuon->at(Index)) { 		// <---<< Triggering Muon ID

							if (BIn_SLeptonPt->at(Index)>SLeptonPtCut) { 		// <---<< Non-Triggering Lepton Pt

								if (abs(BIn_SLeptonEta->at(Index))<SLeptonEtaCut) { 		// <---<< Non-Triggering Lepton Eta

									if ((abs(BIn_SLeptondxy->at(Index))/BIn_SLeptondxyError->at(Index))>SLeptonTIPSCut) { 		// <---<< Non-Triggering Lepton TIPS

										if (BIn_SLeptonNoIsoPFMVAElectronWP90V2->at(Index)) { 		// <---<< Non-Triggering Lepton ID
										//if (BIn_SLeptonIsMediumMuon->at(Index)) {

											if ((BIn_PsiMass->at(Index)>2.80 && BIn_PsiMass->at(Index)<2.95) || (BIn_PsiMass->at(Index)>3.25 && BIn_PsiMass->at(Index)<3.40)) { 		// <---<< Side Bands

												for (int jV=0; jV<NBVtxP; jV++) {
													double BVP = BVtxPStart+(jV*BVtxPStep);
													if (BIn_BVtxProbability->at(Index)>=BVP) {
														DBVtxPHistograms[jV]->Fill(BIn_BVtxProbability->at(Index));
													}
												}
											} 		// <---<< End of "Side Bands"
										} 		// <---<< End of "Non-Triggering Lepton ID"
									} 		// <---<< End of "Non-Triggering Lepton TIPS"
								} 		// <---<< End of "Non-Triggering Lepton Eta"
							} 		// <---<< End of "Non-Triggering Lepton Pt"
						} 		// <---<< End of "Triggering Muon ID"
					} 		// <---<< End of "Triggering Muon TIPS"
				} 		// <---<< End of "Triggering Muon Eta"
			} 		// <---<< End of "Triggering Muon Pt"
		} 		// <---<< End of "Combination Loop"
	} 		// <---<< End of "Main Loop"

	THStack* BVtxPHStack = new THStack("BVtxPHStack","B-Vertex Probability Histogram Stack");

	// Rescaling Histograms //
	double MCScaleFactor = MCBVtxPHistograms[0]->GetEntries();
	for (int ih=0; ih<(int)MCBVtxPHistograms.size(); ih++) {
		MCBVtxPHistograms[ih]->Scale(1/MCScaleFactor);
		int Color = ih+1;
		if (Color==10)
			Color = 46;
		MCBVtxPHistograms[ih]->SetLineColor(Color);
		MCBVtxPHistograms[ih]->SetFillColor(Color);
		MCBVtxPHistograms[ih]->SetFillStyle(3001);
		BVtxPHStack->Add(MCBVtxPHistograms[ih]);
	}

	double DScaleFactor = DBVtxPHistograms[0]->GetEntries();
	for (int ih=0; ih<(int)DBVtxPHistograms.size(); ih++) {
		DBVtxPHistograms[ih]->Scale(1/DScaleFactor);
		int Color = ih+1;
		if (Color==10)
			Color = 46;
		DBVtxPHistograms[ih]->SetLineColor(Color);
		DBVtxPHistograms[ih]->SetFillColor(Color);
		DBVtxPHistograms[ih]->SetFillStyle(3001);
		BVtxPHStack->Add(DBVtxPHistograms[ih]);
	}

	// Creating Canvases //
	double WW = 1930;
	double WH = 1450;

	TCanvas* BVtxProb_canv = new TCanvas("BVertexProbabilities","B-Vertex Probabilities");
		BVtxProb_canv->SetWindowSize(WW, WH);
		gStyle->SetOptStat(0);
			gPad-> SetLogy();
			gPad->SetGridx();
			gPad->SetGridy();
			MCBVtxPHistograms[0]->SetTitle("B-Vertex Probabilities (MC vs. Data)");
			MCBVtxPHistograms[0]->GetXaxis()->SetTitle("VtxProb_{B}");
			MCBVtxPHistograms[0]->GetXaxis()->SetRangeUser(-0.1, 1.0);
			MCBVtxPHistograms[0]->GetXaxis()->SetNdivisions(11);
			MCBVtxPHistograms[0]->GetYaxis()->SetRangeUser(1e-6,(2*max(MCBVtxPHistograms[0]->GetMaximum(),DBVtxPHistograms[0]->GetMaximum())));
			BVtxPHStack->Draw();

	std::cout << endl;
}