// $Id$
//
//    File: DCustomAction_HistTotalPVsBeamE.cc
// Created: Mon May  5 13:46:49 PDT 2014
// Creator: pmatt (on Darwin pmattLaptop 10.8.0 i386)
//

#include "DCustomAction_HistTotalPVsBeamE.h"

void DCustomAction_HistTotalPVsBeamE::Initialize(JEventLoop* locEventLoop)
{
	//Optional: Create histograms and/or modify member variables.
	//Create any histograms/trees/etc. within a ROOT lock. 
		//This is so that when running multithreaded, only one thread is writing to the ROOT file at a time. 

	//When creating a reaction-independent action, only modify member variables within a ROOT lock. 
		//Objects created within a plugin (such as reaction-independent actions) can be accessed by many threads simultaneously. 

	japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
	{
		// Optional: Useful utility functions.
		// locEventLoop->GetSingle(dAnalysisUtilities);

		//Required: Create a folder in the ROOT output file that will contain all of the output ROOT objects (if any) for this action.
			//If another thread has already created the folder, it just changes to it. 
		CreateAndChangeTo_ActionDirectory();

		//total p
		if(gDirectory->Get("TotalP") == NULL) //check to see if already created by another thread
			dHist_TotalP = new TH1D("TotalP", "Total Momentum (GeV/c)", 2000, 0.0, 20.0);
		else //already created by another thread
			dHist_TotalP = static_cast<TH1D*>(gDirectory->Get("TotalP"));

		if(gDirectory->Get("TrueTotalP") == NULL) //check to see if already created by another thread
			dHist_TotalP_True = new TH1D("TrueTotalP", "Total Momentum (GeV/c)", 2000, 0.0, 20.0);
		else //already created by another thread
			dHist_TotalP_True = static_cast<TH1D*>(gDirectory->Get("TrueTotalP"));

		//total p vs beam E
		if(gDirectory->Get("TotalPVsBeamE") == NULL) //check to see if already created by another thread
			dHist_TotalPVsBeamE = new TH2D("TotalPVsBeamE", ";Beam Energy (GeV);Total Momentum (GeV/c)", 600, 0.0, 12.0, 2000, 0.0, 20.0);
		else //already created by another thread
			dHist_TotalPVsBeamE = static_cast<TH2D*>(gDirectory->Get("TotalPVsBeamE"));

		if(gDirectory->Get("TrueTotalPVsBeamE") == NULL) //check to see if already created by another thread
			dHist_TotalPVsBeamE_True = new TH2D("TrueTotalPVsBeamE", ";Beam Energy (GeV);Total Momentum (GeV/c)", 600, 0.0, 12.0, 2000, 0.0, 20.0);
		else //already created by another thread
			dHist_TotalPVsBeamE_True = static_cast<TH2D*>(gDirectory->Get("TrueTotalPVsBeamE"));

		//total pt
		if(gDirectory->Get("TotalPt") == NULL) //check to see if already created by another thread
			dHist_TotalPt = new TH1D("TotalPt", "Total Transverse Momentum (GeV/c)", 5000, 0.0, 5.0);
		else //already created by another thread
			dHist_TotalPt = static_cast<TH1D*>(gDirectory->Get("TotalPt"));

		if(gDirectory->Get("TrueTotalPt") == NULL) //check to see if already created by another thread
			dHist_TotalPt_True = new TH1D("TrueTotalPt", "Total Transverse Momentum (GeV/c)", 5000, 0.0, 5.0);
		else //already created by another thread
			dHist_TotalPt_True = static_cast<TH1D*>(gDirectory->Get("TrueTotalPt"));

		//total py vs px
		if(gDirectory->Get("TotalPyVsPx") == NULL) //check to see if already created by another thread
			dHist_TotalPyVsPx = new TH2D("TotalPyVsPx", ";Total p_{x} (GeV/c);Total p_{y} (GeV/c)", 1000, -2.0, 2.0, 1000, -2.0, 2.0);
		else //already created by another thread
			dHist_TotalPyVsPx = static_cast<TH2D*>(gDirectory->Get("TotalPyVsPx"));

		if(gDirectory->Get("TrueTotalPyVsPx") == NULL) //check to see if already created by another thread
			dHist_TotalPyVsPx_True = new TH2D("TrueTotalPyVsPx", ";Total p_{x} (GeV/c);Total p_{y} (GeV/c)", 1000, -2.0, 2.0, 1000, -2.0, 2.0);
		else //already created by another thread
			dHist_TotalPyVsPx_True = static_cast<TH2D*>(gDirectory->Get("TrueTotalPyVsPx"));

		CreateAndChangeTo_Directory("ParticleFOMCuts", "ParticleFOMCuts");
		{
			//total p
			if(gDirectory->Get("TotalP") == NULL) //check to see if already created by another thread
				dHist_TotalP_ParticleFOMCuts = new TH1D("TotalP", "Total Momentum (GeV/c)", 2000, 0.0, 20.0);
			else //already created by another thread
				dHist_TotalP_ParticleFOMCuts = static_cast<TH1D*>(gDirectory->Get("TotalP"));

			if(gDirectory->Get("TrueTotalP") == NULL) //check to see if already created by another thread
				dHist_TotalP_True_ParticleFOMCuts = new TH1D("TrueTotalP", "Total Momentum (GeV/c)", 2000, 0.0, 20.0);
			else //already created by another thread
				dHist_TotalP_True_ParticleFOMCuts = static_cast<TH1D*>(gDirectory->Get("TrueTotalP"));

			//total p vs beam E
			if(gDirectory->Get("TotalPVsBeamE") == NULL) //check to see if already created by another thread
				dHist_TotalPVsBeamE_ParticleFOMCuts = new TH2D("TotalPVsBeamE", ";Beam Energy (GeV);Total Momentum (GeV/c)", 600, 0.0, 12.0, 2000, 0.0, 20.0);
			else //already created by another thread
				dHist_TotalPVsBeamE_ParticleFOMCuts = static_cast<TH2D*>(gDirectory->Get("TotalPVsBeamE"));

			if(gDirectory->Get("TrueTotalPVsBeamE") == NULL) //check to see if already created by another thread
				dHist_TotalPVsBeamE_True_ParticleFOMCuts = new TH2D("TrueTotalPVsBeamE", ";Beam Energy (GeV);Total Momentum (GeV/c)", 600, 0.0, 12.0, 2000, 0.0, 20.0);
			else //already created by another thread
				dHist_TotalPVsBeamE_True_ParticleFOMCuts = static_cast<TH2D*>(gDirectory->Get("TrueTotalPVsBeamE"));

			//total pt
			if(gDirectory->Get("TotalPt") == NULL) //check to see if already created by another thread
				dHist_TotalPt_ParticleFOMCuts = new TH1D("TotalPt", "Total Transverse Momentum (GeV/c)", 5000, 0.0, 5.0);
			else //already created by another thread
				dHist_TotalPt_ParticleFOMCuts = static_cast<TH1D*>(gDirectory->Get("TotalPt"));

			if(gDirectory->Get("TrueTotalPt") == NULL) //check to see if already created by another thread
				dHist_TotalPt_True_ParticleFOMCuts = new TH1D("TrueTotalPt", "Total Transverse Momentum (GeV/c)", 5000, 0.0, 5.0);
			else //already created by another thread
				dHist_TotalPt_True_ParticleFOMCuts = static_cast<TH1D*>(gDirectory->Get("TrueTotalPt"));

			//total py vs px
			if(gDirectory->Get("TotalPyVsPx") == NULL) //check to see if already created by another thread
				dHist_TotalPyVsPx_ParticleFOMCuts = new TH2D("TotalPyVsPx", ";Total p_{x} (GeV/c);Total p_{y} (GeV/c)", 1000, -2.0, 2.0, 1000, -2.0, 2.0);
			else //already created by another thread
				dHist_TotalPyVsPx_ParticleFOMCuts = static_cast<TH2D*>(gDirectory->Get("TotalPyVsPx"));

			if(gDirectory->Get("TrueTotalPyVsPx") == NULL) //check to see if already created by another thread
				dHist_TotalPyVsPx_True_ParticleFOMCuts = new TH2D("TrueTotalPyVsPx", ";Total p_{x} (GeV/c);Total p_{y} (GeV/c)", 1000, -2.0, 2.0, 1000, -2.0, 2.0);
			else //already created by another thread
				dHist_TotalPyVsPx_True_ParticleFOMCuts = static_cast<TH2D*>(gDirectory->Get("TrueTotalPyVsPx"));
		}
		gDirectory->cd("..");

		CreateAndChangeTo_Directory("HasDetectorHits", "HasDetectorHits");
		{
			//total p
			if(gDirectory->Get("TotalP") == NULL) //check to see if already created by another thread
				dHist_TotalP_HasDetectorHits = new TH1D("TotalP", "Total Momentum (GeV/c)", 2000, 0.0, 20.0);
			else //already created by another thread
				dHist_TotalP_HasDetectorHits = static_cast<TH1D*>(gDirectory->Get("TotalP"));

			if(gDirectory->Get("TrueTotalP") == NULL) //check to see if already created by another thread
				dHist_TotalP_True_HasDetectorHits = new TH1D("TrueTotalP", "Total Momentum (GeV/c)", 2000, 0.0, 20.0);
			else //already created by another thread
				dHist_TotalP_True_HasDetectorHits = static_cast<TH1D*>(gDirectory->Get("TrueTotalP"));

			//total p vs beam E
			if(gDirectory->Get("TotalPVsBeamE") == NULL) //check to see if already created by another thread
				dHist_TotalPVsBeamE_HasDetectorHits = new TH2D("TotalPVsBeamE", ";Beam Energy (GeV);Total Momentum (GeV/c)", 600, 0.0, 12.0, 2000, 0.0, 20.0);
			else //already created by another thread
				dHist_TotalPVsBeamE_HasDetectorHits = static_cast<TH2D*>(gDirectory->Get("TotalPVsBeamE"));

			if(gDirectory->Get("TrueTotalPVsBeamE") == NULL) //check to see if already created by another thread
				dHist_TotalPVsBeamE_True_HasDetectorHits = new TH2D("TrueTotalPVsBeamE", ";Beam Energy (GeV);Total Momentum (GeV/c)", 600, 0.0, 12.0, 2000, 0.0, 20.0);
			else //already created by another thread
				dHist_TotalPVsBeamE_True_HasDetectorHits = static_cast<TH2D*>(gDirectory->Get("TrueTotalPVsBeamE"));

			//total pt
			if(gDirectory->Get("TotalPt") == NULL) //check to see if already created by another thread
				dHist_TotalPt_HasDetectorHits = new TH1D("TotalPt", "Total Transverse Momentum (GeV/c)", 5000, 0.0, 5.0);
			else //already created by another thread
				dHist_TotalPt_HasDetectorHits = static_cast<TH1D*>(gDirectory->Get("TotalPt"));

			if(gDirectory->Get("TrueTotalPt") == NULL) //check to see if already created by another thread
				dHist_TotalPt_True_HasDetectorHits = new TH1D("TrueTotalPt", "Total Transverse Momentum (GeV/c)", 5000, 0.0, 5.0);
			else //already created by another thread
				dHist_TotalPt_True_HasDetectorHits = static_cast<TH1D*>(gDirectory->Get("TrueTotalPt"));

			//total py vs px
			if(gDirectory->Get("TotalPyVsPx") == NULL) //check to see if already created by another thread
				dHist_TotalPyVsPx_HasDetectorHits = new TH2D("TotalPyVsPx", ";Total p_{x} (GeV/c);Total p_{y} (GeV/c)", 1000, -2.0, 2.0, 1000, -2.0, 2.0);
			else //already created by another thread
				dHist_TotalPyVsPx_HasDetectorHits = static_cast<TH2D*>(gDirectory->Get("TotalPyVsPx"));

			if(gDirectory->Get("TrueTotalPyVsPx") == NULL) //check to see if already created by another thread
				dHist_TotalPyVsPx_True_HasDetectorHits = new TH2D("TrueTotalPyVsPx", ";Total p_{x} (GeV/c);Total p_{y} (GeV/c)", 1000, -2.0, 2.0, 1000, -2.0, 2.0);
			else //already created by another thread
				dHist_TotalPyVsPx_True_HasDetectorHits = static_cast<TH2D*>(gDirectory->Get("TrueTotalPyVsPx"));
		}
		gDirectory->cd("..");
	}
	japp->RootUnLock(); //RELEASE ROOT LOCK!!
}

