//Program to read data from root file and make invarient mass hist
// Based on: https://mhance.scipp.ucsc.edu/analysisintro.php

//V2: Now possible to open second tree and add to hists. also added pos/neg hists and posVneg hists
//V3: Invariant mass plot works!!!
//V4: PosNeg plots, pX, pY, pZ
//V5: Plot Total muon, divide canvas

//V7, 10/1: Adding output of PxPyPzE from 4Vec, got rid of DoBothTrees. Added momentum hists.
	// https://root.cern.ch/root/roottalk/roottalk99/2970.html

# include <iostream>
# include <fstream>
# include <math.h>
# include <iomanip>
# include <cmath>
# include <stdlib.h>
# include <cstdlib>
//# include <fstream.h>
# include <string.h>
# include <string>
# include <cstring> 
# include <vector>
//# include <dos.h> //For Sleep() 
#include <typeinfo>

# include "TROOT.h"
# include "TFile.h"
# include "TTree.h"
# include "TBrowser.h"
# include "TH1.h"
# include "TH2.h"
# include "TH3.h"
# include "TRandom.h"

# include <TLorentzVector.h>


using namespace std;
int main(){
	
	/////Controls/////

	//int NEntries = 469384; //469384 for data.root. Click on tree and read hist to find these.
	
	char rootFilePrefix[50] = ""; //Prefix to root file, this is the path to it - all the way up until its name.
	char rootFileName[50] = "mhance.root"; //The root file name itself - this is where the tree is.
	char TreeName1[50] = "OutputTree;34";
	
	char rootOutputFileName1[50] = "pTLepHist.root"; //Name of the output file for nLep hist
	char rootOutputFileName2[50] = "etaLepHist.root"; //Name of the output file for nLep hist
	char rootOutputFileName3[50] = "phiLepHist.root"; //Name of the output file for nLep hist
	
	char rootOutputFileName4[50] = "pTLepHistPos.root"; //Name of the output file for nLep hist
	char rootOutputFileName5[50] = "nLepHistPos.root"; //Name of the output file for nLep hist
	char rootOutputFileName6[50] = "nLepHistPos.root"; //Name of the output file for nLep hist
	
	char rootOutputFileName7[50] = "pTLepHistNeg.root"; //Name of the output file for nLep hist
	char rootOutputFileName8[50] = "nLepHistNeg.root"; //Name of the output file for nLep hist
	char rootOutputFileName9[50] = "nLepHistNeg.root"; //Name of the output file for nLep hist
	
	char rootOutputFileName10[50] = "nLepHist.root"; //Name of the output file for nLep hist
	char rootOutputFileName11[50] = "InvarientMass.root"; //Name of the output file for Invarient Mass hist
	
	char rootOutputFileName12[50] = "pTLepHistPosVNeg.root"; //Name of the output file for nLep hist
	char rootOutputFileName13[50] = "nLepHistPosVNeg.root"; //Name of the output file for nLep hist
	char rootOutputFileName14[50] = "nLepHistPosVNeg.root"; //Name of the output file for nLep hist
	
	char rootOutputFileName15[50] = "PosMomentumPrime.root"; //Name of the output file for nLep hist
	char rootOutputFileName16[50] = "NegMomentumPrime.root"; //Name of the output file for nLep hist
	char rootOutputFileName17[50] = "PosVNegMomentumPrime.root"; //Name of the output file for nLep hist
	char rootOutputFileName18[50] = "TotMomentumPrime.root"; //Name of the output file for nLep hist
	char rootOutputFileName19[50] = "TotMomentum.root"; //Name of the output file for nLep hist
	char rootOutputFileName20[50] = "PosMomentum.root"; //Name of the output file for nLep hist
	char rootOutputFileName21[50] = "NegMomentum.root"; //Name of the output file for nLep hist
	char rootOutputFileName22[50] = "PosVNegMomentum.root"; //Name of the output file for nLep hist
	
	double xMin = -100;
	double xMax = 100;
	double yMin = -300;
	double yMax = 300;
	
	const int xNChannels = 500; //Number of channels for root hist on x axis
	const int yNChannels = 500; //Number of channels for root hist on y axis 
	
	bool SaveAsRootFile = false; //If true, will save canvases to root file instead of printing to screen.
	bool IncludeOnlyMuons = true; //If true, will include flavLep = +-1 as well. ie flavLep<0, >0
	
	///////////////////////////////////////
	
	/////Varibles (Pre-stated)/////
	
	char rootFilePath1[100]; 
	
	int flavLep[3], nLeps35;
	int NumOfPosNegPairs = 0;
	
	float pTLep[3], etaLep[3], phiLep[3];
	
	double pTLepPos, etaLepPos, phiLepPos, pTLepNeg, etaLepNeg, phiLepNeg;
	
	double RapidityPosPrime, RapidityNegPrime, RapidityPos, RapidityNeg;
	
	double InvarientMass;
	
	double c = 1; // 2.99792458e8;
	
	ofstream outFile1;
	
	TH1F *histpTLep = new TH1F("histpTLep", "pTLep hist", xNChannels, xMin, xMax);
	TH1F *histetaLep = new TH1F("histetaLep", "etaLep hist", xNChannels, xMin, xMax);
	TH1F *histphiLep = new TH1F("histphiLep", "phiLep hist", xNChannels, xMin, xMax);
	
	TH1F *histpTLepPos = new TH1F("histpTLepPos", "pTLepPos hist", xNChannels, xMin, xMax);
	TH1F *histetaLepPos = new TH1F("histetaLepPos", "etaLepPos hist", xNChannels, xMin, xMax);
	TH1F *histphiLepPos = new TH1F("histphiLepPos", "phiLepPos hist", xNChannels, xMin, xMax);
	
	TH1F *histpTLepNeg = new TH1F("histpTLepNeg", "pTLepNeg hist", xNChannels, xMin, xMax);
	TH1F *histetaLepNeg = new TH1F("histetaLepNeg", "etaLepNeg hist", xNChannels, xMin, xMax);
	TH1F *histphiLepNeg = new TH1F("histphiLepNeg", "phiLepNeg hist", xNChannels, xMin, xMax);
	
	TH1F *histnLep = new TH1F("histnLep", "nLep hist", 10, -5, 5);
	
	TH1F *histM = new TH1F("histM", "Invarient Mass hist", xNChannels, -200, 200);
	TH1F *histM2 = new TH1F("histM2", "Invarient Mass 2 hist", xNChannels, -200, 200);
	
	TH2F *histPosVNegpTLep = new TH2F("histPosVNegpTLep", "PosVNegpTLep hist", xNChannels, xMin, xMax, yNChannels, yMin, yMax);
	TH2F *histPosVNegetaLep = new TH2F("histPosVNegetaLep", "PosVNegetaLep hist", xNChannels, xMin, xMax, yNChannels, yMin, yMax);
	TH2F *histPosVNegphiLep = new TH2F("histPosVNegphiLep", "PosVNegphiLep hist", xNChannels, xMin, xMax, yNChannels, yMin, yMax);
	
	TH1F *histPosPxPrime = new TH1F("histPosPxPrime", "PosPxPrime hist", xNChannels, xMin, xMax);
	TH1F *histPosPyPrime = new TH1F("histPosPyPrime", "PosPyPrime hist", xNChannels, xMin, xMax);
	TH1F *histPosPzPrime = new TH1F("histPosPzPrime", "PosPzPrime hist", xNChannels, xMin, xMax);
	TH1F *histPosEPrime = new TH1F("histPosEPrime", "PosEPrime hist", xNChannels, xMin, xMax);
	TH1F *histNegPxPrime = new TH1F("histNegPxPrime", "NegPxPrime hist", xNChannels, xMin, xMax);
	TH1F *histNegPyPrime = new TH1F("histNegPyPrime", "NegPyPrime hist", xNChannels, xMin, xMax);
	TH1F *histNegPzPrime = new TH1F("histNegPzPrime", "NegPzPrime hist", xNChannels, xMin, xMax);
	TH1F *histNegEPrime = new TH1F("histNegEPrime", "NegEPrime hist", xNChannels, xMin, xMax);
	
	TH2F *histPosVNegPxPrime = new TH2F("histPosVNegPxPrime", "PosVNegPxPrime hist", xNChannels, xMin, xMax, yNChannels, yMin, yMax);
	TH2F *histPosVNegPyPrime = new TH2F("histPosVNegPyPrime", "PosVNegPyPrime hist", xNChannels, xMin, xMax, yNChannels, yMin, yMax);
	TH2F *histPosVNegPzPrime = new TH2F("histPosVNegPzPrime", "PosVNegPzPrime hist", xNChannels, xMin, xMax, yNChannels, yMin, yMax);
	TH2F *histPosVNegEPrime = new TH2F("histPosVNegEPrime", "PosVNegEPrime hist", xNChannels, xMin, xMax, yNChannels, yMin, yMax);
	
	TH1F *histTotPxPrime = new TH1F("histTotPxPrime", "TotPxPrime hist", xNChannels, xMin, xMax);
	TH1F *histTotPyPrime = new TH1F("histTotPyPrime", "TotPyPrime hist", xNChannels, xMin, xMax);
	TH1F *histTotPzPrime = new TH1F("histTotPzPrime", "TotPzPrime hist", xNChannels, xMin, xMax);
	TH1F *histTotEPrime = new TH1F("histTotEPrime", "TotEPrime hist", xNChannels, xMin, xMax);
	
	TH1F *histTotPx = new TH1F("histTotPx", "TotPx hist", xNChannels, xMin, xMax);
	TH1F *histTotPy = new TH1F("histTotPy", "TotPy hist", xNChannels, xMin, xMax);
	TH1F *histTotPz = new TH1F("histTotPz", "TotPz hist", xNChannels, xMin, xMax);
	TH1F *histTotE = new TH1F("histTotE", "TotE hist", xNChannels, xMin, xMax);
	
	TH1F *histPosPx = new TH1F("histPosPx", "PosPx hist", xNChannels, xMin, xMax);
	TH1F *histPosPy = new TH1F("histPosPy", "PosPy hist", xNChannels, xMin, xMax);
	TH1F *histPosPz = new TH1F("histPosPz", "PosPz hist", xNChannels, xMin, xMax);
	TH1F *histPosE = new TH1F("histPosE", "PosE hist", xNChannels, xMin, xMax);
	TH1F *histNegPx = new TH1F("histNegPx", "NegPx hist", xNChannels, xMin, xMax);
	TH1F *histNegPy = new TH1F("histNegPy", "NegPy hist", xNChannels, xMin, xMax);
	TH1F *histNegPz = new TH1F("histNegPz", "NegPz hist", xNChannels, xMin, xMax);
	TH1F *histNegE = new TH1F("histNegE", "NegE hist", xNChannels, xMin, xMax);
	
	TH2F *histPosVNegPx = new TH2F("histPosVNegPx", "PosVNegPx hist", xNChannels, xMin, xMax, yNChannels, yMin, yMax);
	TH2F *histPosVNegPy = new TH2F("histPosVNegPy", "PosVNegPy hist", xNChannels, xMin, xMax, yNChannels, yMin, yMax);
	TH2F *histPosVNegPz = new TH2F("histPosVNegPz", "PosVNegPz hist", xNChannels, xMin, xMax, yNChannels, yMin, yMax);
	TH2F *histPosVNegE = new TH2F("histPosVNegE", "PosVNegE hist", xNChannels, xMin, xMax, yNChannels, yMin, yMax);
	
	//TH1F *histRapidityPosPrime = new TH1F("histRapidityPosPrime", "RapidityPosPrime hist", xNChannels, xMin, xMax);
	//TH1F *histRapidityNegPrime = new TH1F("histRapidityNegPrime", "RapidityNegPrime hist", xNChannels, xMin, xMax);
	//TH1F *histRapidityPos = new TH1F("histRapidityPos", "RapidityPos hist", xNChannels, xMin, xMax);
	//TH1F *histRapidityNeg = new TH1F("histRapidityNeg", "RapidityNeg  hist", xNChannels, xMin, xMax);
	
	TLorentzVector pos_muon_prime;
	TLorentzVector neg_muon_prime;
	TLorentzVector tot_muon_prime;
	TLorentzVector pos_muon;
	TLorentzVector neg_muon;
	TLorentzVector tot_muon;
	
	TVector3 tot_muon_boost;
	
	TLorentzRotation l;
	TLorentzRotation l2;
	
	/////////////////////////////////////////////////////// Set up vectors to iterate the creation of all canvases //////////////////////////////////////////////////////////////////////////
	
	///// first 12 1D histograms /////
	
	const int numOf1DCanvasesSet1 = 12;
	
	std::vector<TCanvas*> canvases(numOf1DCanvasesSet1);

	std::vector<std::string> titles = {"pTLep Histogram", "etaLep Histogram", "phiLep Histogram", "pTLepPos Histogram", 
	"etaLepPos Histogram", "phiLepPos Histogram",  "pTLepNeg Histogram",  "etaLepNeg Histogram",  "phiLepNeg Histogram", "nLep histogram",
	"Invariant Mass Histogram", "Invariant Mass Histogram 2"};
	
	/* "PosPxPrime hist", "PosPyPrime hist", "PosPzPrime hist", "PosEPrime hist",
	"NegPxPrime hist", "NegPyPrime hist", "NegPzPrime hist", "NegEPrime hist","TotPxPrime hist", "TotPyPrime hist", "TotPzPrime hist", "TotEPrime hist",
	"TotPx hist", "TotPy hist", "TotPz hist", "TotE hist", "PosPx hist", "PosPy hist", "PosPz hist", "PosE hist","NegPx hist", "NegPy hist", "PosPz hist", "PosE hist"};*/

	std::vector<TH1F*> histos = {histpTLep, histetaLep, histphiLep, histpTLepPos, histetaLepPos, histphiLepPos,
	histpTLepNeg, histetaLepNeg, histphiLepNeg, histnLep, histM, histM2};
	
	///// PosVNeg eta, phi, pT 2D histograms /////
	
	const int numOf2DCanvasesSet1 = 3;
	
	std::vector<TCanvas*> canvases2DSet1(numOf2DCanvasesSet1);
	
	std::vector<std::string> titles2DSet1 = {"PosVNegpTLep Histogram", "PosVNegetaLep Histogram", "PosVNegphiLep Histogram"};
	
	std::vector<TH2F*> histos2DSet1 = {histPosVNegpTLep, histPosVNegetaLep, histPosVNegphiLep};
	
	///// 1D division histograms /////
	
	const int numOf1DCanvasesDiv = 6;
	const int numOfDivisions = 4;
	
	std::vector<TCanvas*> canvases1DDiv(numOf1DCanvasesDiv);
	
	std::vector<std::string> titles1DDiv = {"Pos Momentum Histogram", "Neg Momentum Histogram", 
	"Tot Momentum Histogram", "Pos Prime Momentum Histogram", "Neg Prime Momentum Histogram", 
	"Tot Prime Momentum Histogram"};
				std::vector<std::string> PosTitles = {"PosPx Histogram", "PosPy Histogram", "PosPz Histogram",
				"PosE Histogram"};
				std::vector<TH1F*> PosHistos = {histPosPx, histPosPy, histPosPz, histPosE};
				
				std::vector<std::string> NegTitles = {"NegPx Histogram", "NegPy Histogram", "NegPz Histogram",
				"NegE Histogram"};
				std::vector<TH1F*> NegHistos = {histNegPx, histNegPy, histNegPz, histNegE};
				
				std::vector<std::string> TotTitles = {"TotPx Histogram", "TotPy Histogram", "TotPz Histogram",
				"TotE Histogram"};
				std::vector<TH1F*> TotHistos = {histTotPx, histTotPy, histTotPz, histTotE};
				
				std::vector<std::string> PosPrimeTitles = {"PosPxPrime Histogram", "PosPyPrime Histogram", "PosPzPrime Histogram",
				"PosEPrime Histogram"};
				std::vector<TH1F*> PosPrimeHistos = {histPosPxPrime, histPosPyPrime, histPosPzPrime, histPosEPrime};
				
				std::vector<std::string> NegPrimeTitles = {"NegPxPrime Histogram", "NegPyPrime Histogram", "NegPzPrime Histogram",
				"NegEPrime Histogram"};
				std::vector<TH1F*> NegPrimeHistos = {histNegPxPrime, histNegPyPrime, histNegPzPrime, histNegEPrime};
				
				std::vector<std::string> TotPrimeTitles = {"TotPxPrime Histogram", "TotPyPrime Histogram", "TotPzPrime Histogram",
				"TotEPrime Histogram"};
				std::vector<TH1F*> TotPrimeHistos = {histTotPxPrime, histTotPyPrime, histTotPzPrime, histTotEPrime};
			
				std::vector<std::vector<TH1F*>> DivHistos = {PosHistos, NegHistos, TotHistos, PosPrimeHistos, NegPrimeHistos, TotPrimeHistos};
	
	//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	
	//Creating path to root file and getting tree//
	strcpy(rootFilePath1, rootFilePrefix);
	strcat(rootFilePath1, rootFileName);
	cout<<"Root file path 1: "<<rootFilePath1<<endl;
	
	TFile *file = new TFile(rootFilePath1);
	TTree *t1 = (TTree*)file->Get(TreeName1);
	
	t1->SetBranchAddress("pTLep",&pTLep[0]);
	t1->SetBranchAddress("etaLep",&etaLep[0]);
	t1->SetBranchAddress("phiLep",&phiLep[0]);
	t1->SetBranchAddress("flavLep",&flavLep[0]);
	
	t1->SetBranchAddress("nLeps35",&nLeps35);
	
	int NEntries1 = t1->GetEntries();
	cout<<"NEntries 1: "<<NEntries1<<endl;
	
	
	//NEntries1 = 20000; //comment out if the whole data set is wanted
	
	//Loop to get entries and fill hists (ALL LEPS VERSON)//
	if(!IncludeOnlyMuons){
		for(int i = 0; i < NEntries1; i++){
			TLorentzVector pos_muon_prime;
			TLorentzVector neg_muon_prime;
			TLorentzVector tot_muon_prime;
			TLorentzVector pos_muon;
			TLorentzVector neg_muon;
			TLorentzVector tot_muon;
			t1->GetEntry(i);
			if(nLeps35 >= 2){
				histnLep->Fill(nLeps35);
				cout<<"nLeps35: "<<nLeps35<<endl;
				double pTLepPos, etaLepPos, phiLepPos;
				double pos_Px_prime, pos_Py_prime, pos_Pz_prime, pos_E_prime;
				double pTLepNeg, etaLepNeg, phiLepNeg;
				double neg_Px_prime, neg_Py_prime, neg_Pz_prime, neg_E_prime;
				double tot_Px_prime, tot_Py_prime, tot_Pz_prime, tot_E_prime;
				double tot_Px, tot_Py, tot_Pz, tot_E;
				double neg_Px, neg_Py, neg_Pz, neg_E;
				double pos_Px, pos_Py, pos_Pz, pos_E;
				for(int j = 0; j < 2; j++){
					cout<<" i: "<<i<<" j: "<<j<<" flavLep: "<<flavLep[j]<<endl;
					if(flavLep[j] > 0){
						pos_muon_prime.SetPtEtaPhiM(pTLep[j],etaLep[j],phiLep[j],0.105);
						histpTLep->Fill(pTLep[j]);
						histetaLep->Fill(etaLep[j]);
						histphiLep->Fill(phiLep[j]);
						
						histpTLepPos->Fill(pTLep[j]);
						histetaLepPos->Fill(etaLep[j]);
						histphiLepPos->Fill(phiLep[j]);
						
						pTLepPos = pTLep[j];
						etaLepPos = etaLep[j];
						phiLepPos = phiLep[j];
						
						pos_Px_prime = pos_muon_prime.Px();
						pos_Py_prime = pos_muon_prime.Py();
						pos_Pz_prime = pos_muon_prime.Pz();
						pos_E_prime = pos_muon_prime.E();
						histPosPxPrime->Fill(pos_Px_prime);
						histPosPyPrime->Fill(pos_Py_prime);
						histPosPzPrime->Fill(pos_Pz_prime);
						histPosEPrime->Fill(pos_E_prime);
					}
					if(flavLep[j] < 0){
						neg_muon_prime.SetPtEtaPhiM(pTLep[j],etaLep[j],phiLep[j],0.105);
						histpTLep->Fill(pTLep[j]);
						histetaLep->Fill(etaLep[j]);
						histphiLep->Fill(phiLep[j]);
						
						histpTLepNeg->Fill(pTLep[j]);
						histetaLepNeg->Fill(etaLep[j]);
						histphiLepNeg->Fill(phiLep[j]);
						
						pTLepNeg = pTLep[j];
						etaLepNeg = etaLep[j];
						phiLepNeg = phiLep[j];
						
						
						neg_Px_prime = neg_muon_prime.Px();
						neg_Py_prime = neg_muon_prime.Py();
						neg_Pz_prime = neg_muon_prime.Pz();
						neg_E_prime = neg_muon_prime.E();
						histNegPxPrime->Fill(neg_Px_prime);
						histNegPyPrime->Fill(neg_Py_prime);
						histNegPzPrime->Fill(neg_Pz_prime);
						histNegEPrime->Fill(neg_E_prime);
						
					}
					histPosVNegpTLep->Fill(pTLepPos, pTLepNeg);
					histPosVNegetaLep->Fill(etaLepPos, etaLepNeg);
					histPosVNegphiLep->Fill(phiLepPos, phiLepNeg);
					
					histPosVNegPxPrime->Fill(pos_Px_prime, neg_Px_prime);
					histPosVNegPyPrime->Fill(pos_Py_prime, neg_Py_prime);
					histPosVNegPzPrime->Fill(pos_Pz_prime, neg_Pz_prime);
					histPosVNegEPrime->Fill(pos_E_prime, neg_E_prime);
				}
				if((flavLep[0] > 0 && flavLep[1] < 0) || (flavLep[0] < 0 && flavLep[1] > 0)){
					tot_muon_prime = pos_muon_prime + neg_muon_prime; 
					InvarientMass = tot_muon_prime.M();
					histM->Fill(InvarientMass);
					
					TVector3 tot_muon_boost;
					cout<<"tot_muon_boost vect:"<<endl;
					tot_muon_boost = tot_muon_prime.BoostVector();
					tot_muon_boost.Print();
					
					tot_Px_prime = tot_muon_prime.Px();
					tot_Py_prime = tot_muon_prime.Py();
					tot_Pz_prime = tot_muon_prime.Pz();
					tot_E_prime = tot_muon_prime.E();
					histTotPxPrime->Fill(tot_Px_prime);
					histTotPyPrime->Fill(tot_Py_prime);
					histTotPzPrime->Fill(tot_Pz_prime);
					histTotEPrime->Fill(tot_E_prime);
					
					double Z0mass = 91.19;
					double beta = tot_muon_prime.Beta();
					double gamma = tot_muon_prime.Gamma();
					double tot_P_prime = tot_muon_prime.P();
					double V = beta*c;
					double Vx = tot_Px_prime/(InvarientMass*gamma);
					double Vy = tot_Py_prime/(InvarientMass*gamma);
					double Vz = tot_Pz_prime/(InvarientMass*gamma);
					
					TLorentzRotation l;
					TLorentzRotation l2;
					l.Boost(0,0,0);
					l.Boost(tot_muon_boost);
					l2 = l.Inverse();
					neg_muon = l2*neg_muon_prime;
					pos_muon = l2*pos_muon_prime;
					tot_muon = pos_muon + neg_muon; //l2*tot_muon_prime;
					
					tot_Px = tot_muon.Px();
					tot_Py = tot_muon.Py();
					tot_Pz = tot_muon.Pz();
					tot_E = tot_muon.E();
					
					// tot_E = (tot_E_prime/gamma) + V*tot_P_prime;
					// tot_Px = (tot_Px_prime/gamma) + ((Vx*tot_E)/(c*c));
					// tot_Py = (tot_Py_prime/gamma) + ((Vy*tot_E)/(c*c));
					// tot_Pz = (tot_Pz_prime/gamma) + ((Vz*tot_E)/(c*c));
					histTotPx->Fill(tot_Px);
					histTotPy->Fill(tot_Py);
					histTotPz->Fill(tot_Pz);
					histTotE->Fill(tot_E);
				
					neg_Px = neg_muon.Px();
					neg_Py = neg_muon.Py();
					neg_Pz = neg_muon.Pz();
					neg_E = neg_muon.E();
					
					pos_Px = pos_muon.Px();
					pos_Py = pos_muon.Py();
					pos_Pz = pos_muon.Pz();
					pos_E = pos_muon.E();
					
					// This no work in 3D //
					// double neg_P_prime = neg_muon_prime.P();
					// double neg_Px, neg_Py, neg_Pz, neg_E;
					// neg_E = (neg_E_prime/gamma) + V*neg_P_prime;
					// neg_Px = (neg_Px_prime/gamma) + ((Vx*neg_E)/(c*c));
					// neg_Py = (neg_Py_prime/gamma) + ((Vy*neg_E)/(c*c));
					// neg_Pz = (neg_Pz_prime/gamma) + ((Vz*neg_E)/(c*c));
					histNegPx->Fill(neg_Px);
					histNegPy->Fill(neg_Py);
					histNegPz->Fill(neg_Pz);
					histNegE->Fill(neg_E);
					
					// double pos_P_prime = pos_muon_prime.P();
					// double pos_Px, pos_Py, pos_Pz, pos_E;
					// pos_E = (pos_E_prime/gamma) + V*pos_P_prime;
					// pos_Px = (pos_Px_prime/gamma) + ((Vx*pos_E)/(c*c));
					// pos_Py = (pos_Py_prime/gamma) + ((Vy*pos_E)/(c*c));
					// pos_Pz = (pos_Pz_prime/gamma) + ((Vz*pos_E)/(c*c));
					histPosPx->Fill(pos_Px);
					histPosPy->Fill(pos_Py);
					histPosPz->Fill(pos_Pz);
					histPosE->Fill(pos_E);
					
					histPosVNegPx->Fill(pos_Px, neg_Px);
					histPosVNegPy->Fill(pos_Py, neg_Py);
					histPosVNegPz->Fill(pos_Pz, neg_Pz);
					histPosVNegE->Fill(pos_E, neg_E);
					
					pos_muon.SetPxPyPzE(pos_Px, pos_Py, pos_Pz, pos_E);
					neg_muon.SetPxPyPzE(neg_Px, neg_Py, neg_Pz, neg_E);
					tot_muon = pos_muon + neg_muon;
					
					double InvarientMass_2 = tot_muon.M();
					histM2->Fill(InvarientMass_2);
					
					cout<<"Vx: "<<Vx<<" Vy: "<<Vy<<" Vz: "<<Vz<<endl;
					// cout<<"Vactual: "<<V<<" Vsupposed: "<<sqrt(Vx*Vx + Vy*Vy + Vz*Vz)<<endl;
					// cout<<"gammaActual: "<<gamma<<" gammaSupposed: "<<(1/(sqrt(1-beta*beta)))<<endl;
					
					cout<<"pos_Px_prime: "<<pos_Px_prime<<" pos_Py_prime: "<<pos_Py_prime<<" pos_Pz_prime: "<<pos_Pz_prime<<endl;
					cout<<"neg_Px_prime: "<<neg_Px_prime<<" neg_Py_prime: "<<neg_Py_prime<<" neg_Pz_prime: "<<neg_Pz_prime<<endl;
					cout<<"pos_Px: "<<pos_Px<<" pos_Py: "<<pos_Py<<" pos_Pz: "<<pos_Pz<<endl;
					cout<<"neg_Px: "<<neg_Px<<" neg_Py: "<<neg_Py<<" neg_Pz: "<<neg_Pz<<endl;
					
					NumOfPosNegPairs++;
				}
				// if(flavLep[0] < 0 && flavLep[1] > 0){
					// tot_muon = pos_muon + neg_muon; 
					// InvarientMass = tot_muon.M();
					// histM->Fill(InvarientMass);
					
					// double tot_Px_prime, tot_Py_prime, tot_Pz_prime, tot_E_prime;
					// tot_Px_prime = tot_muon.Px();
					// tot_Py_prime = tot_muon.Py();
					// tot_Pz_prime = tot_muon.Pz();
					// tot_E_prime = tot_muon.E();
					// histTotPxPrime->Fill(tot_Px_prime);
					// histTotPyPrime->Fill(tot_Py_prime);
					// histTotPzPrime->Fill(tot_Pz_prime);
					// histTotEPrime->Fill(tot_E_prime);
					
					// double Z0mass = 91.19;
					// double beta = tot_muon.Beta();
					// double gamma = tot_muon.Gamma();
					// double P_prime = tot_muon.P();
					// double V = beta*c;
					// double Vx = tot_Px_prime/(InvarientMass*gamma);
					// double Vy = tot_Py_prime/(InvarientMass*gamma);
					// double Vz = tot_Pz_prime/(InvarientMass*gamma);
					
					// double tot_Px, tot_Py, tot_Pz, tot_E;
					// tot_E = (tot_E_prime/gamma) + V*P_prime;
					// tot_Px = (tot_Px_prime/gamma) + ((Vx*tot_E)/(c*c));
					// tot_Py = (tot_Py_prime/gamma) + ((Vy*tot_E)/(c*c));
					// tot_Pz = (tot_Pz_prime/gamma) + ((Vz*tot_E)/(c*c));
					// histTotPx->Fill(tot_Px);
					// histTotPy->Fill(tot_Py);
					// histTotPz->Fill(tot_Pz);
					// histTotE->Fill(tot_E);
					// cout<<"tot_Px_prime: "<<tot_Px_prime<<" tot_Py_prime: "<<tot_Py_prime<<" tot_Pz_prime: "<<tot_Pz_prime<<endl;
					
					// NumOfPosNegPairs++;
				// }
			}
		}
	}
	
	//ONLY MUONS//
	if(IncludeOnlyMuons){
		for(int i = 0; i < NEntries1; i++){
			TLorentzVector pos_muon_prime;
			TLorentzVector neg_muon_prime;
			TLorentzVector tot_muon_prime;
			TLorentzVector pos_muon;
			TLorentzVector neg_muon;
			TLorentzVector tot_muon;
			t1->GetEntry(i);
			if(nLeps35 >= 2){
				histnLep->Fill(nLeps35);
				cout<<"nLeps35: "<<nLeps35<<endl;
				double pTLepPos, etaLepPos, phiLepPos;
				double pos_Px_prime, pos_Py_prime, pos_Pz_prime, pos_E_prime;
				double pTLepNeg, etaLepNeg, phiLepNeg;
				double neg_Px_prime, neg_Py_prime, neg_Pz_prime, neg_E_prime;
				double tot_Px_prime, tot_Py_prime, tot_Pz_prime, tot_E_prime;
				double tot_Px, tot_Py, tot_Pz, tot_E;
				double neg_Px, neg_Py, neg_Pz, neg_E;
				double pos_Px, pos_Py, pos_Pz, pos_E;
				for(int j = 0; j < 2; j++){
					cout<<" i: "<<i<<" j: "<<j<<" flavLep: "<<flavLep[j]<<endl;
					if(flavLep[j] == 2){
						pos_muon_prime.SetPtEtaPhiM(pTLep[j],etaLep[j],phiLep[j],0.105);
						histpTLep->Fill(pTLep[j]);
						histetaLep->Fill(etaLep[j]);
						histphiLep->Fill(phiLep[j]);
						
						histpTLepPos->Fill(pTLep[j]);
						histetaLepPos->Fill(etaLep[j]);
						histphiLepPos->Fill(phiLep[j]);
						
						//double pTLepPos, etaLepPos, phiLepPos;
						pTLepPos = pTLep[j];
						etaLepPos = etaLep[j];
						phiLepPos = phiLep[j];
						
						//double pos_Px_prime, pos_Py_prime, pos_Pz_prime, pos_E_prime;
						pos_Px_prime = pos_muon_prime.Px();
						pos_Py_prime = pos_muon_prime.Py();
						pos_Pz_prime = pos_muon_prime.Pz();
						pos_E_prime = pos_muon_prime.E();
						histPosPxPrime->Fill(pos_Px_prime);
						histPosPyPrime->Fill(pos_Py_prime);
						histPosPzPrime->Fill(pos_Pz_prime);
						histPosEPrime->Fill(pos_E_prime);
						
						//RapidityPosPrime = pos_muon_prime.Rapidity();
						//histRapidityPosPrime->Fill(RapidityPosPrime);
					}
					if(flavLep[j] == -2){
						neg_muon_prime.SetPtEtaPhiM(pTLep[j],etaLep[j],phiLep[j],0.105);
						histpTLep->Fill(pTLep[j]);
						histetaLep->Fill(etaLep[j]);
						histphiLep->Fill(phiLep[j]);
						
						histpTLepNeg->Fill(pTLep[j]);
						histetaLepNeg->Fill(etaLep[j]);
						histphiLepNeg->Fill(phiLep[j]);
						
						//double pTLepNeg, etaLepNeg, phiLepNeg;
						pTLepNeg = pTLep[j];
						etaLepNeg = etaLep[j];
						phiLepNeg = phiLep[j];
						
						//double neg_Px_prime, neg_Py_prime, neg_Pz_prime, neg_E_prime;
						neg_Px_prime = neg_muon_prime.Px();
						neg_Py_prime = neg_muon_prime.Py();
						neg_Pz_prime = neg_muon_prime.Pz();
						neg_E_prime = neg_muon_prime.E();
						histNegPxPrime->Fill(neg_Px_prime);
						histNegPyPrime->Fill(neg_Py_prime);
						histNegPzPrime->Fill(neg_Pz_prime);
						histNegEPrime->Fill(neg_E_prime);
						
						//RapidityNegPrime = neg_muon_prime.Rapidity();
						//histRapidityNegPrime->Fill(RapidityNegPrime);
						
					}
					histPosVNegpTLep->Fill(pTLepPos, pTLepNeg);
					histPosVNegetaLep->Fill(etaLepPos, etaLepNeg);
					histPosVNegphiLep->Fill(phiLepPos, phiLepNeg);
					
					histPosVNegPxPrime->Fill(pos_Px_prime, neg_Px_prime);
					histPosVNegPyPrime->Fill(pos_Py_prime, neg_Py_prime);
					histPosVNegPzPrime->Fill(pos_Pz_prime, neg_Pz_prime);
					histPosVNegEPrime->Fill(pos_E_prime, neg_E_prime);
				}
				if((flavLep[0] == 2 && flavLep[1] == -2) || (flavLep[0] == -2 && flavLep[1] == 2)){
					tot_muon_prime = pos_muon_prime + neg_muon_prime; 
					InvarientMass = tot_muon_prime.M();
					histM->Fill(InvarientMass);
					
					TVector3 tot_muon_boost;
					cout<<"tot_muon_boost vect:"<<endl;
					tot_muon_boost = tot_muon_prime.BoostVector();
					tot_muon_boost.Print();
					
					//double tot_Px_prime, tot_Py_prime, tot_Pz_prime, tot_E_prime;
					tot_Px_prime = tot_muon_prime.Px();
					tot_Py_prime = tot_muon_prime.Py();
					tot_Pz_prime = tot_muon_prime.Pz();
					tot_E_prime = tot_muon_prime.E();
					histTotPxPrime->Fill(tot_Px_prime);
					histTotPyPrime->Fill(tot_Py_prime);
					histTotPzPrime->Fill(tot_Pz_prime);
					histTotEPrime->Fill(tot_E_prime);
					
					double Z0mass = 91.19;
					double beta = tot_muon_prime.Beta();
					double gamma = tot_muon_prime.Gamma();
					double tot_P_prime = tot_muon_prime.P();
					double V = beta*c;
					double Vx = tot_Px_prime/(InvarientMass*gamma);
					double Vy = tot_Py_prime/(InvarientMass*gamma);
					double Vz = tot_Pz_prime/(InvarientMass*gamma);
					
					TLorentzRotation l;
					TLorentzRotation l2;
					l.Boost(0,0,0);
					l.Boost(tot_muon_boost);
					l2 = l.Inverse();
					neg_muon = l2*neg_muon_prime;
					pos_muon = l2*pos_muon_prime;
					tot_muon = pos_muon + neg_muon; //l2*tot_muon_prime;
					
					//double tot_Px, tot_Py, tot_Pz, tot_E;
					tot_Px = tot_muon.Px();
					tot_Py = tot_muon.Py();
					tot_Pz = tot_muon.Pz();
					tot_E = tot_muon.E();
					
					// tot_E = (tot_E_prime/gamma) + V*tot_P_prime;
					// tot_Px = (tot_Px_prime/gamma) + ((Vx*tot_E)/(c*c));
					// tot_Py = (tot_Py_prime/gamma) + ((Vy*tot_E)/(c*c));
					// tot_Pz = (tot_Pz_prime/gamma) + ((Vz*tot_E)/(c*c));
					histTotPx->Fill(tot_Px);
					histTotPy->Fill(tot_Py);
					histTotPz->Fill(tot_Pz);
					histTotE->Fill(tot_E);
					
					//double neg_Px, neg_Py, neg_Pz, neg_E;
					neg_Px = neg_muon.Px();
					neg_Py = neg_muon.Py();
					neg_Pz = neg_muon.Pz();
					neg_E = neg_muon.E();
					
					//double pos_Px, pos_Py, pos_Pz, pos_E;
					pos_Px = pos_muon.Px();
					pos_Py = pos_muon.Py();
					pos_Pz = pos_muon.Pz();
					pos_E = pos_muon.E();
					
					// This no work in 3D //
					// double neg_P_prime = neg_muon_prime.P();
					// double neg_Px, neg_Py, neg_Pz, neg_E;
					// neg_E = (neg_E_prime/gamma) + V*neg_P_prime;
					// neg_Px = (neg_Px_prime/gamma) + ((Vx*neg_E)/(c*c));
					// neg_Py = (neg_Py_prime/gamma) + ((Vy*neg_E)/(c*c));
					// neg_Pz = (neg_Pz_prime/gamma) + ((Vz*neg_E)/(c*c));
					histNegPx->Fill(neg_Px);
					histNegPy->Fill(neg_Py);
					histNegPz->Fill(neg_Pz);
					histNegE->Fill(neg_E);
					
					// double pos_P_prime = pos_muon_prime.P();
					// double pos_Px, pos_Py, pos_Pz, pos_E;
					// pos_E = (pos_E_prime/gamma) + V*pos_P_prime;
					// pos_Px = (pos_Px_prime/gamma) + ((Vx*pos_E)/(c*c));
					// pos_Py = (pos_Py_prime/gamma) + ((Vy*pos_E)/(c*c));
					// pos_Pz = (pos_Pz_prime/gamma) + ((Vz*pos_E)/(c*c));
					histPosPx->Fill(pos_Px);
					histPosPy->Fill(pos_Py);
					histPosPz->Fill(pos_Pz);
					histPosE->Fill(pos_E);
					
					histPosVNegPx->Fill(pos_Px, neg_Px);
					histPosVNegPy->Fill(pos_Py, neg_Py);
					histPosVNegPz->Fill(pos_Pz, neg_Pz);
					histPosVNegE->Fill(pos_E, neg_E);
					
					pos_muon.SetPxPyPzE(pos_Px, pos_Py, pos_Pz, pos_E);
					neg_muon.SetPxPyPzE(neg_Px, neg_Py, neg_Pz, neg_E);
					tot_muon = pos_muon + neg_muon;
					
					double InvarientMass_2 = tot_muon.M();
					histM2->Fill(InvarientMass_2);
					
					cout<<"Vx: "<<Vx<<" Vy: "<<Vy<<" Vz: "<<Vz<<endl;
					// cout<<"Vactual: "<<V<<" Vsupposed: "<<sqrt(Vx*Vx + Vy*Vy + Vz*Vz)<<endl;
					// cout<<"gammaActual: "<<gamma<<" gammaSupposed: "<<(1/(sqrt(1-beta*beta)))<<endl;
					
					cout<<"pos_Px_prime: "<<pos_Px_prime<<" pos_Py_prime: "<<pos_Py_prime<<" pos_Pz_prime: "<<pos_Pz_prime<<endl;
					cout<<"neg_Px_prime: "<<neg_Px_prime<<" neg_Py_prime: "<<neg_Py_prime<<" neg_Pz_prime: "<<neg_Pz_prime<<endl;
					cout<<"pos_Px: "<<pos_Px<<" pos_Py: "<<pos_Py<<" pos_Pz: "<<pos_Pz<<endl;
					cout<<"neg_Px: "<<neg_Px<<" neg_Py: "<<neg_Py<<" neg_Pz: "<<neg_Pz<<endl;
					
					NumOfPosNegPairs++;
				}
				// if(flavLep[0] < 0 && flavLep[1] > 0){
					// tot_muon = pos_muon + neg_muon; 
					// InvarientMass = tot_muon.M();
					// histM->Fill(InvarientMass);
					
					// double tot_Px_prime, tot_Py_prime, tot_Pz_prime, tot_E_prime;
					// tot_Px_prime = tot_muon.Px();
					// tot_Py_prime = tot_muon.Py();
					// tot_Pz_prime = tot_muon.Pz();
					// tot_E_prime = tot_muon.E();
					// histTotPxPrime->Fill(tot_Px_prime);
					// histTotPyPrime->Fill(tot_Py_prime);
					// histTotPzPrime->Fill(tot_Pz_prime);
					// histTotEPrime->Fill(tot_E_prime);
					
					// double Z0mass = 91.19;
					// double beta = tot_muon.Beta();
					// double gamma = tot_muon.Gamma();
					// double P_prime = tot_muon.P();
					// double V = beta*c;
					// double Vx = tot_Px_prime/(InvarientMass*gamma);
					// double Vy = tot_Py_prime/(InvarientMass*gamma);
					// double Vz = tot_Pz_prime/(InvarientMass*gamma);
					
					// double tot_Px, tot_Py, tot_Pz, tot_E;
					// tot_E = (tot_E_prime/gamma) + V*P_prime;
					// tot_Px = (tot_Px_prime/gamma) + ((Vx*tot_E)/(c*c));
					// tot_Py = (tot_Py_prime/gamma) + ((Vy*tot_E)/(c*c));
					// tot_Pz = (tot_Pz_prime/gamma) + ((Vz*tot_E)/(c*c));
					// histTotPx->Fill(tot_Px);
					// histTotPy->Fill(tot_Py);
					// histTotPz->Fill(tot_Pz);
					// histTotE->Fill(tot_E);
					// cout<<"tot_Px_prime: "<<tot_Px_prime<<" tot_Py_prime: "<<tot_Py_prime<<" tot_Pz_prime: "<<tot_Pz_prime<<endl;
					
					// NumOfPosNegPairs++;
				// }
			}
		}
	}
	
	cout<<"NumOfPosNegPairs: "<<NumOfPosNegPairs<<endl; 
	
	//////Graphics//////
	cout<<"Making histrograms..."<<endl;
	
	/////////////// Loop to generate first eleven 1D canvases
	for (int i = 0; i < numOf1DCanvasesSet1; ++i) {
	   TCanvas *c = new TCanvas(("c" + std::to_string(i)).c_str(), titles[i].c_str(), 200, 10, 900, 700);
	   canvases[i] = c;
			c->SetFillColor(10);
			c->SetGrid();
			c->GetFrame()->SetFillColor(10);
			c->GetFrame()->SetBorderSize(12);
			c->Range(0,0,1,1);

			gStyle->SetOptStat(1);

			c->cd();
			histos[i]->SetXTitle(titles[i].c_str());
			histos[i]->SetYTitle("Counts");
			histos[i]->SetStats(1);
			histos[i]->Draw();
			if(SaveAsRootFile){
				TFile outfile1(rootOutputFileName1, "RECREATE");
				c->Write(rootOutputFileName1);
				outfile1.Close();
		}
		//cout << "Creating canvas: " << ("c" + std::to_string(i)).c_str() << " " << titles[i].c_str() << endl;
	}
	
		/////////////// Loop to generate first set  of 2D canvases
	for (int i = 0; i < numOf2DCanvasesSet1; ++i) {
	   TCanvas *c = new TCanvas(("c" + std::to_string(i+numOf1DCanvasesSet1)).c_str(), titles2DSet1[i].c_str(), 200, 10, 900, 700);
	   canvases2DSet1[i] = c;
			c->SetFillColor(10);
			c->SetGrid();
			c->GetFrame()->SetFillColor(10);
			c->GetFrame()->SetBorderSize(12);
			c->Range(0,0,1,1);

			gStyle->SetOptStat(1);

			c->cd();
			histos2DSet1[i]->SetXTitle(titles2DSet1[i].c_str());
			histos2DSet1[i]->SetYTitle("Counts");
			histos2DSet1[i]->SetStats(1);
			histos2DSet1[i]->Draw();
			if(SaveAsRootFile){
				TFile outfile1(rootOutputFileName1, "RECREATE");
				c->Write(rootOutputFileName1);
				outfile1.Close();
		}
		//cout << "Creating canvas: " << ("c" + std::to_string(i+numOf1DCanvasesSet1)).c_str() << " " << titles2DSet1[i].c_str() << endl;
	}
	
	/////////////// Loop to generate first set 1D canvases with divisions
	for (int i = 0; i < numOf1DCanvasesDiv; ++i) {
	   TCanvas *c = new TCanvas(("c" + std::to_string(i+numOf1DCanvasesSet1+numOf2DCanvasesSet1)).c_str(), titles1DDiv[i].c_str(), 200, 10, 900, 700);  
	   canvases1DDiv[i] = c;
			c->SetFillColor(10);
			c->SetGrid();
			c->GetFrame()->SetFillColor(10);
			c->GetFrame()->SetBorderSize(12);
			c->Range(0,0,1,1);

			gStyle->SetOptStat(1);
			
			for (int j = 0; j < numOfDivisions; ++j){
				c->cd(j+1);
				DivHistos[i][j]->SetXTitle("test");
				DivHistos[i][j]->SetYTitle("Counts");
				DivHistos[i][j]->SetStats(1);
				DivHistos[i][j]->Draw();
				cout << "\t Creating histrogram: " << typeid(DivHistos[i][j]).name() << endl;
			}
			
		cout << "Creating canvas: " << ("c" + std::to_string(i+numOf1DCanvasesSet1+numOf2DCanvasesSet1)).c_str() << " " << titles1DDiv[i].c_str() << endl;
	}

/*
	TCanvas *c15 = new TCanvas("c15","Pos Momentum Prime Histrogram" ,200,10,900,700);
		c15->Divide(2,2);
		c15->SetFillColor(10);
		c15->SetGrid();
		c15->GetFrame()->SetFillColor(10);
		c15->GetFrame()->SetBorderSize(12);
		c15->Range(0,0,1,1);

		gStyle->SetOptStat(1);
			
			c15->cd(1);
			histPosPxPrime->SetXTitle("PosPx");
			histPosPxPrime->SetYTitle("counts");
			histPosPxPrime->SetStats(1);
			histPosPxPrime->Draw();
			
			c15->cd(2);
			histPosPyPrime->SetXTitle("PosPy");
			histPosPyPrime->SetYTitle("counts");
			histPosPyPrime->SetStats(1);
			histPosPyPrime->Draw();
		
			c15->cd(3);
			histPosPzPrime->SetXTitle("PosPz");
			histPosPzPrime->SetYTitle("counts");
			histPosPzPrime->SetStats(1);
			histPosPzPrime->Draw();
			
			c15->cd(4);
			histPosEPrime->SetXTitle("PosE");
			histPosEPrime->SetYTitle("counts");
			histPosEPrime->SetStats(1);
			histPosEPrime->Draw();
			
			if(SaveAsRootFile){
				TFile outfile15(rootOutputFileName15, "RECREATE");
				c15->Write(rootOutputFileName15);
				outfile15.Close();
			}


	TCanvas *c1 = new TCanvas("c1","pTLep Histrogram" ,200,10,900,700);
		c1->SetFillColor(10);
		c1->SetGrid();
		c1->GetFrame()->SetFillColor(10);
		c1->GetFrame()->SetBorderSize(12);
		c1->Range(0,0,1,1);

		gStyle->SetOptStat(1);
			
			c1->cd();
			histpTLep->SetXTitle("pTLep");
			histpTLep->SetYTitle("Counts");
			histpTLep->SetStats(1);
			histpTLep->Draw();
			if(SaveAsRootFile){
				TFile outfile1(rootOutputFileName1, "RECREATE");
				c1->Write(rootOutputFileName1);
				outfile1.Close();
			}
			
			
	TCanvas *c2 = new TCanvas("c2","etaLep Histrogram" ,200,10,900,700);
		c2->SetFillColor(10);
		c2->SetGrid();
		c2->GetFrame()->SetFillColor(10);
		c2->GetFrame()->SetBorderSize(12);
		c2->Range(0,0,1,1);

		gStyle->SetOptStat(1);
			
			c2->cd();
			histetaLep->SetXTitle("etaLep");
			histetaLep->SetYTitle("Counts");
			histetaLep->SetStats(1);
			histetaLep->Draw();
			if(SaveAsRootFile){
				TFile outfile2(rootOutputFileName2, "RECREATE");
				c2->Write(rootOutputFileName2);
				outfile2.Close();
			}


	TCanvas *c3 = new TCanvas("c3","phiLep Histrogram" ,200,10,900,700);
		c3->SetFillColor(10);
		c3->SetGrid();
		c3->GetFrame()->SetFillColor(10);
		c3->GetFrame()->SetBorderSize(12);
		c3->Range(0,0,1,1);

		gStyle->SetOptStat(1);
		
			c3->cd();
			histphiLep->SetXTitle("phiLep");
			histphiLep->SetYTitle("Counts");
			histphiLep->SetStats(1);
			histphiLep->Draw();
			if(SaveAsRootFile){
				TFile outfile3(rootOutputFileName3, "RECREATE");
				c3->Write(rootOutputFileName3);
				outfile3.Close();
			}
			
			
			
	TCanvas *c4 = new TCanvas("c4","pTLepPos Histrogram" ,200,10,900,700);
		c4->SetFillColor(10);
		c4->SetGrid();
		c4->GetFrame()->SetFillColor(10);
		c4->GetFrame()->SetBorderSize(12);
		c4->Range(0,0,1,1);

		gStyle->SetOptStat(1);
			
			c4->cd();
			histpTLepPos->SetXTitle("pTLepPos");
			histpTLepPos->SetYTitle("Counts");
			histpTLepPos->SetStats(1);
			histpTLepPos->Draw();
			if(SaveAsRootFile){
				TFile outfile4(rootOutputFileName4, "RECREATE");
				c4->Write(rootOutputFileName4);
				outfile4.Close();
			}
			
			
	TCanvas *c5 = new TCanvas("c5","etaLepPos Histrogram" ,200,10,900,700);
		c5->SetFillColor(10);
		c5->SetGrid();
		c5->GetFrame()->SetFillColor(10);
		c5->GetFrame()->SetBorderSize(12);
		c5->Range(0,0,1,1);

		gStyle->SetOptStat(1);
			
			c5->cd();
			histetaLepPos->SetXTitle("etaLepPos");
			histetaLepPos->SetYTitle("Counts");
			histetaLepPos->SetStats(1);
			histetaLepPos->Draw();
			if(SaveAsRootFile){
				TFile outfile5(rootOutputFileName5, "RECREATE");
				c5->Write(rootOutputFileName5);
				outfile5.Close();
			}


	TCanvas *c6 = new TCanvas("c6","phiLepPos Histrogram" ,200,10,900,700);
		c6->SetFillColor(10);
		c6->SetGrid();
		c6->GetFrame()->SetFillColor(10);
		c6->GetFrame()->SetBorderSize(12);
		c6->Range(0,0,1,1);

		gStyle->SetOptStat(1);
		
			c6->cd();
			histphiLepPos->SetXTitle("phiLepPos");
			histphiLepPos->SetYTitle("Counts");
			histphiLepPos->SetStats(1);
			histphiLepPos->Draw();
			if(SaveAsRootFile){
				TFile outfile6(rootOutputFileName6, "RECREATE");
				c6->Write(rootOutputFileName6);
				outfile6.Close();
			}
			
	TCanvas *c7 = new TCanvas("c7","pTLepNeg Histrogram" ,200,10,900,700);
		c7->SetFillColor(10);
		c7->SetGrid();
		c7->GetFrame()->SetFillColor(10);
		c7->GetFrame()->SetBorderSize(12);
		c7->Range(0,0,1,1);

		gStyle->SetOptStat(1);
			
			c7->cd();
			histpTLepNeg->SetXTitle("pTLepNeg");
			histpTLepNeg->SetYTitle("Counts");
			histpTLepNeg->SetStats(1);
			histpTLepNeg->Draw();
			if(SaveAsRootFile){
				TFile outfile7(rootOutputFileName7, "RECREATE");
				c7->Write(rootOutputFileName7);
				outfile7.Close();
			}
			
			
	TCanvas *c8 = new TCanvas("c8","etaLepNeg Histrogram" ,200,10,900,700);
		c8->SetFillColor(10);
		c8->SetGrid();
		c8->GetFrame()->SetFillColor(10);
		c8->GetFrame()->SetBorderSize(12);
		c8->Range(0,0,1,1);

		gStyle->SetOptStat(1);
			
			c8->cd();
			histetaLep->SetXTitle("etaLepNeg");
			histetaLep->SetYTitle("Counts");
			histetaLep->SetStats(1);
			histetaLep->Draw();
			if(SaveAsRootFile){
				TFile outfile8(rootOutputFileName8, "RECREATE");
				c8->Write(rootOutputFileName8);
				outfile8.Close();
			}


	TCanvas *c9 = new TCanvas("c9","phiLepNeg Histrogram" ,200,10,900,700);
		c9->SetFillColor(10);
		c9->SetGrid();
		c9->GetFrame()->SetFillColor(10);
		c9->GetFrame()->SetBorderSize(12);
		c9->Range(0,0,1,1);

		gStyle->SetOptStat(1);
		
			c9->cd();
			histphiLepNeg->SetXTitle("phiLepNeg");
			histphiLepNeg->SetYTitle("Counts");
			histphiLepNeg->SetStats(1);
			histphiLepNeg->Draw();
			if(SaveAsRootFile){
				TFile outfile9(rootOutputFileName9, "RECREATE");
				c9->Write(rootOutputFileName9);
				outfile9.Close();
			}
			
			
	TCanvas *c10 = new TCanvas("c10","nLep Histrogram" ,200,10,900,700);
		c10->SetFillColor(10);
		c10->SetGrid();
		c10->GetFrame()->SetFillColor(10);
		c10->GetFrame()->SetBorderSize(12);
		c10->Range(0,0,1,1);

		gStyle->SetOptStat(1);
			
			c10->cd();
			histnLep->SetXTitle("nLepNeg");
			histnLep->SetYTitle("Counts");
			histnLep->SetStats(1);
			histnLep->Draw();
			if(SaveAsRootFile){
				TFile outfile10(rootOutputFileName10, "RECREATE");
				c10->Write(rootOutputFileName10);
				outfile10.Close();
			}
			
	TCanvas *c11 = new TCanvas("c11","Invarient Mass Histrogram" ,200,10,900,700);
		c11->Divide(1,2);
		c11->SetFillColor(10);
		c11->SetGrid();
		c11->GetFrame()->SetFillColor(10);
		c11->GetFrame()->SetBorderSize(12);
		c11->Range(0,0,1,1);

		gStyle->SetOptStat(1);
			
			c11->cd(1);
			histM->SetXTitle("Invarient Mass (GeV)");
			histM->SetYTitle("Counts");
			histM->SetStats(1);
			histM->Draw();
			
			c11->cd(2);
			histM2->SetXTitle("Invarient Mass 2 (GeV)");
			histM2->SetYTitle("Counts");
			histM2->SetStats(1);
			histM2->Draw();
			
			if(SaveAsRootFile){
				TFile outfile11(rootOutputFileName11, "RECREATE");
				c11->Write(rootOutputFileName11);
				outfile11.Close();
			}

	
	TCanvas *c12 = new TCanvas("c12","PosVNegpTLep Histrogram" ,200,10,900,700);
		c12->SetFillColor(10);
		c12->SetGrid();
		c12->GetFrame()->SetFillColor(10);
		c12->GetFrame()->SetBorderSize(12);
		c12->Range(0,0,1,1);

		gStyle->SetOptStat(1);
			
			c12->cd();
			histPosVNegpTLep->SetXTitle("Pos");
			histPosVNegpTLep->SetYTitle("Neg");
			histPosVNegpTLep->SetZTitle("counts");
			histPosVNegpTLep->SetStats(1);
			histPosVNegpTLep->Draw();
			if(SaveAsRootFile){
				TFile outfile12(rootOutputFileName12, "RECREATE");
				c12->Write(rootOutputFileName12);
				outfile12.Close();
			}
	
	TCanvas *c13 = new TCanvas("c13","PosVNegetaLep Histrogram" ,200,10,900,700);
		c13->SetFillColor(10);
		c13->SetGrid();
		c13->GetFrame()->SetFillColor(10);
		c13->GetFrame()->SetBorderSize(12);
		c13->Range(0,0,1,1);

		gStyle->SetOptStat(1);
			
			c13->cd();
			histPosVNegetaLep->SetXTitle("Pos");
			histPosVNegetaLep->SetYTitle("Neg");
			histPosVNegetaLep->SetZTitle("counts");
			histPosVNegetaLep->SetStats(1);
			histPosVNegetaLep->Draw();
			if(SaveAsRootFile){
				TFile outfile13(rootOutputFileName13, "RECREATE");
				c13->Write(rootOutputFileName13);
				outfile13.Close();
			}
			
	TCanvas *c14 = new TCanvas("c14","PosVNegphiLep Histrogram" ,200,10,900,700);
		c14->SetFillColor(10);
		c14->SetGrid();
		c14->GetFrame()->SetFillColor(10);
		c14->GetFrame()->SetBorderSize(12);
		c14->Range(0,0,1,1);

		gStyle->SetOptStat(1);
			
			c14->cd();
			histPosVNegphiLep->SetXTitle("Pos");
			histPosVNegphiLep->SetYTitle("Neg");
			histPosVNegphiLep->SetZTitle("counts");
			histPosVNegphiLep->SetStats(1);
			histPosVNegphiLep->Draw();
			if(SaveAsRootFile){
				TFile outfile14(rootOutputFileName14, "RECREATE");
				c14->Write(rootOutputFileName14);
				outfile14.Close();
			}


	
	/////////// The Momentum Hists ////////////		
		
	TCanvas *c16 = new TCanvas("c16","Neg Momentum Prime Histrogram" ,200,10,900,700);
		c16->Divide(2,2);
		c16->SetFillColor(10);
		c16->SetGrid();
		c16->GetFrame()->SetFillColor(10);
		c16->GetFrame()->SetBorderSize(12);
		c16->Range(0,0,1,1);

		gStyle->SetOptStat(1);
			
			c16->cd(1);
			histNegPxPrime->SetXTitle("NegPx");
			histNegPxPrime->SetYTitle("counts");
			histNegPxPrime->SetStats(1);
			histNegPxPrime->Draw();
			
			c16->cd(2);
			histNegPyPrime->SetXTitle("NegPy");
			histNegPyPrime->SetYTitle("counts");
			histNegPyPrime->SetStats(1);
			histNegPyPrime->Draw();
			
			c16->cd(3);
			histNegPzPrime->SetXTitle("NegPz");
			histNegPzPrime->SetYTitle("counts");
			histNegPzPrime->SetStats(1);
			histNegPzPrime->Draw();
			
			c16->cd(4);
			histNegEPrime->SetXTitle("NegE");
			histNegEPrime->SetYTitle("counts");
			histNegEPrime->SetStats(1);
			histNegEPrime->Draw();
			
			if(SaveAsRootFile){
				TFile outfile16(rootOutputFileName16, "RECREATE");
				c16->Write(rootOutputFileName16);
				outfile16.Close();
			}
			
	TCanvas *c17 = new TCanvas("c17","PosVNegPrime Histrogram" ,200,10,900,700);
		c17->Divide(2,2);
		c17->SetFillColor(10);
		c17->SetGrid();
		c17->GetFrame()->SetFillColor(10);
		c17->GetFrame()->SetBorderSize(12);
		c17->Range(0,0,1,1);

		gStyle->SetOptStat(1);
			
			c17->cd(1);
			histPosVNegPxPrime->SetXTitle("PosPxPrime");
			histPosVNegPxPrime->SetYTitle("NegPxPrime");
			histPosVNegPxPrime->SetStats(1);
			histPosVNegPxPrime->Draw();
			
			c17->cd(2);
			histPosVNegPyPrime->SetXTitle("PosPyPrime");
			histPosVNegPyPrime->SetYTitle("NegPxPrime");
			histPosVNegPyPrime->SetStats(1);
			histPosVNegPyPrime->Draw();
			
			c17->cd(3);
			histPosVNegPzPrime->SetXTitle("PosPzPrime");
			histPosVNegPzPrime->SetYTitle("NegPxPrime");
			histPosVNegPzPrime->SetStats(1);
			histPosVNegPzPrime->Draw();
			
			c17->cd(4);
			histPosVNegEPrime->SetXTitle("PosEPrime");
			histPosVNegEPrime->SetYTitle("NegPxPrime");
			histPosVNegEPrime->SetStats(1);
			histPosVNegEPrime->Draw();
			
			if(SaveAsRootFile){
				TFile outfile17(rootOutputFileName17, "RECREATE");
				c17->Write(rootOutputFileName17);
				outfile17.Close();
			}
			
	TCanvas *c18 = new TCanvas("c18","Tot Histrogram" ,200,10,900,700);
		c18->Divide(2,2);
		c18->SetFillColor(10);
		c18->SetGrid();
		c18->GetFrame()->SetFillColor(10);
		c18->GetFrame()->SetBorderSize(12);
		c18->Range(0,0,1,1);

		gStyle->SetOptStat(1);
			
			c18->cd(1);
			histTotPx->SetXTitle("TotPx");
			histTotPx->SetYTitle("counts");
			histTotPx->SetStats(1);
			histTotPx->Draw();
			
			c18->cd(2);
			histTotPy->SetXTitle("TotPy");
			histTotPy->SetYTitle("counts");
			histTotPy->SetStats(1);
			histTotPy->Draw();
			
			c18->cd(3);
			histTotPz->SetXTitle("TotPz");
			histTotPz->SetYTitle("counts");
			histTotPz->SetStats(1);
			histTotPz->Draw();
			
			c18->cd(4);
			histTotE->SetXTitle("TotE");
			histTotE->SetYTitle("counts");
			histTotE->SetStats(1);
			histTotE->Draw();
			
			if(SaveAsRootFile){
				TFile outfile18(rootOutputFileName18, "RECREATE");
				c18->Write(rootOutputFileName18);
				outfile18.Close();
			}
			
	TCanvas *c19 = new TCanvas("c19","Tot Prime Histrogram" ,200,10,900,700);
		c19->Divide(2,2);
		c19->SetFillColor(10);
		c19->SetGrid();
		c19->GetFrame()->SetFillColor(10);
		c19->GetFrame()->SetBorderSize(12);
		c19->Range(0,0,1,1);

		gStyle->SetOptStat(1);
			
			c19->cd(1);
			histTotPxPrime->SetXTitle("TotPxPrime");
			histTotPxPrime->SetYTitle("counts");
			histTotPxPrime->SetStats(1);
			histTotPxPrime->Draw();
			
			c19->cd(2);
			histTotPyPrime->SetXTitle("TotPyPrime");
			histTotPyPrime->SetYTitle("counts");
			histTotPyPrime->SetStats(1);
			histTotPyPrime->Draw();
			
			c19->cd(3);
			histTotPzPrime->SetXTitle("TotPzPrime");
			histTotPzPrime->SetYTitle("counts");
			histTotPzPrime->SetStats(1);
			histTotPzPrime->Draw();
			
			c19->cd(4);
			histTotEPrime->SetXTitle("TotEPrime");
			histTotEPrime->SetYTitle("counts");
			histTotEPrime->SetStats(1);
			histTotEPrime->Draw();
			
			if(SaveAsRootFile){
				TFile outfile19(rootOutputFileName19, "RECREATE");
				c19->Write(rootOutputFileName19);
				outfile19.Close();
			}
			
	TCanvas *c20 = new TCanvas("c20","Pos Momentum Histrogram" ,200,10,900,700);
		c20->Divide(2,2);
		c20->SetFillColor(10);
		c20->SetGrid();
		c20->GetFrame()->SetFillColor(10);
		c20->GetFrame()->SetBorderSize(12);
		c20->Range(0,0,1,1);

		gStyle->SetOptStat(1);
			
			c20->cd(1);
			histPosPx->SetXTitle("PosPx");
			histPosPx->SetYTitle("counts");
			histPosPx->SetStats(1);
			histPosPx->Draw();
			
			c20->cd(2);
			histPosPy->SetXTitle("PosPy");
			histPosPy->SetYTitle("counts");
			histPosPy->SetStats(1);
			histPosPy->Draw();
		
			c20->cd(3);
			histPosPz->SetXTitle("PosPz");
			histPosPz->SetYTitle("counts");
			histPosPz->SetStats(1);
			histPosPz->Draw();
			
			c20->cd(4);
			histPosE->SetXTitle("PosE");
			histPosE->SetYTitle("counts");
			histPosE->SetStats(1);
			histPosE->Draw();
			
			if(SaveAsRootFile){
				TFile outfile20(rootOutputFileName20, "RECREATE");
				c20->Write(rootOutputFileName20);
				outfile20.Close();
			}
		
		
	TCanvas *c21 = new TCanvas("c21","Neg Momentum Histrogram" ,200,10,900,700);
		c21->Divide(2,2);
		c21->SetFillColor(10);
		c21->SetGrid();
		c21->GetFrame()->SetFillColor(10);
		c21->GetFrame()->SetBorderSize(12);
		c21->Range(0,0,1,1);

		gStyle->SetOptStat(1);
			
			c21->cd(1);
			histNegPx->SetXTitle("NegPx");
			histNegPx->SetYTitle("counts");
			histNegPx->SetStats(1);
			histNegPx->Draw();
			
			c21->cd(2);
			histNegPy->SetXTitle("NegPy");
			histNegPy->SetYTitle("counts");
			histNegPy->SetStats(1);
			histNegPy->Draw();
			
			c21->cd(3);
			histNegPz->SetXTitle("NegPz");
			histNegPz->SetYTitle("counts");
			histNegPz->SetStats(1);
			histNegPz->Draw();
			
			c21->cd(4);
			histNegE->SetXTitle("NegE");
			histNegE->SetYTitle("counts");
			histNegE->SetStats(1);
			histNegE->Draw();
			
			if(SaveAsRootFile){
				TFile outfile21(rootOutputFileName21, "RECREATE");
				c21->Write(rootOutputFileName21);
				outfile21.Close();
			}
			
	TCanvas *c22 = new TCanvas("c22","PosVNeg Histrogram" ,200,10,900,700);
		c22->Divide(2,2);
		c22->SetFillColor(10);
		c22->SetGrid();
		c22->GetFrame()->SetFillColor(10);
		c22->GetFrame()->SetBorderSize(12);
		c22->Range(0,0,1,1);

		gStyle->SetOptStat(1);
			
			c22->cd(1);
			histPosVNegPx->SetXTitle("PosPx");
			histPosVNegPx->SetYTitle("NegPx");
			histPosVNegPx->SetStats(1);
			histPosVNegPx->Draw();
			
			c22->cd(2);
			histPosVNegPy->SetXTitle("PosPy");
			histPosVNegPy->SetYTitle("NegPx");
			histPosVNegPy->SetStats(1);
			histPosVNegPy->Draw();
			
			c22->cd(3);
			histPosVNegPz->SetXTitle("PosPz");
			histPosVNegPz->SetYTitle("NegPx");
			histPosVNegPz->SetStats(1);
			histPosVNegPz->Draw();
			
			c22->cd(4);
			histPosVNegE->SetXTitle("PosE");
			histPosVNegE->SetYTitle("NegPx");
			histPosVNegE->SetStats(1);
			histPosVNegE->Draw();
			
			if(SaveAsRootFile){
				TFile outfile22(rootOutputFileName22, "RECREATE");
				c22->Write(rootOutputFileName22);
				outfile22.Close();
			}
*/

	return 0;
}