#define Selector_TrackEff_cxx
// The class definition in Selector_TrackEff.h has been generated automatically
// by the ROOT utility TTree::MakeSelector(). This class is derived
// from the ROOT class TSelector. For more information on the TSelector
// framework see $ROOTSYS/README/README.SELECTOR or the ROOT User Manual.

// The following methods are defined in this file:
//    Begin():        called every time a loop on the tree starts,
//                    a convenient place to create your histograms.
//    SlaveBegin():   called after Begin(), when on PROOF called only on the
//                    slave servers.
//    Process():      called for each event, in this function you decide what
//                    to read and fill your histograms.
//    SlaveTerminate: called at the end of the loop on the tree, when on PROOF
//                    called only on the slave servers.
//    Terminate():    called at the end of the loop on the tree,
//                    a convenient place to draw/fit your histograms.
//
// To use this file, try the following session on your Tree T:
//
// Root > T->Process("Selector_TrackEff.C")
// Root > T->Process("Selector_TrackEff.C","some options")
// Root > T->Process("Selector_TrackEff.C+")
//

#include "Selector_TrackEff.h"

bool gDebugFlag = false;

void Selector_TrackEff::Begin(TTree * /*tree*/)
{
	// The Begin() function is called at the start of the query.
	// When running with PROOF Begin() is only called on the client.
	// The tree argument is deprecated (on PROOF 0 is passed).

	TString option = GetOption();
}

void Selector_TrackEff::SlaveBegin(TTree * /*tree*/)
{
	// The SlaveBegin() function is called after the Begin() function.
	// When running with PROOF SlaveBegin() is called on each slave server.
	// The tree argument is deprecated (on PROOF 0 is passed).

	TString option = GetOption();
}

void Selector_TrackEff::Init(TTree *tree)
{
	Setup_Branches(tree);

	//Get reaction, run information from tree metadata
	TChain* locChain = dynamic_cast<TChain*>(tree);
	TList* locUserInfo = (locChain != NULL) ? locChain->GetTree()->GetUserInfo() : tree->GetUserInfo();
	TMap* locMiscInfoMap = (TMap*)locUserInfo->FindObject("MiscInfoMap");

	//Run Number
	{
		TObjString* locObjString = (TObjString*)locMiscInfoMap->GetValue("RunNumber");
		istringstream locIStream;
		locIStream.str(locObjString->GetName());
		locIStream >> dRunNumber;
	}
	Set_BeamBunchPeriod();

	if(dSetupFlag)
		return; //no need to do the rest

	//Missing PID
	{
		TObjString* locObjString = (TObjString*)locMiscInfoMap->GetValue("MissingPID_PDG");
		int locPDGPID;
		istringstream locIStream;
		locIStream.str(locObjString->GetName());
		locIStream >> locPDGPID;
		dMissingPID = PDGtoPType(locPDGPID);
	}

	//Is 4pi, 3pi
	{
		TObjString* locObjString = (TObjString*)locMiscInfoMap->GetValue("ReactionName");
		dIs4PiFlag = (locObjString->GetString().Contains("4pi") || locObjString->GetString().Contains("4Pi"));
		dIs3PiFlag = (locObjString->GetString().Contains("3pi") || locObjString->GetString().Contains("3Pi"));;
	}

	//output file
	string locOutputFileName = string("trackeff_") + string(ParticleType(dMissingPID));
	if(dIs4PiFlag)
		locOutputFileName += "_4pi.root";
	else if(dIs3PiFlag)
		locOutputFileName += "_3pi.root";
	else
		locOutputFileName += "_2pi.root";

	if(gProofServ == NULL)
	{
		dOutputFile = new TFile(locOutputFileName.c_str(), "RECREATE");
		dProofFile = nullptr;
	}
	else
	{
		dProofFile = new TProofOutputFile(locOutputFileName.c_str(), TProofOutputFile::kMerge);
		dOutputFile = dProofFile->OpenFile("RECREATE");
	}
	dOutputFile->cd();

	cout << "runno, beam period, missing pid, is 4pi, outputfilename = " << dRunNumber << ", " << dBeamBunchPeriod << ", " << ParticleType(dMissingPID) << ", " << dIs4PiFlag << ", " << locOutputFileName << endl;

	//setup histograms
	Setup_Histograms();

	//read in cuts
	Setup_Cuts();

	dSetupFlag = true;
}

void Selector_TrackEff::Set_BeamBunchPeriod(void)
{
	//CCDB environment must be setup!!

	//Pipe the current constant into this function
	ostringstream locCommandStream;
	locCommandStream << "ccdb dump PHOTON_BEAM/RF/beam_period -r " << dRunNumber;
	FILE* locInputFile = gSystem->OpenPipe(locCommandStream.str().c_str(), "r");
	if(locInputFile == NULL)
	{
		dBeamBunchPeriod = -1.0;
		return;
	}

	//get the first line
	char buff[1024]; // I HATE char buffers
	if(fgets(buff, sizeof(buff), locInputFile) == NULL)
	{
		dBeamBunchPeriod = -1.0;
		gSystem->ClosePipe(locInputFile);
		return;
	}

	//get the second line (where the # is)
	if(fgets(buff, sizeof(buff), locInputFile) == NULL)
	{
		dBeamBunchPeriod = -1.0;
		gSystem->ClosePipe(locInputFile);
		return;
	}
	istringstream locStringStream(buff);

	//extract it
	if(!(locStringStream >> dBeamBunchPeriod))
		dBeamBunchPeriod = -1.0;

	//Close the pipe
	gSystem->ClosePipe(locInputFile);
}

/******************************************************************** SETUP HISTOGRAMS ********************************************************************/

void Selector_TrackEff::Setup_Histograms(void)
{
	//RF beam bunches //min/max #peaks range
	dBeamBunchRange_RFSignal = make_pair(0, 0); //e.g. 0, 0 for center, 0, 2 for middle 5
	dBeamBunchRange_RFSideband = make_pair(1, 1); //e.g. 1, 1 for +/-1 (or 2, 4 for +/-(2 -> 4))

	dMinVertexZ = 50.0;
	dMaxVertexZ = 80.0;
	dVertexZBinSize = 30.0;

	if(dMissingPID == Proton)
	{
		dNumMassBins = 3500;
		dNum2DMassBins = 1750;
		dMinMissingMass = -0.5;
		dMaxMissingMass = 3.0;
	}
	else
	{
		dNumMassBins = 2000;
		dNum2DMassBins = 2000;
		dMinMissingMass = -1.0;
		dMaxMissingMass = 1.0;
	}

	dMinBeamEnergy = 3.0;
	dMaxBeamEnergy = 12.0;
	dNumBeamEnergyBins = 900;
	dNum2DBeamEnergyBins = 450;

	dFitBeamEnergyBinSize = 1.0;

	//DERIVED QUANTITIES (DO NOT EDIT THE BELOW!)
	int locNumBunches_RFSideband = 2*(dBeamBunchRange_RFSideband.second - dBeamBunchRange_RFSideband.first + 1);
	int locNumBunches_RFSignal = 2*(dBeamBunchRange_RFSignal.second - dBeamBunchRange_RFSignal.first) + 1;
	dRFSidebandWeight = -1.0*double(locNumBunches_RFSignal)/double(locNumBunches_RFSideband);

	dNumFitBeamEnergyBins = (dMaxBeamEnergy - dMinBeamEnergy) / dFitBeamEnergyBinSize;

	dNumVertexZBins = int((dMaxVertexZ - dMinVertexZ)/dVertexZBinSize + 0.01);

	//Everything is done within E-gamma bins:

	//Resolution distros: sideband-subtract RF, determine +/- 3sigma cuts
	//Missing mass distros (cases: track found/missing): sideband subtract rf, determine +/- 3sigma cuts
	//Kinematic distros (cases: track found/missing, mass signal/sideband): sideband subtract rf
		//Then: Kinematic distros (cases: track found/missing): sideband subtract missing mass
		//Take ratio of found/missing for each kinematic distro plot in each e-gamma bin

	Create_MassHists();
	//Create_MassCuts();
	//Create_MatchCuts();
	Create_MatchingHists();
	Create_ResolutionHists();
	Create_MissingP3SigmaHists();
	Create_EfficiencyHists(true);
	Create_EfficiencyHists(false);

	dComboCounter = 0;
	dFinalizeFlag = false;
}

