//Program to read CERN data from trees and make hists. (So far)
//V3: Invariant mass plot works!!!
//V4: PosNeg plots, pX, pY, pZ

# 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; //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 TreeName2[50] = "OutputTree;33";
	
	char rootOutputFileName1[50]  = "pTLepHist.root";
	char rootOutputFileName2[50]  = "etaLepHist.root";
	char rootOutputFileName3[50]  = "phiLepHist.root";
	
	char rootOutputFileName4[50]  = "pTLepHistPos.root";
	char rootOutputFileName5[50]  = "etaLepHistPos.root";
	char rootOutputFileName6[50]  = "phiLepHistPos.root";
	
	char rootOutputFileName7[50]  = "pTLepHistNeg.root";
	char rootOutputFileName8[50]  = "etaLepHistNeg.root";
	char rootOutputFileName9[50]  = "phiLepHistNeg.root";
	
	char rootOutputFileName10[50]  = "nLepHist.root";
	char rootOutputFileName11[50] = "InvariantMass.root";
	
	char rootOutputFileName12[50]  = "pTLepPosNeg.root";
	char rootOutputFileName13[50] = "etaLepPosNeg.root";
	char rootOutputFileName14[50]  = "phiLepPosNeg.root";

	
	double xMin = -300;
	double xMax = 300;
	double yMin = -300;
	double yMax = 300;
	// double zMin = 0;
	// double zMax = 10000;
	
	const int xNChannels = 500; //Number of channels for root hist
	const int yNChannels = 500; 
	
	bool SaveFiles = true;
	
	bool IncludeMuons = false;
	bool DoBothTrees = false;
	
	bool OutputToFile = false;
	
	///////////////////////////////////////
	
	/////Varibles (Pre-stated)/////
	
	char rootFilePath[100]; 
	
	int flavLep[3], nLeps35;
	int NumOfPosNegPairs = 0;
	
	Float_t pTLep[3], etaLep[3], phiLep[3];
	
	double InvariantMass;
	
	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 *histpTLepPosNeg = new TH2F("histpTLepPosNeg", "pTLepPosNeg hist", xNChannels, xMin, xMax, yNChannels, yMin, yMax);
	TH2F *histetaLepPosNeg = new TH2F("histetaLepPosNeg", "etaLepPosNeg hist", xNChannels, xMin, xMax, yNChannels, yMin, yMax);
	TH2F *histphiLepPosNeg = new TH2F("histphiLepPosNeg", "phiLepPosNeg hist", xNChannels, xMin, xMax, yNChannels, yMin, yMax);
	
	TH2F *histpX = new TH2F("histpX", "pX hist", xNChannels, xMin, xMax, yNChannels, yMin, yMax);
	TH2F *histpY = new TH2F("histpY", "pY hist", xNChannels, xMin, xMax, yNChannels, yMin, yMax);
	TH2F *histpZ = new TH2F("histpZ", "pZ hist", xNChannels, xMin, xMax, yNChannels, yMin, yMax);
	TH2F *histE = new TH2F("histE", "E hist", xNChannels, xMin, xMax, yNChannels, yMin, yMax);
	
	TH1F *histTotPx = new TH1F("histTotPx", "TotpX hist", xNChannels, xMin, xMax, yNChannels, yMin, yMax);
	TH1F *histTotPy = new TH1F("histTotPy", "TotpY hist", xNChannels, xMin, xMax, yNChannels, yMin, yMax);
	TH1F *histTotPz = new TH1F("histTotPz", "TotpZ hist", xNChannels, xMin, xMax, yNChannels, yMin, yMax);
	TH1F *histTotE = new TH1F("histTotE", "TotE hist", xNChannels, xMin, xMax, yNChannels, yMin, yMax);

	TLorentzVector pos_muon;
	TLorentzVector neg_muon;
	TLorentzVector tot_muon;
	
	//////////////////////////////////////
	
	//Creating path to root file and getting tree//
	strcpy(rootFilePath, rootFilePrefix);
	strcat(rootFilePath, rootFileName);
	cout<<"Root file path: "<<rootFilePath<<endl;
	
	TFile *file = new TFile(rootFilePath);
	TTree *t1 = (TTree*)file->Get(TreeName1);
	
	//TFile *file = new TFile(rootFilePath1);
