// $Id$
//
//    File: DCustomAction_TrackPStudies.cc
// Created: Mon May  5 19:51:25 CDT 2014
// Creator: pmatt (on Darwin pmattLaptop 10.8.0 i386)
//

#include "DCustomAction_TrackPStudies.h"

void DCustomAction_TrackPStudies::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. 

	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();

		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("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("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("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"));
	}
	japp->RootUnLock(); //RELEASE ROOT LOCK!!
}

bool DCustomAction_TrackPStudies::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
{
	//Optional: check whether the user wanted to use the kinematic fit results when performing this action
	//bool locUseKinFitResultsFlag = Get_UseKinFitResultsFlag();

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

	deque<const DKinematicData*> locChargedParticles, locNeutralParticles, locAllParticles;
	locParticleCombo->Get_DetectedFinalChargedParticles_Measured(locChargedParticles);
	locParticleCombo->Get_DetectedFinalNeutralParticles_Measured(locNeutralParticles);
	locParticleCombo->Get_DetectedFinalParticles_Measured(locAllParticles);

	DVector3 locTotalMomentum;
	for(size_t loc_i = 0; loc_i < locChargedParticles.size(); ++loc_i)
	{
		const DChargedTrackHypothesis* locChargedTrackHypothesis = static_cast<const DChargedTrackHypothesis*>(locChargedParticles[loc_i]);
		locTotalMomentum += locChargedTrackHypothesis->momentum();
	}
	for(size_t loc_i = 0; loc_i < locNeutralParticles.size(); ++loc_i)
	{
		const DNeutralParticleHypothesis* locNeutralParticleHypothesis = static_cast<const DNeutralParticleHypothesis*>(locNeutralParticles[loc_i]);
		locTotalMomentum += locNeutralParticleHypothesis->momentum();
	}

	double locBeamEnergy = locParticleCombo->Get_ParticleComboStep(0)->Get_InitialParticle_Measured()->momentum().Mag();

	//Optional: Fill histograms
	japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
	{
		dHist_TotalP->Fill(locTotalMomentum.Mag());
		dHist_TotalPVsBeamE->Fill(locBeamEnergy, locTotalMomentum.Mag());

		dHist_TotalPt->Fill(locTotalMomentum.Perp());
		dHist_TotalPyVsPx->Fill(locTotalMomentum.Px(), locTotalMomentum.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!)
}