void Selector_TrackEff::Create_MassHists(void)
{
	string locHistName, locHistTitle;
	string locHistParticleName = ParticleName_ROOT(dMissingPID);

	//Beam RF DeltaT
	locHistName = "BeamRFDeltaT";
	locHistTitle = locHistParticleName + string(";#Deltat (Beam - RF) (ns)");
	dHist_BeamRFDeltaT = new TH1I(locHistName.c_str(), locHistTitle.c_str(), 1200, -6.0, 6.0);

	gDirectory->mkdir("BeamEnergy", "BeamEnergy");
	gDirectory->cd("BeamEnergy");

	//Beam Energy
	locHistName = "BeamEnergy";
	locHistTitle = locHistParticleName + string(";Beam Energy (GeV)");
	dHist_BeamEnergy = new TH1D(locHistName.c_str(), locHistTitle.c_str(), dNumBeamEnergyBins, dMinBeamEnergy, dMaxBeamEnergy);

	//Beam Energy: Found
	locHistName = "BeamEnergy_Found";
	locHistTitle = locHistParticleName + string(" Center RF Bin, Track Found;Beam Energy (GeV)");
	dHist_BeamEnergy_Found = new TH1D(locHistName.c_str(), locHistTitle.c_str(), dNumBeamEnergyBins, dMinBeamEnergy, dMaxBeamEnergy);

	//Beam Energy: Missing
	locHistName = "BeamEnergy_Missing";
	locHistTitle = locHistParticleName + string(" Center RF Bin, Track Missing;Beam Energy (GeV)");
	dHist_BeamEnergy_Missing = new TH1D(locHistName.c_str(), locHistTitle.c_str(), dNumBeamEnergyBins, dMinBeamEnergy, dMaxBeamEnergy);

	gDirectory->cd("..");
	gDirectory->mkdir("MissingMass", "MissingMass");
	gDirectory->cd("MissingMass");

//	if(dMissingPID != Proton)
//	{
		//Missing Mass Squared
		locHistName = "MissingMass";
		locHistTitle = locHistParticleName + string(";Missing Mass Squared (GeV/c^{2})^{2}");
		dHist_MissingMass = new TH1D(locHistName.c_str(), locHistTitle.c_str(), dNumMassBins, dMinMissingMass, dMaxMissingMass);

		// Missing Mass Squared vs Beam Energy
		locHistName = "MissingMassVsBeamEnergy";
		locHistTitle = locHistParticleName + string(";Beam Energy (GeV); Missing Mass Squared (GeV/c^{2})^{2}");
		dHist_MissingMassVsBeamEnergy = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DBeamEnergyBins, dMinBeamEnergy, dMaxBeamEnergy, dNum2DMassBins, dMinMissingMass, dMaxMissingMass);

		// Missing Mass Squared vs Beam Energy: Found
		locHistName = "MissingMassVsBeamEnergy_Found";
		locHistTitle = locHistParticleName + string(" Found;Beam Energy (GeV); Missing Mass Squared (GeV/c^{2})^{2}");
		dHist_MissingMassVsBeamEnergy_Found = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DBeamEnergyBins, dMinBeamEnergy, dMaxBeamEnergy, dNum2DMassBins, dMinMissingMass, dMaxMissingMass);

		// Missing Mass Squared vs Beam Energy: Missing
		locHistName = "MissingMassVsBeamEnergy_Missing";
		locHistTitle = locHistParticleName + string(" Missing;Beam Energy (GeV); Missing Mass Squared (GeV/c^{2})^{2}");
		dHist_MissingMassVsBeamEnergy_Missing = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DBeamEnergyBins, dMinBeamEnergy, dMaxBeamEnergy, dNum2DMassBins, dMinMissingMass, dMaxMissingMass);
/*	}
	else
	{
		//Missing Mass
		locHistName = "MissingMass";
		locHistTitle = locHistParticleName + string(";Missing Mass (GeV/c^{2})");
		dHist_MissingMass = new TH1D(locHistName.c_str(), locHistTitle.c_str(), dNumMassBins, dMinMissingMass, dMaxMissingMass);

		// Missing Mass vs Beam Energy
		locHistName = "MissingMassVsBeamEnergy";
		locHistTitle = locHistParticleName + string(";Beam Energy (GeV); Missing Mass (GeV/c^{2})");
		dHist_MissingMassVsBeamEnergy = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DBeamEnergyBins, dMinBeamEnergy, dMaxBeamEnergy, dNum2DMassBins, dMinMissingMass, dMaxMissingMass);

		// Missing Mass vs Beam Energy: Found
		locHistName = "MissingMassVsBeamEnergy_Found";
		locHistTitle = locHistParticleName + string(" Found;Beam Energy (GeV); Missing Mass (GeV/c^{2})");
		dHist_MissingMassVsBeamEnergy_Found = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DBeamEnergyBins, dMinBeamEnergy, dMaxBeamEnergy, dNum2DMassBins, dMinMissingMass, dMaxMissingMass);

		// Missing Mass vs Beam Energy: Missing
		locHistName = "MissingMassVsBeamEnergy_Missing";
		locHistTitle = locHistParticleName + string(" Missing;Beam Energy (GeV); Missing Mass (GeV/c^{2})");
		dHist_MissingMassVsBeamEnergy_Missing = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DBeamEnergyBins, dMinBeamEnergy, dMaxBeamEnergy, dNum2DMassBins, dMinMissingMass, dMaxMissingMass);
	}
*/

	// Total Missing Energy vs Beam Energy
	locHistName = "TotalMissingEnergyVsBeamEnergy";
	locHistTitle = locHistParticleName + string(";Beam Energy (GeV); Total Missing Energy (Initial - Final) (GeV)");
	dHist_TotalMissingEnergyVsBeamEnergy = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DBeamEnergyBins, dMinBeamEnergy, dMaxBeamEnergy, dNum2DMissingBeamEBins, dMinMissingBeamE, dMaxMissingBeamE);

	// Total Missing Pt vs Beam Energy
	locHistName = "TotalMissingPtVsBeamEnergy";
	locHistTitle = locHistParticleName + string(";Beam Energy (GeV); Total Missing Transverse Momentum (GeV/c)");
	dHist_TotalMissingPtVsBeamEnergy = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DBeamEnergyBins, dMinBeamEnergy, dMaxBeamEnergy, dNum2DMissingPtBins, 0.0, dMaxMissingPt);

	// Total Missing Pt vs Beam Energy
	locHistName = "TotalMissingPyVsPx";
	locHistTitle = locHistParticleName + string("; Total Missing Px (GeV/c); Total Missing Py (GeV/c)");
	dHist_TotalMissingPyVsPx = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DMissingPxyBins, dMinMissingPxy, dMaxMissingPxy, dNum2DMissingPxyBins, dMinMissingPxy, dMaxMissingPxy);

	// Reconstructed Missing Mass Squared vs Beam Energy
	locHistName = "ReconMissingMassSquaredVsBeamEnergy";
	locHistTitle = locHistParticleName + string(" Reconstructed;Beam Energy (GeV); Missing Mass Squared (GeV/c^{2})^{2}");
	dHist_ReconMissingMassSquaredVsBeamEnergy = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DBeamEnergyBins, dMinBeamEnergy, dMaxBeamEnergy, 4000, -1.5, 0.5);

	gDirectory->cd("..");
}

void Selector_TrackEff::Create_MatchingHists(void)
{
	string locHistName, locHistTitle;
	string locHistParticleName = ParticleName_ROOT(dMissingPID);

	gDirectory->mkdir("Matching", "Matching");
	gDirectory->cd("Matching");

	//FOM
	locHistName = "MatchingFOM";
	locHistTitle = locHistParticleName + string(";Missing Match FOM");
	dHist_MatchingFOM = new TH1I(locHistName.c_str(), locHistTitle.c_str(), dNumMatchFOMBins, 0, 1.0);

	int locNumLogBins = 0;
	double* locLogBinning = Generate_LogYBinning(locNumLogBins);

	//FOM vs DeltaP
	locHistName = "MatchingFOMVsDeltaPOverP";
	locHistTitle = locHistParticleName + string(";#Deltap/p (Measured - Missing);Missing Match FOM");
	dHist_MatchingFOMVsDeltaPOverP = new TH2I(locHistName.c_str(), locHistTitle.c_str(), dNum2DDeltaPOverPBins, dMinDeltaPOverP, dMaxDeltaPOverP, dNum2DMatchFOMBins, 0, 1.0);
	locHistName = "MatchingFOMVsDeltaPOverP_LogY";
	dHist_MatchingFOMVsDeltaPOverP_LogY = new TH2I(locHistName.c_str(), locHistTitle.c_str(), dNum2DDeltaPOverPBins, dMinDeltaPOverP, dMaxDeltaPOverP, locNumLogBins, locLogBinning);

	//FOM vs DeltaTheta
	locHistName = "MatchingFOMVsDeltaTheta";
	locHistTitle = locHistParticleName + string(";#Delta#theta#circ;Missing Match FOM");
	dHist_MatchingFOMVsDeltaTheta = new TH2I(locHistName.c_str(), locHistTitle.c_str(), dNum2DDeltaThetaBins, dMinDeltaTheta, dMaxDeltaTheta, dNum2DMatchFOMBins, 0, 1.0);
	locHistName = "MatchingFOMVsDeltaTheta_LogY";
	dHist_MatchingFOMVsDeltaTheta_LogY = new TH2I(locHistName.c_str(), locHistTitle.c_str(), dNum2DDeltaThetaBins, dMinDeltaTheta, dMaxDeltaTheta, locNumLogBins, locLogBinning);

	//FOM vs DeltaPhi
	locHistName = "MatchingFOMVsDeltaPhi";
	locHistTitle = locHistParticleName + string(";#Delta#phi#circ;Missing Match FOM");
	dHist_MatchingFOMVsDeltaPhi = new TH2I(locHistName.c_str(), locHistTitle.c_str(), dNum2DDeltaPhiBins, dMinDeltaPhi, dMaxDeltaPhi, dNum2DMatchFOMBins, 0, 1.0);
	locHistName = "MatchingFOMVsDeltaPhi_LogY";
	dHist_MatchingFOMVsDeltaPhi_LogY = new TH2I(locHistName.c_str(), locHistTitle.c_str(), dNum2DDeltaPhiBins, dMinDeltaPhi, dMaxDeltaPhi, locNumLogBins, locLogBinning);

	gDirectory->cd("..");
}