bool DCustomAction_HistTotalPVsBeamE::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
{
	//Expect locParticleCombo to be NULL since this is a reaction-independent action.

	//Optional: Quit the action if it has already been executed this event (else may result in double-counting when filling histograms)
	if(Get_NumPreviousParticleCombos() != 0)
		return true;

	vector<const DBeamPhoton*> locBeamPhotons;
	locEventLoop->Get(locBeamPhotons);

	vector<const DChargedTrack*> locChargedTracks;
	locEventLoop->Get(locChargedTracks);

	vector<const DNeutralParticle*> locNeutralParticles;
	locEventLoop->Get(locNeutralParticles);

	vector<const DMCThrown*> locMCThrowns;
	locEventLoop->Get(locMCThrowns, "FinalState");

	const DMCThrownMatching* locMCThrownMatching = NULL;
	locEventLoop->GetSingle(locMCThrownMatching);

	DVector3 locMomentum, locMomentum_True, locMomentum_ParticleFOMCuts, locMomentum_True_ParticleFOMCuts, locMomentum_HasDetectorHits, locMomentum_True_HasDetectorHits;
	for(size_t loc_i = 0; loc_i < locChargedTracks.size(); ++loc_i)
	{
		const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTracks[loc_i]->Get_BestFOM();
		double locTrackingFOM = TMath::Prob(locChargedTrackHypothesis->dChiSq_Track, locChargedTrackHypothesis->dNDF_Track);

		bool locHasDetectorHit = true;
		do
		{
			if(locChargedTrackHypothesis->dSCHitMatchParams.dTrackTimeBased != NULL)
				break;
			if(locChargedTrackHypothesis->dTOFHitMatchParams.dTrackTimeBased != NULL)
				break;
			if(locChargedTrackHypothesis->dBCALShowerMatchParams.dTrackTimeBased != NULL)
				break;
			if(locChargedTrackHypothesis->dFCALShowerMatchParams.dTrackTimeBased != NULL)
				break;
			locHasDetectorHit = false;
		}
		while(false);

		locMomentum += locChargedTrackHypothesis->momentum();
		if(locTrackingFOM >= 5.73303E-7)
		{
			locMomentum_ParticleFOMCuts += locChargedTrackHypothesis->momentum();
			if(locHasDetectorHit)
				locMomentum_HasDetectorHits += locChargedTrackHypothesis->momentum();
		}

		if(locMCThrownMatching == NULL)
			continue;
		double locMatchFOM = 0.0;
		const DMCThrown* locMCThrown = locMCThrownMatching->Get_MatchingMCThrown(locChargedTracks[loc_i], locMatchFOM);
		if((locMCThrown == NULL) || (locMatchFOM < dMinThrownMatchFOM))
			continue;
		Particle_t locPID = locMCThrown->PID();

		locChargedTrackHypothesis = locChargedTracks[loc_i]->Get_Hypothesis(locPID);
		if(locChargedTrackHypothesis == NULL)
			locChargedTrackHypothesis = locChargedTracks[loc_i]->Get_BestFOM();

		locMomentum_True += locChargedTrackHypothesis->momentum();
		if(locTrackingFOM >= 5.73303E-7)
		{
			locMomentum_True_ParticleFOMCuts += locChargedTrackHypothesis->momentum();
			if(locHasDetectorHit)
				locMomentum_True_HasDetectorHits += locChargedTrackHypothesis->momentum();
		}
	}

	for(size_t loc_i = 0; loc_i < locNeutralParticles.size(); ++loc_i)
	{
		const DNeutralParticleHypothesis* locNeutralParticleHypothesis = locNeutralParticles[loc_i]->Get_BestFOM();
		locMomentum += locNeutralParticleHypothesis->momentum();
		if(locNeutralParticleHypothesis->dFOM >= 5.73303E-7)
		{
			locMomentum_ParticleFOMCuts += locNeutralParticleHypothesis->momentum();
			locMomentum_HasDetectorHits += locNeutralParticleHypothesis->momentum();
		}

		if(locMCThrownMatching == NULL)
			continue;
		double locMatchFOM = 0.0;
		const DMCThrown* locMCThrown = locMCThrownMatching->Get_MatchingMCThrown(locNeutralParticles[loc_i], locMatchFOM);
		if((locMCThrown == NULL) || (locMatchFOM < dMinThrownMatchFOM))
			continue;
		Particle_t locPID = locMCThrown->PID();

		locNeutralParticleHypothesis = locNeutralParticles[loc_i]->Get_Hypothesis(locPID);
		if(locNeutralParticleHypothesis == NULL)
			locNeutralParticleHypothesis = locNeutralParticles[loc_i]->Get_BestFOM();

		locMomentum_True += locNeutralParticleHypothesis->momentum();
		if(locNeutralParticleHypothesis->dFOM >= 5.73303E-7)
		{
			locMomentum_True_ParticleFOMCuts += locNeutralParticleHypothesis->momentum();
			locMomentum_True_HasDetectorHits += locNeutralParticleHypothesis->momentum();
		}
	}

	unsigned int locNumReconstructedTrueParticles = 0;
	if(locMCThrownMatching != NULL)
	{
		double locMatchFOM = 0.0;
		for(size_t loc_i = 0; loc_i < locChargedTracks.size(); ++loc_i)
		{
			const DMCThrown* locMCThrown = locMCThrownMatching->Get_MatchingMCThrown(locChargedTracks[loc_i], locMatchFOM);
			if((locMCThrown == NULL) || (locMatchFOM < dMinThrownMatchFOM))
				continue;
			++locNumReconstructedTrueParticles;
		}
		for(size_t loc_i = 0; loc_i < locNeutralParticles.size(); ++loc_i)
		{
			const DMCThrown* locMCThrown = locMCThrownMatching->Get_MatchingMCThrown(locNeutralParticles[loc_i], locMatchFOM);
			if((locMCThrown == NULL) || (locMatchFOM < dMinThrownMatchFOM))
				continue;
			++locNumReconstructedTrueParticles;
		}
	}

	double locBeamEnergy = locBeamPhotons[0]->energy();
	bool locTrueOKFlag = ((!locMCThrowns.empty()) && (locNumReconstructedTrueParticles >= (locMCThrowns.size() - 1)));

	//Optional: Fill histograms
	japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
	{
		// Fill any histograms here
		dHist_TotalP->Fill(locMomentum.Mag());
		dHist_TotalPVsBeamE->Fill(locBeamEnergy, locMomentum.Mag());
		dHist_TotalPt->Fill(locMomentum.Perp());
		dHist_TotalPyVsPx->Fill(locMomentum.Px(), locMomentum.Py());

		dHist_TotalP_ParticleFOMCuts->Fill(locMomentum_ParticleFOMCuts.Mag());
		dHist_TotalPVsBeamE_ParticleFOMCuts->Fill(locBeamEnergy, locMomentum_ParticleFOMCuts.Mag());
		dHist_TotalPt_ParticleFOMCuts->Fill(locMomentum_ParticleFOMCuts.Perp());
		dHist_TotalPyVsPx_ParticleFOMCuts->Fill(locMomentum_ParticleFOMCuts.Px(), locMomentum_ParticleFOMCuts.Py());

		dHist_TotalP_HasDetectorHits->Fill(locMomentum_HasDetectorHits.Mag());
		dHist_TotalPVsBeamE_HasDetectorHits->Fill(locBeamEnergy, locMomentum_HasDetectorHits.Mag());
		dHist_TotalPt_HasDetectorHits->Fill(locMomentum_HasDetectorHits.Perp());
		dHist_TotalPyVsPx_HasDetectorHits->Fill(locMomentum_HasDetectorHits.Px(), locMomentum_HasDetectorHits.Py());

		if(locTrueOKFlag)
		{
			dHist_TotalP_True->Fill(locMomentum_True.Mag());
			dHist_TotalPVsBeamE_True->Fill(locBeamEnergy, locMomentum_True.Mag());
			dHist_TotalPt_True->Fill(locMomentum_True.Perp());
			dHist_TotalPyVsPx_True->Fill(locMomentum_True.Px(), locMomentum_True.Py());

			dHist_TotalP_True_ParticleFOMCuts->Fill(locMomentum_True_ParticleFOMCuts.Mag());
			dHist_TotalPVsBeamE_True_ParticleFOMCuts->Fill(locBeamEnergy, locMomentum_True_ParticleFOMCuts.Mag());
			dHist_TotalPt_True_ParticleFOMCuts->Fill(locMomentum_True_ParticleFOMCuts.Perp());
			dHist_TotalPyVsPx_True_ParticleFOMCuts->Fill(locMomentum_True_ParticleFOMCuts.Px(), locMomentum_True_ParticleFOMCuts.Py());

			dHist_TotalP_True_HasDetectorHits->Fill(locMomentum_True_HasDetectorHits.Mag());
			dHist_TotalPVsBeamE_True_HasDetectorHits->Fill(locBeamEnergy, locMomentum_True_HasDetectorHits.Mag());
			dHist_TotalPt_True_HasDetectorHits->Fill(locMomentum_True_HasDetectorHits.Perp());
			dHist_TotalPyVsPx_True_HasDetectorHits->Fill(locMomentum_True_HasDetectorHits.Px(), locMomentum_True_HasDetectorHits.Py());
		}
	}
	japp->RootUnLock(); //RELEASE ROOT LOCK!!

	return true; //return false if you want to use this action to apply a cut (and it fails the cut!)
}