//Program to read data from root file and make Invariant Mass hist
// Based on: https://mhance.scipp.ucsc.edu/analysisintro.php
//V3: Invariant mass plot works!!!
//V4: PosNeg plots, pX, pY, pZ
//V5: Plot Total muon, divide canvas
//https://root.cern.ch/root/roottalk/roottalk99/2970.html
//V6: Implemented output of PxPyPzE and momentum histograms


# 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 <dos.h> //For Sleep() 

# 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; // No longer used, entries are grabbed automatically now
	
	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] = "InvariantMass.root"; //Name of the output file for Invariant 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] = "PosMomentum.root"; //Name of the output file for nLep hist
	char rootOutputFileName16[50] = "NegMomentum.root"; //Name of the output file for nLep hist
	char rootOutputFileName17[50] = "PosVNegMomentum.root"; //Name of the output file for nLep hist
	char rootOutputFileName18[50] = "TotMomentum.root"; //Name of the output file for nLep hist
	char rootOutputFileName19[50] = "TotMomentumPrime.root"; //Name of the output file for nLep hist
	
	double xMin = -300;
	double xMax = 300;
	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 = true; //If true, will save canvases to root file instead of printing to screen.
	bool IncludeOnlyMuons = false; //If true, will include flavLep = +-1 as well. ie flavLep<0, >0
	
	///////////////////////////////////////
	
	/////Variables/////
	
	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 InvariantMass;
	
	double c = 1;// 2.99792458e8; // speed of light!!!
	
	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", "Invariant Mass hist", xNChannels, 0, 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 *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);
	
	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 *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 *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 *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);
	
	TLorentzVector pos_muon;
	TLorentzVector neg_muon;
	TLorentzVector tot_muon;
	
	//////////////////////////////////////
	
	//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;
	
	
	//NEntries = 2; // for testing purposes
	
	// This includes leptons ONLY//
	if(!IncludeOnlyMuons){
		for(int i = 0; i < NEntries1; i++){
			t1->GetEntry(i);
			if(nLeps35 >= 2){
				histnLep->Fill(nLeps35);
				cout<<"nLeps35: "<<nLeps35<<endl;
				for(int j = 0; j < 2; j++){
					cout<<" i: "<<i<<" j: "<<j<<" flavLep: "<<flavLep[j]<<endl;
					double pTLepPos, etaLepPos, phiLepPos;
					double pos_Px, pos_Py, pos_Pz, pos_E;
					double pTLepNeg, etaLepNeg, phiLepNeg;
					double neg_Px, neg_Py, neg_Pz, neg_E;
					if(flavLep[j] > 0){
						pos_muon.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 = pos_muon.Px();
						pos_Py = pos_muon.Py();
						pos_Pz = pos_muon.Pz();
						pos_E = pos_muon.E();
						histPosPx->Fill(pos_Px);
						histPosPy->Fill(pos_Py);
						histPosPz->Fill(pos_Pz);
						histPosE->Fill(pos_E);
					}
					if(flavLep[j] < 0){
						neg_muon.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 = neg_muon.Px();
						neg_Py = neg_muon.Py();
						neg_Pz = neg_muon.Pz();
						neg_E = neg_muon.E();
						histNegPx->Fill(neg_Px);
						histNegPy->Fill(neg_Py);
						histNegPz->Fill(neg_Pz);
						histNegE->Fill(neg_E); */
						
						double neg_Px_prime, neg_Py_prime, neg_Pz_prime, neg_E_prime;
						neg_Px_prime = neg_muon.Px();
						neg_Py_prime = neg_muon.Py();
						neg_Pz_prime = neg_muon.Pz();
						neg_E_prime = neg_muon.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);
					
					histPosVNegPx->Fill(pos_Px, neg_Px);
					histPosVNegPy->Fill(pos_Py, neg_Py);
					histPosVNegPz->Fill(pos_Pz, neg_Pz);
					histPosVNegE->Fill(pos_E, neg_E);
				}
				if(flavLep[0] > 0 && flavLep[1] < 0){
					tot_muon = pos_muon + neg_muon; 
					InvariantMass = tot_muon.M();
					histM->Fill(InvariantMass);
					
					double tot_Px_prime, tot_Py_prime, tot_Pz_prime, tot_E_prime;
					double tot_Px, tot_Py, tot_Pz, tot_E;
					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(); // returns beta (velocity / speed of light)
					double gamma = tot_muon.Gamma(); // returns 1/Sqrt(1-beta*beta)
					double P_prime = tot_muon.P(); // returns magnitude of three vector component
					double V = beta*c; 
					double Vx = tot_Px_prime/(InvariantMass*gamma); // calculate vx, knowing p = m*v/gamma
					double Vy = tot_Py_prime/(InvariantMass*gamma); // calculate vy, knowing p = m*v/gamma
					double Vz = tot_Pz_prime/(InvariantMass*gamma); // calculate vz, knowing p = m*v/gamma
	
					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 << endl;
					cout<<"tot_Py_prime: " << tot_Py_prime << endl;
					cout<<"tot_Pz_prime: " << tot_Pz_prime << endl;
					cout<<"tot_E_prime: " << tot_Px_prime << endl;
					
					NumOfPosNegPairs++;
				}
				if(flavLep[0] < 0 && flavLep[1] > 0){
					tot_muon = pos_muon + neg_muon; 
					InvariantMass = tot_muon.M();
					histM->Fill(InvariantMass);
					
					double tot_Px_prime, tot_Py_prime, tot_Pz_prime, tot_E_prime;
					double tot_Px, tot_Py, tot_Pz, tot_E;
					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(); // returns beta (velocity / speed of light)
					double gamma = tot_muon.Gamma(); // returns 1/Sqrt(1-beta*beta)
					double P_prime = tot_muon.P(); // returns magnitude of three vector component
					double V = beta*c; 
					double Vx = tot_Px_prime/(InvariantMass*gamma); // calculate vx, knowing p = m*v/gamma
					double Vy = tot_Py_prime/(InvariantMass*gamma); // calculate vy, knowing p = m*v/gamma
					double Vz = tot_Pz_prime/(InvariantMass*gamma); // calculate vz, knowing p = m*v/gamma
	
					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);
					
					NumOfPosNegPairs++;
				}
				
					
					double neg_P_prime = neg_muon.P(); // returns magnitude of three vector components
					double neg_Px, neg_Py, neg_Pz, neg_E;
					neg_E = (neg_E_prime/gamma) + V*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.P(); // returns magnitude of three vector components
					double pos_Px, pos_Py, pos_Pz, pos_E;
					pos_E = (pos_E_prime/gamma) + V*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);
			}
		}
	}
	
	//ONLY MUONS//
	if(IncludeOnlyMuons){
		for(int i = 0; i < NEntries1; i++){
			t1->GetEntry(i);
			if(nLeps35 >= 2){
				histnLep->Fill(nLeps35);
				cout<<"nLeps35: "<<nLeps35<<endl;
				for(int j = 0; j < 2; j++){
					cout<<" i: "<<i<<" j: "<<j<<" flavLep: "<<flavLep[j]<<endl;
					double pTLepPos, etaLepPos, phiLepPos;
					double pos_Px, pos_Py, pos_Pz, pos_E;
					double pTLepNeg, etaLepNeg, phiLepNeg;
					double neg_Px, neg_Py, neg_Pz, neg_E;
					if(flavLep[j] == 2){
						pos_muon.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 = pos_muon.Px();
						pos_Py = pos_muon.Py();
						pos_Pz = pos_muon.Pz();
						pos_E = pos_muon.E();
						histPosPx->Fill(pos_Px);
						histPosPy->Fill(pos_Py);
						histPosPz->Fill(pos_Pz);
						histPosE->Fill(pos_E);
					}
					if(flavLep[j] == -2){
						neg_muon.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 = neg_muon.Px();
						neg_Py = neg_muon.Py();
						neg_Pz = neg_muon.Pz();
						neg_E = neg_muon.E();
						histNegPx->Fill(neg_Px);
						histNegPy->Fill(neg_Py);
						histNegPz->Fill(neg_Pz);
						histNegE->Fill(neg_E);
					}
					histPosVNegpTLep->Fill(pTLepPos, pTLepNeg);
					histPosVNegetaLep->Fill(etaLepPos, etaLepNeg);
					histPosVNegphiLep->Fill(phiLepPos, phiLepNeg);
					
					histPosVNegPx->Fill(pos_Px, neg_Px);
					histPosVNegPy->Fill(pos_Py, neg_Py);
					histPosVNegPz->Fill(pos_Pz, neg_Pz);
					histPosVNegE->Fill(pos_E, neg_E);
				}
				
				if(flavLep[0] == 2 && flavLep[1] == -2){
					tot_muon = pos_muon + neg_muon; 
					InvariantMass = tot_muon.M();
					histM->Fill(InvariantMass);
					
					double tot_Px_prime, tot_Py_prime, tot_Pz_prime, tot_E_prime;
					double tot_Px, tot_Py, tot_Pz, tot_E;
					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(); // returns beta (velocity / speed of light)
					double gamma = tot_muon.Gamma(); // returns 1/Sqrt(1-beta*beta)
					double P_prime = tot_muon.P(); // returns magnitude of three vector component
					double V = beta*c; 
					double Vx = tot_Px_prime/(Z0mass*gamma); // calculate vx, knowing p = m*v/gamma
					double Vy = tot_Py_prime/(Z0mass*gamma); // calculate vy, knowing p = m*v/gamma
					double Vz = tot_Pz_prime/(Z0mass*gamma); // calculate vz, knowing p = m*v/gamma
	
					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);
					
					NumOfPosNegPairs++;
				}
				if(flavLep[0] == -2 && flavLep[1] == 2){
					tot_muon = pos_muon + neg_muon; 
					InvariantMass = tot_muon.M();
					histM->Fill(InvariantMass);
					
					double tot_Px_prime, tot_Py_prime, tot_Pz_prime, tot_E_prime;
					double tot_Px, tot_Py, tot_Pz, tot_E;
					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(); // returns beta (velocity / speed of light)
					double gamma = tot_muon.Gamma(); // returns 1/Sqrt(1-beta*beta)
					double P_prime = tot_muon.P(); // returns magnitude of three vector component
					double V = beta*c; 
					double Vx = tot_Px_prime/(Z0mass*gamma); // calculate vx, knowing p = m*v/gamma
					double Vy = tot_Py_prime/(Z0mass*gamma); // calculate vy, knowing p = m*v/gamma
					double Vz = tot_Pz_prime/(Z0mass*gamma); // calculate vz, knowing p = m*v/gamma
	
					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);
					
					NumOfPosNegPairs++;
				}
			}
		}
	}
	
	cout<<"NumOfPosNegPairs: "<<NumOfPosNegPairs<<endl; 
	
	//////Graphics//////
	cout<<"Making histrograms..."<<endl;
	
	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","Invariant Mass Histrogram" ,200,10,900,700);
		c11->SetFillColor(10);
		c11->SetGrid();
		c11->GetFrame()->SetFillColor(10);
		c11->GetFrame()->SetBorderSize(12);
		c11->Range(0,0,1,1);

		gStyle->SetOptStat(1);
			
			c11->cd();
			histM->SetXTitle("Invariant Mass (GeV)");
			histM->SetYTitle("Counts");
			histM->SetStats(1);
			histM->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();
			}
			
	// MOMENTUM HISTOGRAMS
	// V5: Dividng canvases to include 4 plots: Px, Py, Pz, and E
	
	TCanvas *c15 = new TCanvas("c15","Pos Momentum 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 *c16 = new TCanvas("c16","Neg Momentum 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);
			histNegPx->SetXTitle("NegPx");
			histNegPx->SetYTitle("counts");
			histNegPx->SetStats(1);
			histNegPx->Draw();
			
			c16->cd(2);
			histNegPy->SetXTitle("NegPy");
			histNegPy->SetYTitle("counts");
			histNegPy->SetStats(1);
			histNegPy->Draw();
			
			c16->cd(3);
			histNegPz->SetXTitle("NegPz");
			histNegPz->SetYTitle("counts");
			histNegPz->SetStats(1);
			histNegPz->Draw();
			
			c16->cd(4);
			histNegE->SetXTitle("NegE");
			histNegE->SetYTitle("counts");
			histNegE->SetStats(1);
			histNegE->Draw();
			
			if(SaveAsRootFile){
				TFile outfile16(rootOutputFileName16, "RECREATE");
				c16->Write(rootOutputFileName16);
				outfile16.Close();
			}
			
	TCanvas *c17 = new TCanvas("c17","PosVNeg 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);
			histPosVNegPx->SetXTitle("PosPx");
			histPosVNegPx->SetYTitle("NegPx");
			histPosVNegPx->SetStats(1);
			histPosVNegPx->Draw();
			
			c17->cd(2);
			histPosVNegPy->SetXTitle("PosPy");
			histPosVNegPy->SetYTitle("NegPx");
			histPosVNegPy->SetStats(1);
			histPosVNegPy->Draw();
			
			c17->cd(3);
			histPosVNegPz->SetXTitle("PosPz");
			histPosVNegPz->SetYTitle("NegPx");
			histPosVNegPz->SetStats(1);
			histPosVNegPz->Draw();
			
			c17->cd(4);
			histPosVNegE->SetXTitle("PosE");
			histPosVNegE->SetYTitle("NegPx");
			histPosVNegE->SetStats(1);
			histPosVNegE->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();
			}
		
	// to get the total histograms for the new frame
	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();
			}

	return 0;
}