double* Selector_TrackEff::Generate_LogYBinning(int& locNumBins)
{
	int locLowest10Power = -40, locHighest10Power = 0;
	unsigned int locNumBinsPerPower = 20;

	int locNumPowerRanges = locHighest10Power - locLowest10Power; //Num powers of 10, minus 1
	locNumBins = locNumBinsPerPower*locNumPowerRanges;

	double* locBinArray = new double[locNumBins + 1];
	locBinArray[0] = pow(10.0, locLowest10Power);
	for(int loc_j = 0; loc_j < locNumPowerRanges; ++loc_j)
	{
		double locCurrent10Power = double(locLowest10Power + loc_j);
		for(unsigned int loc_k = 0; loc_k < locNumBinsPerPower; ++loc_k)
		{
		        double locMultiplier = (double(loc_k) + 1.0) / locNumBinsPerPower;
			locBinArray[loc_j*locNumBinsPerPower + loc_k + 1] = pow(10,locMultiplier) * pow(10.0, locCurrent10Power);
		}
	}

	return locBinArray;
}

void Selector_TrackEff::Create_ResolutionHists(void)
{
	string locHistName, locHistTitle;
	string locHistParticleName = ParticleName_ROOT(dMissingPID);

	//Dir Name
	string locDirectoryName = "Resolution";

	//Make dir
	gDirectory->mkdir(locDirectoryName.c_str(), locDirectoryName.c_str());
	gDirectory->cd(locDirectoryName.c_str());

	for(int loc_i = 0; loc_i < dNumFitBeamEnergyBins; ++loc_i)
	{
		ostringstream locDirNameStream;
		locDirNameStream << "EnergyBin_" << loc_i;
		gDirectory->mkdir(locDirNameStream.str().c_str(), locDirNameStream.str().c_str());
		gDirectory->cd(locDirNameStream.str().c_str());

		// DeltaP Vs P
		locHistName = "DeltaPVsP";
		locHistTitle = locHistParticleName + string(";Missing p (GeV/c);#Deltap (GeV/c) (Measured - Missing)");
		dHistMap_Resolution_DeltaPVsP[loc_i] = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DPBins, dMinP, dMaxP, dNum2DDeltaPBins, dMinDeltaP, dMaxDeltaP);

		// DeltaP Vs Theta
		locHistName = string("DeltaPVsTheta");
		locHistTitle = locHistParticleName + string(";Missing #theta#circ;#Deltap (GeV/c) (Measured - Missing)");
		dHistMap_Resolution_DeltaPVsTheta[loc_i] = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DDeltaPBins, dMinDeltaP, dMaxDeltaP);

		// DeltaP Vs Phi
		locHistName = string("DeltaPVsPhi");
		locHistTitle = locHistParticleName + string(";Missing #phi#circ;#Deltap (GeV/c) (Measured - Missing)");
		dHistMap_Resolution_DeltaPVsPhi[loc_i] = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DDeltaPBins, dMinDeltaP, dMaxDeltaP);

		// DeltaP/P Vs P
		locHistName = "DeltaPOverPVsP";
		locHistTitle = locHistParticleName + string(";Missing p (GeV/c);#Deltap/p (Measured - Missing)");
		dHistMap_Resolution_DeltaPOverPVsP[loc_i] = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DPBins, dMinP, dMaxP, dNum2DDeltaPOverPBins, dMinDeltaPOverP, dMaxDeltaPOverP);

		// DeltaP/P Vs Theta
		locHistName = string("DeltaPOverPVsTheta");
		locHistTitle = locHistParticleName + string(";Missing #theta#circ;#Deltap/p (Measured - Missing)");
		dHistMap_Resolution_DeltaPOverPVsTheta[loc_i] = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DDeltaPOverPBins, dMinDeltaPOverP, dMaxDeltaPOverP);

		// DeltaP/P Vs Phi
		locHistName = string("DeltaPOverPVsPhi");
		locHistTitle = locHistParticleName + string(";Missing #phi#circ;#Deltap/p (Measured - Missing)");
		dHistMap_Resolution_DeltaPOverPVsPhi[loc_i] = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DDeltaPOverPBins, dMinDeltaPOverP, dMaxDeltaPOverP);

		// DeltaTheta Vs P
		locHistName = string("DeltaThetaVsP");
		locHistTitle = locHistParticleName + string(";Missing p (GeV/c);#Delta#theta#circ (Measured - Missing)");
		dHistMap_Resolution_DeltaThetaVsP[loc_i] = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DPBins, dMinP, dMaxP, dNum2DDeltaThetaBins, dMinDeltaTheta, dMaxDeltaTheta);

		// DeltaTheta Vs Theta
		locHistName = string("DeltaThetaVsTheta");
		locHistTitle = locHistParticleName + string(";Missing #theta#circ;#Delta#theta#circ (Measured - Missing)");
		dHistMap_Resolution_DeltaThetaVsTheta[loc_i] = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DDeltaThetaBins, dMinDeltaTheta, dMaxDeltaTheta);

		// DeltaTheta Vs Phi
		locHistName = string("DeltaThetaVsPhi");
		locHistTitle = locHistParticleName + string(";Missing #phi#circ;#Delta#theta#circ (Measured - Missing)");
		dHistMap_Resolution_DeltaThetaVsPhi[loc_i] = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DDeltaThetaBins, dMinDeltaTheta, dMaxDeltaTheta);

		// DeltaPhi Vs P
		locHistName = string("DeltaPhiVsP");
		locHistTitle = locHistParticleName + string(";Missing p (GeV/c);#Delta#phi#circ (Measured - Missing)");
		dHistMap_Resolution_DeltaPhiVsP[loc_i] = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DPBins, dMinP, dMaxP, dNum2DDeltaPhiBins, dMinDeltaPhi, dMaxDeltaPhi);

		// DeltaPhi Vs Theta
		locHistName = string("DeltaPhiVsTheta");
		locHistTitle = locHistParticleName + string(";Missing #theta#circ;#Delta#phi#circ (Measured - Missing)");
		dHistMap_Resolution_DeltaPhiVsTheta[loc_i] = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DDeltaPhiBins, dMinDeltaPhi, dMaxDeltaPhi);

		// DeltaPhi Vs Theta
		locHistName = string("DeltaPhiVsPhi");
		locHistTitle = locHistParticleName + string(";Missing #phi#circ;#Delta#phi#circ (Measured - Missing)");
		dHistMap_Resolution_DeltaPhiVsPhi[loc_i] = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DDeltaPhiBins, dMinDeltaPhi, dMaxDeltaPhi);

		gDirectory->cd("..");
	}
	gDirectory->cd("..");
}

void Selector_TrackEff::Create_MissingP3SigmaHists(void)
{
	string locHistName, locHistTitle;
	string locHistParticleName = ParticleName_ROOT(dMissingPID);

	//Dir Name
	string locDirectoryName = "MissingP3Sigma";

	//Make dir
	gDirectory->mkdir(locDirectoryName.c_str(), locDirectoryName.c_str());
	gDirectory->cd(locDirectoryName.c_str());

	for(int loc_i = 0; loc_i < dNumFitBeamEnergyBins; ++loc_i)
	{
		ostringstream locDirNameStream;
		locDirNameStream << "EnergyBin_" << loc_i;
		gDirectory->mkdir(locDirNameStream.str().c_str(), locDirNameStream.str().c_str());
		gDirectory->cd(locDirNameStream.str().c_str());

		// MissingPSigmaVsP
		locHistName = "MissingPSigmaVsP";
		locHistTitle = locHistParticleName + string(";Missing p (GeV/c);Missing p #sigma (GeV/c)");
		dHistMap_MissingPSigmaVsP[loc_i] = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DPBins, dMinP, dMaxP, dNum2DMissingPSigmaBins, 0.0, dMaxMissingPSigma);

		// MissingPSigmaVsTheta
		locHistName = string("MissingPSigmaVsTheta");
		locHistTitle = locHistParticleName + string(";Missing #theta#circ;Missing p #sigma (GeV/c)");
		dHistMap_MissingPSigmaVsTheta[loc_i] = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DMissingPSigmaBins, 0.0, dMaxMissingPSigma);

		// MissingPSigmaVsPhi
		locHistName = string("MissingPSigmaVsPhi");
		locHistTitle = locHistParticleName + string(";Missing #phi#circ;Missing p #sigma (GeV/c)");
		dHistMap_MissingPSigmaVsPhi[loc_i] = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DMissingPSigmaBins, 0.0, dMaxMissingPSigma);

		// MissingPSigmaOverPVsP
		locHistName = "MissingPSigmaOverPVsP";
		locHistTitle = locHistParticleName + string(";Missing p (GeV/c);Missing p #sigma/p");
		dHistMap_MissingPSigmaOverPVsP[loc_i] = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DPBins, dMinP, dMaxP, dNum2DMissingPSigmaOverPBins, 0.0, dMaxMissingPSigmaOverP);

		// MissingPSigmaOverPVsTheta
		locHistName = string("MissingPSigmaOverPVsTheta");
		locHistTitle = locHistParticleName + string(";Missing #theta#circ;Missing p #sigma/p");
		dHistMap_MissingPSigmaOverPVsTheta[loc_i] = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DMissingPSigmaOverPBins, 0.0, dMaxMissingPSigmaOverP);

		// MissingPSigmaOverPVsPhi
		locHistName = string("MissingPSigmaOverPVsPhi");
		locHistTitle = locHistParticleName + string(";Missing #phi#circ;Missing p #sigma/p");
		dHistMap_MissingPSigmaOverPVsPhi[loc_i] = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DMissingPSigmaOverPBins, 0.0, dMaxMissingPSigmaOverP);

		// MissingThetaSigmaVsP
		locHistName = "MissingThetaSigmaVsP";
		locHistTitle = locHistParticleName + string(";Missing p (GeV/c);Missing #theta #sigma (#circ)");
		dHistMap_MissingThetaSigmaVsP[loc_i] = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DPBins, dMinP, dMaxP, dNum2DMissingThetaSigmaBins, 0.0, dMaxMissingThetaSigma);

		// MissingThetaSigmaVsTheta
		locHistName = string("MissingThetaSigmaVsTheta");
		locHistTitle = locHistParticleName + string(";Missing #theta#circ;Missing #theta #sigma (#circ)");
		dHistMap_MissingThetaSigmaVsTheta[loc_i] = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DMissingThetaSigmaBins, 0.0, dMaxMissingThetaSigma);

		// MissingThetaSigmaVsPhi
		locHistName = string("MissingThetaSigmaVsPhi");
		locHistTitle = locHistParticleName + string(";Missing #phi#circ;Missing #theta #sigma (#circ)");
		dHistMap_MissingThetaSigmaVsPhi[loc_i] = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DMissingThetaSigmaBins, 0.0, dMaxMissingThetaSigma);

		// MissingPhiSigmaVsP
		locHistName = "MissingPhiSigmaVsP";
		locHistTitle = locHistParticleName + string(";Missing p (GeV/c);Missing #phi #sigma (#circ)");
		dHistMap_MissingPhiSigmaVsP[loc_i] = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DPBins, dMinP, dMaxP, dNum2DMissingPhiSigmaBins, 0.0, dMaxMissingPhiSigma);

		// MissingPhiSigmaVsTheta
		locHistName = string("MissingPhiSigmaVsTheta");
		locHistTitle = locHistParticleName + string(";Missing #theta#circ;Missing #phi #sigma (#circ)");
		dHistMap_MissingPhiSigmaVsTheta[loc_i] = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DMissingPhiSigmaBins, 0.0, dMaxMissingPhiSigma);

		// MissingPhiSigmaVsPhi
		locHistName = string("MissingPhiSigmaVsPhi");
		locHistTitle = locHistParticleName + string(";Missing #phi#circ;Missing #phi #sigma (#circ)");
		dHistMap_MissingPhiSigmaVsPhi[loc_i] = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DMissingPhiSigmaBins, 0.0, dMaxMissingPhiSigma);

		gDirectory->cd("..");
	}
	gDirectory->cd("..");
}

