// $Id$
//
//    File: DCustomAction_InvariantMass_pippim4gamma.cc
// Created: Tue Jul  2 21:57:46 EDT 2013
// Creator: pmatt (on Darwin pmattLaptop 10.8.0 i386)
//

#include "DCustomAction_InvariantMass_pippim4gamma.h"

void DCustomAction_InvariantMass_pippim4gamma::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!!
	{
		//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();

		//	Create a histogram. 
		string locHistTitle = "#it{#gamma}#it{p}#rightarrow#it{#omega}#it{#pi}^{0}";
		locHistTitle += ";#it{#omega}#it{#pi}^{0} Invariant Mass (GeV/c^{2})";
		locHistTitle += ";# Combos / 2 MeV/c^{2}";
		if(gDirectory->Get("InvariantMass") == NULL) //check to see if already created by another thread
			dInvariantMassHist = new TH1D("InvariantMass", locHistTitle.c_str(), 2000, 0.0, 4.0);
		else //already created by another thread
			dInvariantMassHist = static_cast<TH1D*>(gDirectory->Get("InvariantMass"));
	}
	japp->RootUnLock(); //RELEASE ROOT LOCK!!
}

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

// 	const DParticleComboStep* locParticleComboStep0 = locParticleCombo->Get_ParticleComboStep(0);
// 	const DParticleComboStep* locParticleComboStep1 = locParticleCombo->Get_ParticleComboStep(1);
// 	const DParticleComboStep* locParticleComboStep2 = locParticleCombo->Get_ParticleComboStep(2);
// 	const DParticleComboStep* locParticleComboStep3 = locParticleCombo->Get_ParticleComboStep(3);
// 
// 	//first get measured particle objects, and check to see if combination of particles is unique
// 
// 	const DKinematicData* locPiPlus1 = locParticleComboStep1->Get_FinalParticle_Measured(0);
// 	const DKinematicData* locPiMinus1 = locParticleComboStep1->Get_FinalParticle_Measured(1);
// 
// 	const DKinematicData* locPhoton1 = locParticleComboStep2->Get_FinalParticle_Measured(0);
// 	const DKinematicData* locPhoton2 = locParticleComboStep2->Get_FinalParticle_Measured(1);
// 
// 	const DKinematicData* locPhoton3 = locParticleComboStep3->Get_FinalParticle_Measured(0);
// 	const DKinematicData* locPhoton4 = locParticleComboStep3->Get_FinalParticle_Measured(1);
// 
// 	set<const DKinematicData*> locCurrentParticles;
// 	locCurrentParticles.insert(locPiMinus1);
// 	locCurrentParticles.insert(locPiPlus1);
// 	locCurrentParticles.insert(locPhoton1);
// 	locCurrentParticles.insert(locPhoton2);
// 	locCurrentParticles.insert(locPhoton3);
// 	locCurrentParticles.insert(locPhoton4);
// 
// 	//if new event: clear past particles, else check if duplicate
// 	if(Get_NumPreviousParticleCombos() == 0)
// 		dPastParticles.clear();
// 	else //have had previous combos for this event, check to make sure particles used to compute this quantity aren't duplicate
// 	{
// 		for(size_t loc_i = 0; loc_i < dPastParticles.size(); ++loc_i)
// 		{
// 			if(locCurrentParticles != dPastParticles[loc_i])
// 				continue;
// 			return true; //duplicate combo of particles, don't fill histogram!
// 		}
// 	}
// 	dPastParticles.push_back(locCurrentParticles);
// 
// 	//calculate invariant mass
// 	DLorentzVector locP4;
// //	if(!locUseKinFitResultsFlag || (locParticleCombo->Get_KinFitResults() == NULL)) //measured or kinfit failed to converge
// 	{
// 		locP4 += locPiMinus1->lorentzMomentum();
// 		locP4 += locPiPlus1->lorentzMomentum();
// 		locP4 += locPhoton1->lorentzMomentum() + locPhoton2->lorentzMomentum() + locPhoton3->lorentzMomentum() + locPhoton4->lorentzMomentum();
// 	}
// /*	else //kinfit: get kinfit objects first
// 	{
// 		locPiPlus1 = locParticleComboStep0->Get_FinalParticle(2);
// 		locPiPlus2 = locParticleComboStep1->Get_FinalParticle(0);
// 		locPiMinus2 = locParticleComboStep1->Get_FinalParticle(1);
// 		const DKinematicData* locPiZero = locParticleComboStep1->Get_FinalParticle(2);
// 		locP4 += locPiMinus2->lorentzMomentum();
// 		locP4 += locPiPlus1->lorentzMomentum() + locPiPlus2->lorentzMomentum();
// 		locP4 += locPiZero->lorentzMomentum();
// 	}
// 	*/
// 	double locInvariantMass = locP4.M();
// 
// 	//FILL HISTOGRAMS
// 	//Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
// 	//Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
//  	japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
// 	{
// 		// Fill any histograms here
// 		dInvariantMassHist->Fill(locInvariantMass);
// 	}
//  	japp->RootUnLock(); //RELEASE ROOT LOCK!!


	//Get the step containing the particles
//	const DParticleComboStep* locParticleComboStep0 = locParticleCombo->Get_ParticleComboStep(0);
	const DParticleComboStep* locParticleComboStep0 = locParticleCombo->Get_ParticleComboStep(0);
//	const DParticleComboStep* locParticleComboStep1 = locParticleCombo->Get_ParticleComboStep(1);
	//Get the particles, using the kinematic-fit or measured data as requested
	const DKinematicData* locPiPlus = locUseKinFitResultsFlag ? locParticleComboStep0->Get_FinalParticle(0) : locParticleComboStep0->Get_FinalParticle_Measured(0);
	const DKinematicData* locPiMinus = locUseKinFitResultsFlag ? locParticleComboStep0->Get_FinalParticle(1) : locParticleComboStep0->Get_FinalParticle_Measured(1);
	const DKinematicData* locPhoton1 = locUseKinFitResultsFlag ? locParticleComboStep0->Get_FinalParticle(3) : locParticleComboStep0->Get_FinalParticle_Measured(3);
	const DKinematicData* locPhoton2 = locUseKinFitResultsFlag ? locParticleComboStep0->Get_FinalParticle(4) : locParticleComboStep0->Get_FinalParticle_Measured(4);
	const DKinematicData* locPhoton3 = locUseKinFitResultsFlag ? locParticleComboStep0->Get_FinalParticle(5) : locParticleComboStep0->Get_FinalParticle_Measured(5);
	const DKinematicData* locPhoton4 = locUseKinFitResultsFlag ? locParticleComboStep0->Get_FinalParticle(6) : locParticleComboStep0->Get_FinalParticle_Measured(6);
	//const DKinematicData* locPiMinus = locUseKinFitResultsFlag ? locParticleComboStep->Get_FinalParticle(2) : locParticleComboStep->Get_FinalParticle_Measured(2);

	//Calculate the 4-Vector
	DLorentzVector locP4Sum = locPiPlus->lorentzMomentum() + locPiMinus->lorentzMomentum() + locPhoton1->lorentzMomentum() + locPhoton2->lorentzMomentum() + locPhoton3->lorentzMomentum() + locPhoton4->lorentzMomentum();

	//Optional: Fill histograms
	japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
	{
		// Fill any histograms here
		dInvariantMassHist->Fill(locP4Sum.M());
	}
	japp->RootUnLock(); //RELEASE ROOT LOCK!!

	return true; //return false if you want to use this action to apply a cut
}