//	TTree *t1 = (TTree*)file->Get(TreeName2);
	
	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 NEntries = t1->GetEntries();
	cout<<"NEntries: "<<NEntries<<endl;
	
	if (!IncludeMuons) {
		for(int i = 0; i < NEntries; i++) {
			t1->GetEntry(i);
			if(nLeps35 >= 0) {
				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 pTLepNeg, etaLepNeg, phiLepNeg;
					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_muon.Px();
						//pos_muon.Py();
						//pos_muon.Pz();
					}
					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_muon.Px();
						//neg_muon.Py();
						//neg_muon.Pz();
					}
					histpX->Fill(pos_muon.Px(),neg_muon.Px());
					histpY->Fill(pos_muon.Py(),neg_muon.Py());
					histpZ->Fill(pos_muon.Pz(),neg_muon.Pz());
					histE->Fill(pos_muon.E(),neg_muon.E());
					
					histpTLepPosNeg->Fill(pTLepPos, pTLepNeg);
					histetaLepPosNeg->Fill(etaLepPos, etaLepNeg);
					histphiLepPosNeg->Fill(phiLepPos, phiLepNeg);
				}
				if (flavLep[0] > 0 && flavLep[1] < 0){
					tot_muon = pos_muon + neg_muon;
					InvariantMass = tot_muon.M();
					histM->Fill(InvariantMass);
					
					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();
					histTotPx->Fill(tot_Px);
					histTotPy->Fill(tot_Py);		
					histTotPz->Fill(tot_Pz);	
					histTotE->Fill(tot_E);
					
					NumOfPosNegPairs++;
				}
				if (flavLep[0] < 0 && flavLep[1] > 0){
					tot_muon = pos_muon + neg_muon;
					InvariantMass = tot_muon.M();
					histM->Fill(InvariantMass);
					
					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();
					histTotPx->Fill(tot_Px);
					histTotPy->Fill(tot_Py);		
					histTotPz->Fill(tot_Pz);	
					histTotE->Fill(tot_E);
					
					NumOfPosNegPairs++;
				}
			}
		}
	}
	
	if (IncludeMuons) {
		for(int i = 0; i < NEntries; 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 pTLepNeg, etaLepNeg, phiLepNeg;
					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];
					}
					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];
					}
					histpX->Fill(pos_muon.Px(),neg_muon.Px());
					histpY->Fill(pos_muon.Py(),neg_muon.Py());
					histpZ->Fill(pos_muon.Pz(),neg_muon.Pz());
					histE->Fill(pos_muon.E(),neg_muon.E());
					
					histpTLepPosNeg->Fill(pTLepPos, pTLepNeg);
					histetaLepPosNeg->Fill(etaLepPos, etaLepNeg);
					histphiLepPosNeg->Fill(phiLepPos, phiLepNeg);
				}
				if (flavLep[0] == 2 && flavLep[1] == -2){
					tot_muon = pos_muon + neg_muon;
					InvariantMass = tot_muon.M();
					histM->Fill(InvariantMass);
					
					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();
					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, 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();
					histTotPx->Fill(tot_Px);
					histTotPy->Fill(tot_Py);		
					histTotPz->Fill(tot_Pz);	
					histTotE->Fill(tot_E);
					
					NumOfPosNegPairs++;
				}
			}
		}
	}
	
	cout << "NumOfPosNegPairs: " << NumOfPosNegPairs << endl;
	
	/*
	//Loop to get entries and fill hists//
	for(int i = 0; i < NEntries; i++){
		t1->GetEntry(i);
		cout<<"i: "<<i<<"  "<<"nLep: "<<nLeps35<<endl;
		if(nLeps35 == 2){
			for(int j = 0; j < nLeps35; j++){
				//t1->GetEvent(j);
				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]);
				}
				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]);
				}
				//cout<<"i: "<<i<<"  "<<"j: "<<j<<"  "<<"ptLep: "<<pTLep[j]<<endl;
			}
			histnLep->Fill(nLeps35);
		} 
	}
	*/
	 
	
	//////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();
			
			
	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();


	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();
			
			
			
	// TCanvas *c4 = new TCanvas("c1","pTLep 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);
			
			// c1->cd();
			// histpTLep->SetXTitle("pTLep");
			// histpTLep->SetYTitle("Counts");
			// histpTLep->SetStats(1);
			// histpTLep->Draw();
			
			
	// 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();


	// 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();
			
			
	// TCanvas *c7 = new TCanvas("c7","nLep Histrogram" ,100,5,400,300);
		// c7->SetFillColor(10);
		// c7->SetGrid();
		// c7->GetFrame()->SetFillColor(10);
		// c7->GetFrame()->SetBorderSize(12);
		// c7->Range(0,0,1,1);

		// gStyle->SetOptStat(1);
			
			// c7->cd();
			// histnLep->SetXTitle("nLep");
			// histnLep->SetYTitle("Counts");
			// histnLep->SetStats(1);
			// histnLep->Draw();
			
	TCanvas *c11 = new TCanvas("c11","Invariant Mass Histogram" ,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 (SaveFiles) {
				TFile outputfile11(rootOutputFileName11, "RECREATE");
				c11->Write(rootOutputFileName11);
				outputfile11.Close();
			}
			
	TCanvas *c12 = new TCanvas("c12","pTLepPosNeg Histogram" ,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();
		histpTLepPosNeg->SetXTitle("Pos");
		histpTLepPosNeg->SetYTitle("Neg");
		histpTLepPosNeg->SetZTitle("counts");
		histpTLepPosNeg->SetStats(1);
		histpTLepPosNeg->Draw();
		if (SaveFiles) {
			TFile outputfile12(rootOutputFileName12, "RECREATE");
			c12->Write(rootOutputFileName12);
			outputfile12.Close();
			}
			
	TCanvas *c13 = new TCanvas("c13","etaLepPosNeg Histogram" ,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();
		histetaLepPosNeg->SetXTitle("Pos");
		histetaLepPosNeg->SetYTitle("Neg");
		histetaLepPosNeg->SetZTitle("counts");
		histetaLepPosNeg->SetStats(1);
		histetaLepPosNeg->Draw();
		if (SaveFiles) {
			TFile outputfile13(rootOutputFileName13, "RECREATE");
			c13->Write(rootOutputFileName13);
			outputfile13.Close();
			}
			
	TCanvas *c14 = new TCanvas("c14","phiLepPosNeg Histogram" ,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();
		histphiLepPosNeg->SetXTitle("Pos");
		histphiLepPosNeg->SetYTitle("Neg");
		histphiLepPosNeg->SetZTitle("counts");
		histphiLepPosNeg->SetStats(1);
		histphiLepPosNeg->Draw();
		if (SaveFiles) {
			TFile outputfile14(rootOutputFileName14, "RECREATE");
			c14->Write(rootOutputFileName14);
			outputfile14.Close();
			}
			
	TCanvas *c15 = new TCanvas("c15","pXPosNeg Histogram" ,200,10,900,700);
	c15->SetFillColor(10);
	c15->SetGrid();
	c15->GetFrame()->SetFillColor(10);
	c15->GetFrame()->SetBorderSize(12);
	c15->Range(0,0,1,1);

	gStyle->SetOptStat(1);
		
		c15->cd();
		histpX->SetXTitle("Pos");
		histpX->SetYTitle("Neg");
		histpX->SetZTitle("counts");
		histpX->SetStats(1);
		histpX->Draw();
		if (!SaveFiles) {
			//TFile outputfile15(rootOutputFileName15, "RECREATE");
			//c15->Write(rootOutputFileName15);
			//outputfile15.Close();
			}
			
	TCanvas *c16 = new TCanvas("c16","pYPosNeg Histogram" ,200,10,900,700);
	c16->SetFillColor(10);
	c16->SetGrid();
	c16->GetFrame()->SetFillColor(10);
	c16->GetFrame()->SetBorderSize(12);
	c16->Range(0,0,1,1);

	gStyle->SetOptStat(1);
		
		c16->cd();
		histpY->SetXTitle("Pos");
		histpY->SetYTitle("Neg");
		histpY->SetZTitle("counts");
		histpY->SetStats(1);
		histpY->Draw();
		if (!SaveFiles) {
			//TFile outputfile16(rootOutputFileName16, "RECREATE");
			//c16->Write(rootOutputFileName16);
			//outputfile16.Close();
			}
			
	TCanvas *c17 = new TCanvas("c17","pZPosNeg Histogram" ,200,10,900,700);
	c17->SetFillColor(10);
	c17->SetGrid();
	c17->GetFrame()->SetFillColor(10);
	c17->GetFrame()->SetBorderSize(12);
	c17->Range(0,0,1,1);

	gStyle->SetOptStat(1);
		
		c17->cd();
		histpZ->SetXTitle("Pos");
		histpZ->SetYTitle("Neg");
		histpZ->SetZTitle("counts");
		histpZ->SetStats(1);
		histpZ->Draw();
		if (!SaveFiles) {
			//TFile outputfile17(rootOutputFileName17, "RECREATE");
			//c19->Write(rootOutputFileName17);
			//outputfile17.Close();
			}
			
	TCanvas *c18 = new TCanvas("c18","EPosNeg Histogram" ,200,10,900,700);
	c18->SetFillColor(10);
	c18->SetGrid();
	c18->GetFrame()->SetFillColor(10);
	c18->GetFrame()->SetBorderSize(12);
	c18->Range(0,0,1,1);

	gStyle->SetOptStat(1);
		
		c18->cd();
		histE->SetXTitle("Pos");
		histE->SetYTitle("Neg");
		histE->SetZTitle("counts");
		histE->SetStats(1);
		histE->Draw();
		if (!SaveFiles) {
			//TFile outputfile17(rootOutputFileName17, "RECREATE");
			//c18->Write(rootOutputFileName17);
			//outputfile17.Close();
			}
			
	TCanvas *c19 = new TCanvas("c19","4PosNeg Histogram" ,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_1->cd();
		histpX->SetXTitle("PosPx");
		histpX->SetYTitle("NegPx");
		histpX->SetZTitle("counts");
		histpX->SetStats(1);
		histpX->Draw();
		
		c19_2->cd();
		histpY->SetXTitle("PosPy");
		histpY->SetYTitle("NegPy");
		histpY->SetZTitle("counts");
		histpY->SetStats(1);
		histpY->Draw();
		
		c19_3->cd();
		histpZ->SetXTitle("PosPz");
		histpZ->SetYTitle("NegPz");
		histpZ->SetZTitle("counts");
		histpZ->SetStats(1);
		histpZ->Draw();
		
		c19_4->cd();
		histE->SetXTitle("PosE");
		histE->SetYTitle("NegE");
		histE->SetZTitle("counts");
		histE->SetStats(1);
		histE->Draw();
		
	TCanvas *c20 = new TCanvas("c20","Tot Histogram" ,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_1->cd();
		histpX->SetXTitle("PosPx");
		histpX->SetYTitle("NegPx");
//		histpX->SetZTitle("counts");
		histpX->SetStats(1);
		histpX->Draw();
		
		c20_2->cd();
		histpY->SetXTitle("PosPy");
		histpY->SetYTitle("NegPy");
//		histpY->SetZTitle("counts");
		histpY->SetStats(1);
		histpY->Draw();
		
		c20_3->cd();
		histpZ->SetXTitle("PosPz");
		histpZ->SetYTitle("NegPz");
//		histpZ->SetZTitle("counts");
		histpZ->SetStats(1);
		histpZ->Draw();
		
		c20_4->cd();
		histE->SetXTitle("PosE");
		histE->SetYTitle("NegE");
//		histE->SetZTitle("counts");
		histE->SetStats(1);
		histE->Draw();

	return 0;
}