void Selector_TrackEff::Create_EfficiencyHists(bool locSignalPeakFlag)
{
	string locHistName, locHistTitle;
	string locHistParticleName = ParticleName_ROOT(dMissingPID);

	//Dir Name
	string locDirectoryName = "Efficiency";
	if(locSignalPeakFlag)
		locDirectoryName += "_SignalPeak";
	else
		locDirectoryName += "_SignalSideBands";

	gDirectory->mkdir(locDirectoryName.c_str(), locDirectoryName.c_str());
	gDirectory->cd(locDirectoryName.c_str());

	for(int loc_i = 0; loc_i < dNumFitBeamEnergyBins; ++loc_i)
	{
		ostringstream locDirNameStream;
		locDirNameStream << "EnergyBin_" << loc_i;
		gDirectory->mkdir(locDirNameStream.str().c_str(), locDirNameStream.str().c_str());
		gDirectory->cd(locDirNameStream.str().c_str());

		//Resize vectors
		dHistMap_TrackFound_PVsTheta[locSignalPeakFlag][loc_i].resize(dNumVertexZBins);
		dHistMap_TrackMissing_PVsTheta[locSignalPeakFlag][loc_i].resize(dNumVertexZBins);

		dHistMap_TrackFound_PhiVsTheta[locSignalPeakFlag][loc_i].resize(dNumVertexZBins);
		dHistMap_TrackMissing_PhiVsTheta[locSignalPeakFlag][loc_i].resize(dNumVertexZBins);

		for(int locVertexZBin = 0; locVertexZBin < dNumVertexZBins; ++locVertexZBin)
		{
			double locBinMinVertZ = dMinVertexZ + double(locVertexZBin)*dVertexZBinSize;
			double locBinMaxVertZ = locBinMinVertZ + dVertexZBinSize;

			ostringstream locVertexZBinStream, locVertexZRangeStream;
			locVertexZBinStream << locVertexZBin;
			locVertexZRangeStream << locBinMinVertZ << " #leq Vertex-Z (cm) < " << locBinMaxVertZ;

			//Reconstruction
			locHistName = string("PVsTheta_TrackFound_VertexZBin") + locVertexZBinStream.str();
			locHistTitle = locHistParticleName + string(", ") + locVertexZRangeStream.str() + string(", Track Found;Missing #theta#circ;Missing p (GeV/c)");
			dHistMap_TrackFound_PVsTheta[locSignalPeakFlag][loc_i][locVertexZBin] = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);

			locHistName = string("PVsTheta_TrackMissing_VertexZBin") + locVertexZBinStream.str();
			locHistTitle = locHistParticleName + string(", ") + locVertexZRangeStream.str() + string(", Track Missing;Missing #theta#circ;Missing p (GeV/c)");
			dHistMap_TrackMissing_PVsTheta[locSignalPeakFlag][loc_i][locVertexZBin] = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);

			locHistName = string("PhiVsTheta_TrackFound_VertexZBin") + locVertexZBinStream.str();
			locHistTitle = locHistParticleName + string(", ") + locVertexZRangeStream.str() + string(", Track Found;Missing #theta#circ;Missing #phi#circ");
			dHistMap_TrackFound_PhiVsTheta[locSignalPeakFlag][loc_i][locVertexZBin] = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPhiBins, dMinPhi, dMaxPhi);

			locHistName = string("PhiVsTheta_TrackMissing_VertexZBin") + locVertexZBinStream.str();
			locHistTitle = locHistParticleName + string(", ") + locVertexZRangeStream.str() + string(", Track Missing;Missing #theta#circ;Missing #phi#circ");
			dHistMap_TrackMissing_PhiVsTheta[locSignalPeakFlag][loc_i][locVertexZBin] = new TH2D(locHistName.c_str(), locHistTitle.c_str(), dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPhiBins, dMinPhi, dMaxPhi);
		}

		gDirectory->cd("..");
	}

	gDirectory->cd("..");
}

/*********************************************************************** SETUP CUTS ***********************************************************************/

void Selector_TrackEff::Setup_Cuts(void)
{
	return;
	Setup_ResolutionCuts();
	Setup_MassCuts();
}

void Selector_TrackEff::Setup_ResolutionCuts(void)
{
	//resolution file
	string locResolutionFileName = string("find_res_merged/trackeff_") + string(ParticleType(dMissingPID));
	if(dMissingPID != Proton)
	{
		if(dIs4PiFlag)
			locResolutionFileName += "_4pi.root";
		else
			locResolutionFileName += "_2pi.root";
	}
	else
		locResolutionFileName += ".root";
	locResolutionFileName += ".res_subtracted.root.fit.root.smoothed.root";

	vector<string> locHistTypes = {"DeltaPOverPVsP", "DeltaPOverPVsTheta", "DeltaPOverPVsPhi", "DeltaThetaVsP", "DeltaThetaVsTheta",
			"DeltaThetaVsPhi", "DeltaPhiVsP", "DeltaPhiVsTheta", "DeltaPhiVsPhi"};

	TFile* locResolutionFile = new TFile(locResolutionFileName.c_str(), "READ");
	if(!locResolutionFile->IsZombie())
	{
		for(int loc_i = 0; loc_i < dNumFitBeamEnergyBins; ++loc_i)
		{
			ostringstream locDirNameStream;
			locDirNameStream << "Resolution/EnergyBin_" << loc_i;
			for(auto locHistName : locHistTypes)
			{
				//means
				string locGraphName = locHistName + string("_Means");
				string locGraphPath = locDirNameStream.str() + string("/") + locGraphName;
				TGraphErrors* locGraphErrors = (TGraphErrors*)locResolutionFile->Get(locGraphPath.c_str());
				if(locGraphErrors != nullptr)
				{
					string locFuncName = string(locGraphErrors->GetName()) + string("_Func");
					TF1* locFunc = (TF1*)locGraphErrors->GetListOfFunctions()->FindObject(locFuncName.c_str());
					if(locFunc != nullptr)
						dResolutionFuncMap_Means[locHistName][loc_i] = locFunc;
				}

				//sigmas
				locGraphName = locHistName + string("_Sigmas");
				locGraphPath = locDirNameStream.str() + string("/") + locGraphName;
				locGraphErrors = (TGraphErrors*)locResolutionFile->Get(locGraphPath.c_str());
				if(locGraphErrors != nullptr)
				{
					string locFuncName = string(locGraphErrors->GetName()) + string("_Func");
					TF1* locFunc = (TF1*)locGraphErrors->GetListOfFunctions()->FindObject(locFuncName.c_str());
					if(locFunc != nullptr)
						dResolutionFuncMap_Sigmas[locHistName][loc_i] = locFunc;
				}
			}
		}
		locResolutionFile->Close();
		delete locResolutionFile;
	}

	//Fill in gaps if fits missing
	for(auto locHistName : locHistTypes)
	{
		//overall
		if((dResolutionFuncMap_Means.find(locHistName) == dResolutionFuncMap_Means.end()) || (dResolutionFuncMap_Sigmas.find(locHistName) == dResolutionFuncMap_Sigmas.end()))
		{
			cout << "assign all res cuts NULL: " << locHistName << endl;
			for(int loc_i = 0; loc_i < dNumFitBeamEnergyBins; ++loc_i)
			{
				dResolutionFuncMap_Means[locHistName][loc_i] = nullptr;
				dResolutionFuncMap_Sigmas[locHistName][loc_i] = nullptr;
			}
			continue;
		}

		//by bin
		for(int loc_i = 0; loc_i < dNumFitBeamEnergyBins; ++loc_i)
		{
			//means
			if(dResolutionFuncMap_Means[locHistName].find(loc_i) == dResolutionFuncMap_Means[locHistName].end())
			{
				//search for nearest non-null
				int locBestDistance = 9999;
				TF1* locBestFunc = nullptr;
				for(auto& locMapPair : dResolutionFuncMap_Means[locHistName])
				{
					int locDistance = locMapPair.first - loc_i;
					if(locDistance > locBestDistance)
						continue;
					locBestDistance = locDistance;
					locBestFunc = locMapPair.second;
				}
				dResolutionFuncMap_Means[locHistName][loc_i] = locBestFunc;
				cout << "set mean " << locHistName << loc_i << " func to: " << locBestFunc << endl;
			}

			//sigmas
			if(dResolutionFuncMap_Sigmas[locHistName].find(loc_i) == dResolutionFuncMap_Sigmas[locHistName].end())
			{
				//search for nearest non-null
				int locBestDistance = 9999;
				TF1* locBestFunc = nullptr;
				for(auto& locMapPair : dResolutionFuncMap_Sigmas[locHistName])
				{
					int locDistance = locMapPair.first - loc_i;
					if(locDistance > locBestDistance)
						continue;
					locBestDistance = locDistance;
					locBestFunc = locMapPair.second;
				}
				dResolutionFuncMap_Sigmas[locHistName][loc_i] = locBestFunc;
				cout << "set sigma " << locHistName << loc_i << " func to: " << locBestFunc << endl;
			}
		}
	}
}

void Selector_TrackEff::Setup_MassCuts(void)
{
	//resolution file
	string locMassFileName = string("mass_cuts/trackeff_") + string(ParticleType(dMissingPID));
	if(dMissingPID != Proton)
	{
		if(dIs4PiFlag)
			locMassFileName += "_4pi.root";
		else
			locMassFileName += "_2pi.root";
	}
	else
		locMassFileName += ".root";
	locMassFileName += ".mass_fit.root.smoothed.root";

	TFile* locMassFile = new TFile(locMassFileName.c_str(), "READ");
	if(locMassFile->IsZombie())
	{
		delete locMassFile;
		dMassFuncMap_Means[true] = nullptr;
		dMassFuncMap_Sigmas[true] = nullptr;
		dMassFuncMap_Means[false] = nullptr;
		dMassFuncMap_Sigmas[false] = nullptr;
		return;
	}

	vector<string> locHistTypes = {"MissingMassVsBeamEnergy_Found", "MissingMassVsBeamEnergy_Missing"};
	for(size_t loc_i = 0; loc_i < locHistTypes.size(); ++loc_i)
	{
		string locHistName = locHistTypes[loc_i];
		bool locFoundFlag = (loc_i == 0);

		//means
		string locGraphName = locHistName + string("_Means");
		string locGraphPath = string("MissingMass/") + locGraphName;
		TGraphErrors* locGraphErrors = (TGraphErrors*)locMassFile->Get(locGraphPath.c_str());
		if(locGraphErrors != nullptr)
		{
			string locFuncName = string(locGraphErrors->GetName()) + string("_Func");
			dMassFuncMap_Means[locFoundFlag] = (TF1*)locGraphErrors->GetListOfFunctions()->FindObject(locFuncName.c_str());
		}
		else
			dMassFuncMap_Means[locFoundFlag] = nullptr;

		//sigmas
		locGraphName = locHistName + string("_Sigmas");
		locGraphPath = string("MissingMass/") + locGraphName;
		locGraphErrors = (TGraphErrors*)locMassFile->Get(locGraphPath.c_str());
		if(locGraphErrors != nullptr)
		{
			string locFuncName = string(locGraphErrors->GetName()) + string("_Func");
			dMassFuncMap_Sigmas[locFoundFlag] = (TF1*)locGraphErrors->GetListOfFunctions()->FindObject(locFuncName.c_str());
		}
		else
			dMassFuncMap_Sigmas[locFoundFlag] = nullptr;
	}

	locMassFile->Close();
}

/********************************************************************* PROCESS EVENTS *********************************************************************/

Bool_t Selector_TrackEff::Process(Long64_t entry)
{
	// The Process() function is called for each entry in the tree (or possibly
	// keyed object in the case of PROOF) to be processed. The entry argument
	// specifies which entry in the currently loaded tree is to be processed.
	// It can be passed to either Selector_TrackEff::GetEntry() or TBranch::GetEntry()
	// to read either all or the required parts of the data. When processing
	// keyed objects with PROOF, the object is already loaded and is available
	// via the fObject pointer.
	//
	// This function should contain the "body" of the analysis. It can contain
	// simple or elaborate selection criteria, run algorithms on the data
	// of the event and typically fill histograms.
	//
	// The processing can be stopped by calling Abort().
	//
	// Use fStatus to set the return value of TTree::Process().
	//
	// The return value is currently not used.

	++dComboCounter;
	if(dComboCounter % 100000 == 0)
		cout << "Event: " << dComboCounter << endl;

	//Get the data for this combo
	GetEntry(entry);

	//Require that the fit converged ... is sort of a resolution cut
	if(KinFitNDF == 0)
		return kTRUE;

	//Find if in RF signal or sideband
	bool locRFCenterFlag = fabs(BeamRFDeltaT) > dBeamBunchPeriod;
	double locRFWeight = Get_RFWeight();

	//Fill matching hists
	if(locRFCenterFlag)
		Fill_MatchingHists();

	//Get beam energy bin
	int locBeamEnergyBin = Select_BeamEnergyBin();
	if(locBeamEnergyBin == -1)
		return kTRUE;

	Int_t locBestTrackMatchIndex = Fill_ResolutionHists(locRFWeight, locBeamEnergyBin);
	bool locTrackFoundFlag = (locBestTrackMatchIndex != -1);
//	if(locTrackFoundFlag)
		Fill_MissingP3SigmaHists(locRFWeight, locBeamEnergyBin);
	Fill_MassHists(locRFWeight, locBestTrackMatchIndex);

	bool locSignalPeakFlag = Select_MassBin(locTrackFoundFlag);
	Fill_EffHists(locRFWeight, locTrackFoundFlag, locSignalPeakFlag, locBeamEnergyBin);

	return kTRUE;
}

double Selector_TrackEff::Get_RFWeight(void)
{
	if(std::isnan(BeamRFDeltaT))
		return 0.0;

	//check if signal peak
	double locMaxDeltaT = dBeamBunchPeriod*(double(dBeamBunchRange_RFSignal.first) + 0.5);
	if(fabs(BeamRFDeltaT) <= locMaxDeltaT)
		return 1.0;

	//check if sideband peak
	double locMinDeltaT = dBeamBunchPeriod*(double(dBeamBunchRange_RFSideband.first) - 0.5);
	locMaxDeltaT = dBeamBunchPeriod*(double(dBeamBunchRange_RFSideband.first) + 0.5);
	if((fabs(BeamRFDeltaT) < locMinDeltaT) || (fabs(BeamRFDeltaT) > locMaxDeltaT))
		return dRFSidebandWeight;

	//not used
	return 0.0;
}

bool Selector_TrackEff::Select_MassBin(bool locTrackFoundFlag)
{
	if(dIs3PiFlag)
		return ((MissingMassSquared >= -0.1) && (MissingMassSquared <= 0.15));

	//true is signal peak region, false is mass sidebands
	if((dMassFuncMap_Means[locTrackFoundFlag] == nullptr) || (dMassFuncMap_Sigmas[locTrackFoundFlag] == nullptr))
		return true;

	double locMassValue = (dMissingPID == Proton) ? sqrt(MissingMassSquared) : MissingMassSquared;
	double locMean = dMassFuncMap_Means[locTrackFoundFlag]->Eval(BeamEnergy);
	double loc3Sigma = 3.0*dMassFuncMap_Sigmas[locTrackFoundFlag]->Eval(BeamEnergy);

	return ((locMassValue >= (locMean - loc3Sigma)) && (locMassValue <= (locMean + loc3Sigma)));
}

int Selector_TrackEff::Select_BeamEnergyBin(void)
{
	int locBeamEnergyBin = (BeamEnergy - dMinBeamEnergy)/dFitBeamEnergyBinSize;
	if(locBeamEnergyBin < 0)
		return -1;
	if(locBeamEnergyBin >= int(dNumFitBeamEnergyBins))
		return -1;
	return locBeamEnergyBin;
}

bool Selector_TrackEff::Cut_DeltaPOverP(size_t locEBin, double locP, double locTheta, double locPhi, double locDeltaPOverP)
{
	return (fabs(locDeltaPOverP) <= 0.6);

	//input & func angles in degrees
	//only reject if it fails ALL cuts
		//e.g. if in area of poor resolution (theta = 15), may fail the vs phi cut (averaged over all angles)
	if((dResolutionFuncMap_Means["DeltaPOverPVsP"][locEBin] == nullptr) || (dResolutionFuncMap_Sigmas["DeltaPOverPVsP"][locEBin] == nullptr))
	{
		cout << "p missing, beamE, ebin = " << BeamEnergy << ", " << locEBin << endl;
		return true;
	}
	if((dResolutionFuncMap_Means["DeltaPOverPVsTheta"][locEBin] == nullptr) || (dResolutionFuncMap_Sigmas["DeltaPOverPVsTheta"][locEBin] == nullptr))
	{
		cout << "theta missing, beamE, ebin = " << BeamEnergy << ", " << locEBin << endl;
		return true;
	}
	if((dResolutionFuncMap_Means["DeltaPOverPVsPhi"][locEBin] == nullptr) || (dResolutionFuncMap_Sigmas["DeltaPOverPVsPhi"][locEBin] == nullptr))
	{
		cout << "phi missing, beamE, ebin = " << BeamEnergy << ", " << locEBin << endl;
		return true;
	}

	//vs p
	TF1* locMeanFunc = dResolutionFuncMap_Means["DeltaPOverPVsP"][locEBin];
	TF1* locSigmaFunc = dResolutionFuncMap_Sigmas["DeltaPOverPVsP"][locEBin];
	if((locP >= locMeanFunc->GetXmin()) && (locP <= locMeanFunc->GetXmax()) && (locP >= locSigmaFunc->GetXmin()) && (locP <= locSigmaFunc->GetXmax()))
	{
		double locMean = locMeanFunc->Eval(locP);
		double loc3Sigma = 3.0*locSigmaFunc->Eval(locP);
//		if(gDebugFlag && (fabs(locDeltaPOverP) > 0.8))
//			cout << "p, ebin, vsp mean, vsp 3sigma, delta = " << locP << ", " << locEBin << ", " << locMean << ", " << loc3Sigma << ", " << locDeltaPOverP << endl;
		if((locDeltaPOverP >= (locMean - loc3Sigma)) && (locDeltaPOverP <= (locMean + loc3Sigma)))
			return true;
	}

	//vs theta
	locMeanFunc = dResolutionFuncMap_Means["DeltaPOverPVsTheta"][locEBin];
	locSigmaFunc = dResolutionFuncMap_Sigmas["DeltaPOverPVsTheta"][locEBin];
	if((locTheta >= locMeanFunc->GetXmin()) && (locTheta <= locMeanFunc->GetXmax()) && (locTheta >= locSigmaFunc->GetXmin()) && (locTheta <= locSigmaFunc->GetXmax()))
	{
		double locMean = locMeanFunc->Eval(locTheta);
		double loc3Sigma = 3.0*locSigmaFunc->Eval(locTheta);
//		if(gDebugFlag && (fabs(locDeltaPOverP) > 0.8))
//			cout << "theta, ebin, vs theta mean, vs theta 3sigma, delta = " << locTheta << ", " << locEBin << ", " << locMean << ", " << loc3Sigma << ", " << locDeltaPOverP << endl;
		if((locDeltaPOverP >= (locMean - loc3Sigma)) && (locDeltaPOverP <= (locMean + loc3Sigma)))
			return true;
	}

	//vs phi
	locMeanFunc = dResolutionFuncMap_Means["DeltaPOverPVsPhi"][locEBin];
	locSigmaFunc = dResolutionFuncMap_Sigmas["DeltaPOverPVsPhi"][locEBin];
	if((locPhi >= locMeanFunc->GetXmin()) && (locPhi <= locMeanFunc->GetXmax()) && (locPhi >= locSigmaFunc->GetXmin()) && (locPhi <= locSigmaFunc->GetXmax()))
	{
		double locMean = locMeanFunc->Eval(locPhi);
		double loc3Sigma = 3.0*locSigmaFunc->Eval(locPhi);
//		if(gDebugFlag && (fabs(locDeltaPOverP) > 0.8))
//			cout << "phi, ebin, vs phi mean, vs phi 3sigma, delta = " << locPhi << ", " << locEBin << ", " << locMean << ", " << loc3Sigma << ", " << locDeltaPOverP << endl;
		if((locDeltaPOverP >= (locMean - loc3Sigma)) && (locDeltaPOverP <= (locMean + loc3Sigma)))
			return true;
	}

	return false;
}

bool Selector_TrackEff::Cut_DeltaTheta(size_t locEBin, double locP, double locTheta, double locPhi, double locDeltaTheta)
{
	return (fabs(locDeltaTheta) <= 30.0);

	//input & func angles in degrees
	//only reject if it fails ALL cuts
		//e.g. if in area of poor resolution (theta = 15), may fail the vs phi cut (averaged over all angles)

	if((dResolutionFuncMap_Means["DeltaThetaVsP"][locEBin] == nullptr) || (dResolutionFuncMap_Sigmas["DeltaThetaVsP"][locEBin] == nullptr))
		return true;
	if((dResolutionFuncMap_Means["DeltaThetaVsTheta"][locEBin] == nullptr) || (dResolutionFuncMap_Sigmas["DeltaThetaVsTheta"][locEBin] == nullptr))
		return true;
	if((dResolutionFuncMap_Means["DeltaThetaVsPhi"][locEBin] == nullptr) || (dResolutionFuncMap_Sigmas["DeltaThetaVsPhi"][locEBin] == nullptr))
		return true;

	//vs p
	TF1* locMeanFunc = dResolutionFuncMap_Means["DeltaThetaVsP"][locEBin];
	TF1* locSigmaFunc = dResolutionFuncMap_Sigmas["DeltaThetaVsP"][locEBin];
	if((locP >= locMeanFunc->GetXmin()) && (locP <= locMeanFunc->GetXmax()) && (locP >= locSigmaFunc->GetXmin()) && (locP <= locSigmaFunc->GetXmax()))
	{
		double locMean = locMeanFunc->Eval(locP);
		double loc3Sigma = 3.0*locSigmaFunc->Eval(locP);
		if((locDeltaTheta >= (locMean - loc3Sigma)) && (locDeltaTheta <= (locMean + loc3Sigma)))
			return true;
	}

	//vs theta
	locMeanFunc = dResolutionFuncMap_Means["DeltaThetaVsTheta"][locEBin];
	locSigmaFunc = dResolutionFuncMap_Sigmas["DeltaThetaVsTheta"][locEBin];
	if((locTheta >= locMeanFunc->GetXmin()) && (locTheta <= locMeanFunc->GetXmax()) && (locTheta >= locSigmaFunc->GetXmin()) && (locTheta <= locSigmaFunc->GetXmax()))
	{
		double locMean = locMeanFunc->Eval(locTheta);
		double loc3Sigma = 3.0*locSigmaFunc->Eval(locTheta);
		if((locDeltaTheta >= (locMean - loc3Sigma)) && (locDeltaTheta <= (locMean + loc3Sigma)))
			return true;
	}

	//vs phi
	locMeanFunc = dResolutionFuncMap_Means["DeltaThetaVsPhi"][locEBin];
	locSigmaFunc = dResolutionFuncMap_Sigmas["DeltaThetaVsPhi"][locEBin];
	if((locPhi >= locMeanFunc->GetXmin()) && (locPhi <= locMeanFunc->GetXmax()) && (locPhi >= locSigmaFunc->GetXmin()) && (locPhi <= locSigmaFunc->GetXmax()))
	{
		double locMean = locMeanFunc->Eval(locPhi);
		double loc3Sigma = 3.0*locSigmaFunc->Eval(locPhi);
		if((locDeltaTheta >= (locMean - loc3Sigma)) && (locDeltaTheta <= (locMean + loc3Sigma)))
			return true;
	}

	return false;
}

bool Selector_TrackEff::Cut_DeltaPhi(size_t locEBin, double locP, double locTheta, double locPhi, double locDeltaPhi)
{
	if(locTheta < 5.0)
		return true;
	return (fabs(locDeltaPhi) <= 30.0);

	//input & func angles in degrees
	//only reject if it fails ALL cuts
		//e.g. if in area of poor resolution (theta = 15), may fail the vs phi cut (averaged over all angles)

	if((dResolutionFuncMap_Means["DeltaPhiVsP"][locEBin] == nullptr) || (dResolutionFuncMap_Sigmas["DeltaPhiVsP"][locEBin] == nullptr))
		return true;
	if((dResolutionFuncMap_Means["DeltaPhiVsTheta"][locEBin] == nullptr) || (dResolutionFuncMap_Sigmas["DeltaPhiVsTheta"][locEBin] == nullptr))
		return true;
	if((dResolutionFuncMap_Means["DeltaPhiVsPhi"][locEBin] == nullptr) || (dResolutionFuncMap_Sigmas["DeltaPhiVsPhi"][locEBin] == nullptr))
		return true;

	//vs p
	TF1* locMeanFunc = dResolutionFuncMap_Means["DeltaPhiVsP"][locEBin];
	TF1* locSigmaFunc = dResolutionFuncMap_Sigmas["DeltaPhiVsP"][locEBin];
	if((locP >= locMeanFunc->GetXmin()) && (locP <= locMeanFunc->GetXmax()) && (locP >= locSigmaFunc->GetXmin()) && (locP <= locSigmaFunc->GetXmax()))
	{
		double locMean = locMeanFunc->Eval(locP);
		double loc3Sigma = 3.0*locSigmaFunc->Eval(locP);
		if((locDeltaPhi >= (locMean - loc3Sigma)) && (locDeltaPhi <= (locMean + loc3Sigma)))
			return true;
	}

	//vs theta
	locMeanFunc = dResolutionFuncMap_Means["DeltaPhiVsTheta"][locEBin];
	locSigmaFunc = dResolutionFuncMap_Sigmas["DeltaPhiVsTheta"][locEBin];
	if((locTheta >= locMeanFunc->GetXmin()) && (locTheta <= locMeanFunc->GetXmax()) && (locTheta >= locSigmaFunc->GetXmin()) && (locTheta <= locSigmaFunc->GetXmax()))
	{
		double locMean = locMeanFunc->Eval(locTheta);
		double loc3Sigma = 3.0*locSigmaFunc->Eval(locTheta);
		if((locDeltaPhi >= (locMean - loc3Sigma)) && (locDeltaPhi <= (locMean + loc3Sigma)))
			return true;
	}

	//vs phi
	locMeanFunc = dResolutionFuncMap_Means["DeltaPhiVsPhi"][locEBin];
	locSigmaFunc = dResolutionFuncMap_Sigmas["DeltaPhiVsPhi"][locEBin];
	if((locPhi >= locMeanFunc->GetXmin()) && (locPhi <= locMeanFunc->GetXmax()) && (locPhi >= locSigmaFunc->GetXmin()) && (locPhi <= locSigmaFunc->GetXmax()))
	{
		double locMean = locMeanFunc->Eval(locPhi);
		double loc3Sigma = 3.0*locSigmaFunc->Eval(locPhi);
		if((locDeltaPhi >= (locMean - loc3Sigma)) && (locDeltaPhi <= (locMean + loc3Sigma)))
			return true;
	}

	return false;
}

Int_t Selector_TrackEff::Fill_ResolutionHists(double locRFWeight, int locBeamEnergyBin)
{
	//returns index of best track match (-1 if fails cut)
	double locMissingP = MissingP3->Mag();
	double locMissingTheta = MissingP3->Theta()*180.0/TMath::Pi();
	double locMissingPhi = MissingP3->Phi()*180.0/TMath::Pi();

	Int_t locBestTrackIndex = -1;
	Float_t locBestMatchFOM = -1.0f;
	for(UInt_t locTrackIndex = 0; locTrackIndex < NumUnusedTimeBased; ++locTrackIndex)
	{
		TVector3& locTrackP3 = *static_cast<TVector3*>(ReconP3->At(locTrackIndex));
		double locMeasuredP = locTrackP3.Mag();
		double locMeasuredTheta = locTrackP3.Theta()*180.0/TMath::Pi();
		double locMeasuredPhi = locTrackP3.Phi()*180.0/TMath::Pi();

		double locDeltaP = locMeasuredP - locMissingP;
		double locDeltaPOverP = locDeltaP/locMissingP;
		double locDeltaTheta = locMeasuredTheta - locMissingTheta;
		double locDeltaPhi = locMeasuredPhi - locMissingPhi;
		while(locDeltaPhi > 180.0)
			locDeltaPhi -= 360.0;
		while(locDeltaPhi < -180.0)
			locDeltaPhi += 360.0;

		bool locPassDeltaPhiCutFlag = Cut_DeltaPhi(locBeamEnergyBin, locMissingP, locMissingTheta, locMissingPhi, locDeltaPhi);
		bool locPassDeltaThetaCutFlag = Cut_DeltaTheta(locBeamEnergyBin, locMissingP, locMissingTheta, locMissingPhi, locDeltaTheta);
		bool locPassDeltaPOverPCutFlag = Cut_DeltaPOverP(locBeamEnergyBin, locMissingP, locMissingTheta, locMissingPhi, locDeltaPOverP);

		if((locRFWeight > 0.8) && (locBeamEnergyBin == 8) && (locMissingP >= 0.8) && (locMissingP <= 1.2) && (locDeltaPOverP > 0.3) && (locDeltaPOverP < 0.4))
		{
			cout << "Missing PID, Run, event# = " << dMissingPID << ", " << dRunNumber << ", " << EventNumber << endl;
			cout << "BeamE, MissingP, MissingTheta, MissingPhi, DeltaPOverP, ReconP, ReconTheta, ReconPhi = " << BeamEnergy << ", " << locMissingP << ", " << locMissingTheta;
			cout << ", " << locMissingPhi << ", " << locDeltaPOverP << ", " << locMeasuredP << ", " << locMeasuredTheta << ", " << locDeltaPhi << endl;
		}

		//Resolution
		if(locPassDeltaPhiCutFlag && locPassDeltaThetaCutFlag && (!gDebugFlag || locPassDeltaPOverPCutFlag))
		{
			dHistMap_Resolution_DeltaPVsP[locBeamEnergyBin]->Fill(locMissingP, locDeltaP, locRFWeight);
			dHistMap_Resolution_DeltaPVsTheta[locBeamEnergyBin]->Fill(locMissingTheta, locDeltaP, locRFWeight);
			dHistMap_Resolution_DeltaPVsPhi[locBeamEnergyBin]->Fill(locMissingPhi, locDeltaP, locRFWeight);

			dHistMap_Resolution_DeltaPOverPVsP[locBeamEnergyBin]->Fill(locMissingP, locDeltaPOverP, locRFWeight);
			dHistMap_Resolution_DeltaPOverPVsTheta[locBeamEnergyBin]->Fill(locMissingTheta, locDeltaPOverP, locRFWeight);
			dHistMap_Resolution_DeltaPOverPVsPhi[locBeamEnergyBin]->Fill(locMissingPhi, locDeltaPOverP, locRFWeight);
		}

		if(locPassDeltaPOverPCutFlag && locPassDeltaThetaCutFlag && (!gDebugFlag || locPassDeltaPhiCutFlag))
		{
			dHistMap_Resolution_DeltaPhiVsP[locBeamEnergyBin]->Fill(locMissingP, locDeltaPhi, locRFWeight);
			dHistMap_Resolution_DeltaPhiVsTheta[locBeamEnergyBin]->Fill(locMissingTheta, locDeltaPhi, locRFWeight);
			dHistMap_Resolution_DeltaPhiVsPhi[locBeamEnergyBin]->Fill(locMissingPhi, locDeltaPhi, locRFWeight);
		}

		if(locPassDeltaPhiCutFlag && locPassDeltaPOverPCutFlag && (!gDebugFlag || locPassDeltaThetaCutFlag))
		{
			dHistMap_Resolution_DeltaThetaVsP[locBeamEnergyBin]->Fill(locMissingP, locDeltaTheta, locRFWeight);
			dHistMap_Resolution_DeltaThetaVsTheta[locBeamEnergyBin]->Fill(locMissingTheta, locDeltaTheta, locRFWeight);
			dHistMap_Resolution_DeltaThetaVsPhi[locBeamEnergyBin]->Fill(locMissingPhi, locDeltaTheta, locRFWeight);
		}

		if(locPassDeltaPOverPCutFlag && locPassDeltaPhiCutFlag && locPassDeltaThetaCutFlag)
		{
			if((ReconMatchFOM[locTrackIndex] > locBestMatchFOM) || (locBestTrackIndex == -1))
			{
				locBestTrackIndex = locTrackIndex;
				locBestMatchFOM = ReconMatchFOM[locTrackIndex];
			}
		}
	}

	return locBestTrackIndex;
}

void Selector_TrackEff::Fill_MissingP3SigmaHists(double locRFWeight, int locBeamEnergyBin)
{
	double locMissingP = MissingP3->Mag();
	double locMissingTheta = MissingP3->Theta()*180.0/TMath::Pi();
	double locMissingPhi = MissingP3->Phi()*180.0/TMath::Pi();

	//Transform errors from pxyz to pthetaphi

	//Build error matrix
	TMatrixFSym locMissingP3Matrix(3);
	locMissingP3Matrix(0, 0) = MissingP3_CovPxPx;
	locMissingP3Matrix(0, 1) = MissingP3_CovPxPy;
	locMissingP3Matrix(0, 2) = MissingP3_CovPxPz;
	locMissingP3Matrix(1, 0) = MissingP3_CovPxPy;
	locMissingP3Matrix(1, 1) = MissingP3_CovPyPy;
	locMissingP3Matrix(1, 2) = MissingP3_CovPyPz;
	locMissingP3Matrix(2, 0) = MissingP3_CovPxPz;
	locMissingP3Matrix(2, 1) = MissingP3_CovPyPz;
	locMissingP3Matrix(2, 2) = MissingP3_CovPzPz;

	//calc p error
	TMatrixF locPJacobian(1, 3);
	locPJacobian(0, 0) = MissingP3->Px()/locMissingP;
	locPJacobian(0, 1) = MissingP3->Py()/locMissingP;
	locPJacobian(0, 2) = MissingP3->Pz()/locMissingP;
	TMatrixFSym locTempMatrix(locMissingP3Matrix);
	float locMissingPSigma = (locTempMatrix.Similarity(locPJacobian))(0, 0);

	//calc theta error
	TMatrixF locThetaJacobian(1, 3);
	locThetaJacobian(0, 0) = MissingP3->Px()*MissingP3->Pz()/(locMissingP*locMissingP*MissingP3->Perp());
	locThetaJacobian(0, 1) = MissingP3->Py()*MissingP3->Pz()/(locMissingP*locMissingP*MissingP3->Perp());
	locThetaJacobian(0, 2) = -1.0*MissingP3->Perp()/(locMissingP*locMissingP);
	TMatrixFSym locTempMatrix2(locMissingP3Matrix);
	float locMissingThetaSigma = 180.0*(locTempMatrix2.Similarity(locThetaJacobian))(0, 0)/TMath::Pi();

	//calc phi error
	TMatrixF locPhiJacobian(1, 3);
	locPhiJacobian(0, 0) = -1.0*MissingP3->Py()/(MissingP3->Perp2());
	locPhiJacobian(0, 1) = MissingP3->Px()/(MissingP3->Perp2());
	locPhiJacobian(0, 2) = 0.0;
	TMatrixFSym locTempMatrix3(locMissingP3Matrix);
	float locMissingPhiSigma = 180.0*(locTempMatrix3.Similarity(locPhiJacobian))(0, 0)/TMath::Pi();

	dHistMap_MissingPSigmaVsP[locBeamEnergyBin]->Fill(locMissingP, locMissingPSigma, locRFWeight);
	dHistMap_MissingPSigmaVsTheta[locBeamEnergyBin]->Fill(locMissingTheta, locMissingPSigma, locRFWeight);
	dHistMap_MissingPSigmaVsPhi[locBeamEnergyBin]->Fill(locMissingPhi, locMissingPSigma, locRFWeight);

	dHistMap_MissingPSigmaOverPVsP[locBeamEnergyBin]->Fill(locMissingP, locMissingPSigma/locMissingP, locRFWeight);
	dHistMap_MissingPSigmaOverPVsTheta[locBeamEnergyBin]->Fill(locMissingTheta, locMissingPSigma/locMissingP, locRFWeight);
	dHistMap_MissingPSigmaOverPVsPhi[locBeamEnergyBin]->Fill(locMissingPhi, locMissingPSigma/locMissingP, locRFWeight);

	dHistMap_MissingThetaSigmaVsP[locBeamEnergyBin]->Fill(locMissingP, locMissingThetaSigma, locRFWeight);
	dHistMap_MissingThetaSigmaVsTheta[locBeamEnergyBin]->Fill(locMissingTheta, locMissingThetaSigma, locRFWeight);
	dHistMap_MissingThetaSigmaVsPhi[locBeamEnergyBin]->Fill(locMissingPhi, locMissingThetaSigma, locRFWeight);

	dHistMap_MissingPhiSigmaVsP[locBeamEnergyBin]->Fill(locMissingP, locMissingPhiSigma, locRFWeight);
	dHistMap_MissingPhiSigmaVsTheta[locBeamEnergyBin]->Fill(locMissingTheta, locMissingPhiSigma, locRFWeight);
	dHistMap_MissingPhiSigmaVsPhi[locBeamEnergyBin]->Fill(locMissingPhi, locMissingPhiSigma, locRFWeight);
}

void Selector_TrackEff::Fill_MassHists(double locRFWeight, Int_t locBestTrackMatchIndex)
{
	//Total Missing P4
	if(locBestTrackMatchIndex != -1)
	{
		TVector3& locTrackP3 = *static_cast<TVector3*>(ReconP3->At(locBestTrackMatchIndex));
		TLorentzVector locReconstructedMissingP4 = TLorentzVector((*MissingP3) - locTrackP3, MeasuredMissingE[locBestTrackMatchIndex]);

		dHist_ReconMissingMassSquaredVsBeamEnergy->Fill(BeamEnergy, locReconstructedMissingP4.M2(), locRFWeight);
		dHist_TotalMissingEnergyVsBeamEnergy->Fill(BeamEnergy, locReconstructedMissingP4.E(), locRFWeight);
		dHist_TotalMissingPyVsPx->Fill(locReconstructedMissingP4.Px(), locReconstructedMissingP4.Py(), locRFWeight);
		dHist_TotalMissingPtVsBeamEnergy->Fill(BeamEnergy, locReconstructedMissingP4.Perp(), locRFWeight);
	}

	dHist_BeamRFDeltaT->Fill(BeamRFDeltaT);

	dHist_BeamEnergy->Fill(BeamEnergy, locRFWeight);
	dHist_MissingMass->Fill(MissingMassSquared, locRFWeight);
	dHist_MissingMassVsBeamEnergy->Fill(BeamEnergy, MissingMassSquared, locRFWeight);

	if(locBestTrackMatchIndex != -1)
	{
		dHist_BeamEnergy_Found->Fill(BeamEnergy, locRFWeight);
		dHist_MissingMassVsBeamEnergy_Found->Fill(BeamEnergy, MissingMassSquared, locRFWeight);
	}
	else //track missing
	{
		dHist_BeamEnergy_Missing->Fill(BeamEnergy, locRFWeight);
		dHist_MissingMassVsBeamEnergy_Missing->Fill(BeamEnergy, MissingMassSquared, locRFWeight);
	}
}

void Selector_TrackEff::Fill_MatchingHists(void)
{
	double locMissingP = MissingP3->Mag();
	double locMissingTheta = MissingP3->Theta()*180.0/TMath::Pi();
	double locMissingPhi = MissingP3->Phi()*180.0/TMath::Pi();

	for(UInt_t locTrackIndex = 0; locTrackIndex < NumUnusedTimeBased; ++locTrackIndex)
	{
		TVector3& locTrackP3 = *static_cast<TVector3*>(ReconP3->At(locTrackIndex));

		double locMeasuredP = locTrackP3.Mag();
		double locMeasuredTheta = locTrackP3.Theta()*180.0/TMath::Pi();
		double locMeasuredPhi = locTrackP3.Phi()*180.0/TMath::Pi();

		double locDeltaPOverP = (locMeasuredP - locMissingP)/locMissingP;
		double locDeltaTheta = locMeasuredTheta - locMissingTheta;
		double locDeltaPhi = locMeasuredPhi - locMissingPhi;
		while(locDeltaPhi > 180.0)
			locDeltaPhi -= 360.0;
		while(locDeltaPhi < -180.0)
			locDeltaPhi += 360.0;

		dHist_MatchingFOM->Fill(ReconMatchFOM[locTrackIndex]);

		dHist_MatchingFOMVsDeltaPOverP->Fill(locDeltaPOverP, ReconMatchFOM[locTrackIndex]);
		dHist_MatchingFOMVsDeltaTheta->Fill(locDeltaTheta, ReconMatchFOM[locTrackIndex]);
		dHist_MatchingFOMVsDeltaPhi->Fill(locDeltaPhi, ReconMatchFOM[locTrackIndex]);

		dHist_MatchingFOMVsDeltaPOverP_LogY->Fill(locDeltaPOverP, ReconMatchFOM[locTrackIndex]);
		dHist_MatchingFOMVsDeltaTheta_LogY->Fill(locDeltaTheta, ReconMatchFOM[locTrackIndex]);
		dHist_MatchingFOMVsDeltaPhi_LogY->Fill(locDeltaPhi, ReconMatchFOM[locTrackIndex]);
	}
}

void Selector_TrackEff::Fill_EffHists(double locRFWeight, bool locTrackFoundFlag, bool locSignalPeakFlag, int locBeamEnergyBin)
{
	int locVertexZBin = (ComboVertexZ - dMinVertexZ)/dVertexZBinSize;
	if((locVertexZBin < 0) || (locVertexZBin >= dNumVertexZBins))
		return;

	double locMissingP = MissingP3->Mag();
	double locMissingTheta = MissingP3->Theta()*180.0/TMath::Pi();
	double locMissingPhi = MissingP3->Phi()*180.0/TMath::Pi();

	//Efficiency
	if(locTrackFoundFlag)
	{
		//Found
		dHistMap_TrackFound_PVsTheta[locSignalPeakFlag][locBeamEnergyBin][locVertexZBin]->Fill(locMissingTheta, locMissingP, locRFWeight);
		dHistMap_TrackFound_PhiVsTheta[locSignalPeakFlag][locBeamEnergyBin][locVertexZBin]->Fill(locMissingTheta, locMissingPhi, locRFWeight);
	}
	else
	{
		//Missing
		dHistMap_TrackMissing_PVsTheta[locSignalPeakFlag][locBeamEnergyBin][locVertexZBin]->Fill(locMissingTheta, locMissingP, locRFWeight);
		dHistMap_TrackMissing_PhiVsTheta[locSignalPeakFlag][locBeamEnergyBin][locVertexZBin]->Fill(locMissingTheta, locMissingPhi, locRFWeight);
	}
}

		map<int, TH2D*> dHistMap_MissingPSigmaVsP;
		map<int, TH2D*> dHistMap_MissingPSigmaVsTheta;
		map<int, TH2D*> dHistMap_MissingPSigmaVsPhi;
		map<int, TH2D*> dHistMap_MissingThetaSigmaVsP;
		map<int, TH2D*> dHistMap_MissingThetaSigmaVsTheta;
		map<int, TH2D*> dHistMap_MissingThetaSigmaVsPhi;
		map<int, TH2D*> dHistMap_MissingPhiSigmaVsP;
		map<int, TH2D*> dHistMap_MissingPhiSigmaVsTheta;
		map<int, TH2D*> dHistMap_MissingPhiSigmaVsPhi;
   Float_t         MissingP3_CovPxPx;
   Float_t         MissingP3_CovPxPy;
   Float_t         MissingP3_CovPxPz;
   Float_t         MissingP3_CovPyPy;
   Float_t         MissingP3_CovPyPz;
   Float_t         MissingP3_CovPzPz;


void Selector_TrackEff::SlaveTerminate()
{
	// The SlaveTerminate() function is called after all entries or objects
	// have been processed. When running with PROOF SlaveTerminate() is called
	// on each slave server.
	Finalize();
}

void Selector_TrackEff::Terminate()
{
	// The Terminate() function is the last function to be called during
	// a query. It always runs on the client, it can be used to present
	// the results graphically or save the results to file.
	if(dFinalizeFlag)
		return;
	Finalize();
}

void Selector_TrackEff::Finalize(void)
{
	dOutputFile->Write();
	dOutputFile->Close();
	if(dProofFile != NULL)
		fOutput->Add(dProofFile);
	dFinalizeFlag = true;
}