// $Id$
//
//    File: JEventProcessor_3pip.cc
// Created: Mon Dec  1 20:41:57 EST 2014
// Creator: mstaib (on Linux ifarm1101 2.6.32-220.7.1.el6.x86_64 x86_64)
//

#include "JEventProcessor_3pip.h"
using namespace jana;


// Routine used to create our JEventProcessor
#include <JANA/JApplication.h>
#include <JANA/JFactory.h>
extern "C"{
void InitPlugin(JApplication *app){
	InitJANAPlugin(app);
	app->AddProcessor(new JEventProcessor_3pip());
}
} // "C"


//------------------
// JEventProcessor_3pip (Constructor)
//------------------
JEventProcessor_3pip::JEventProcessor_3pip()
{

}

//------------------
// ~JEventProcessor_3pip (Destructor)
//------------------
JEventProcessor_3pip::~JEventProcessor_3pip()
{

}

//------------------
// init
//------------------
jerror_t JEventProcessor_3pip::init(void)
{
	// This is called once at program startup. If you are creating
	// and filling historgrams in this plugin, you should lock the
	// ROOT mutex like this:
	//
	// japp->RootWriteLock();
	//  ... fill historgrams or trees ...
	// japp->RootUnLock();
	//
	japp->RootWriteLock();
	if (gDirectory->Get("h3PiInvMass") == NULL)
		h3PiInvMass = new TH1I("h3PiInvMass", "Inv. Mass of #pi^{+} #pi^{-} #pi^{0} system; Invariant Mass [GeV]; Entries per 5 MeV", 1000, 0, 5.0);
	if (gDirectory->Get("h2PiInvMass") == NULL)
                h2PiInvMass = new TH1I("h2PiInvMass", "Inv. Mass of #pi^{+} #pi^{-}; Invariant Mass [GeV]; Entries per 5 MeV", 1000, 0, 5);
	if (gDirectory->Get("hPi0InvMass") == NULL)
		hPi0InvMass = new TH1I("hPi0InvMass", "Inv Mass of 2#gamma system; Invariant Mass [GeV]; Entries per 5 MeV", 200, 0, 1);;
	japp->RootUnLock();
	return NOERROR;
}

//------------------
// brun
//------------------
jerror_t JEventProcessor_3pip::brun(JEventLoop *eventLoop, int runnumber)
{
	// This is called whenever the run number changes
	return NOERROR;
}

//------------------
// evnt
//------------------
jerror_t JEventProcessor_3pip::evnt(JEventLoop *loop, int eventnumber)
{
	// This is called for every event. Use of common resources like writing
	// to a file or filling a histogram should be mutex protected. Using
	// loop->Get(...) to get reconstructed objects (and thereby activating the
	// reconstruction algorithm) should be done outside of any mutex lock
	// since multiple threads may call this method at the same time.
	// Here's an example:
	//
	// vector<const MyDataClass*> mydataclasses;
	// loop->Get(mydataclasses);
	//
	// japp->RootWriteLock();
	//  ... fill historgrams or trees ...
	// japp->RootUnLock();
	vector < const DNeutralParticle * > neutralParticleVector;
	vector < const DTrackTimeBased * > timeBasedTrackVector;

	loop->Get(neutralParticleVector);
	loop->Get(timeBasedTrackVector);

	if (neutralParticleVector.size() == 0 || timeBasedTrackVector.size() == 0) return NOERROR;

	vector < const DNeutralParticleHypothesis *> neutralParticleHypothesisPassedCuts;
	vector < const DTrackTimeBased * > timeBasedTrackProton;
	vector < const DTrackTimeBased * > timeBasedTrackPiPlus;
	vector < const DTrackTimeBased * > timeBasedTrackPiMinus;

	//for ( vector < const DNeutralParticle * >::const_iterator iter = neutralParticleVector.begin(); iter != neutralParticleVector.end(); iter++){
	//	const DNeutralParticleHypothesis *thisGamma = (*iter)->Get_Hypothesis(Gamma);
	//	if (thisGamma->energy() > 0.5) neutralParticleHypothesisPassedCuts.push_back(thisGamma);
	//}

	for (vector < const DTrackTimeBased * >::const_iterator iter = timeBasedTrackVector.begin(); iter != timeBasedTrackVector.end(); iter++){
		const DTrackTimeBased * thisTrack = (*iter);
		if (thisTrack->FOM < 10E-25 ) continue;
		if (thisTrack->PID() == Proton) timeBasedTrackProton.push_back(*iter);
		else if (thisTrack->PID() == PiPlus) timeBasedTrackPiPlus.push_back(*iter);
		else if (thisTrack->PID() == PiMinus) timeBasedTrackPiMinus.push_back(*iter);
	}

	// For Now lets just do the ones that have exactly 2 positive and one negative track
	
	if( timeBasedTrackProton.size() != 2 || timeBasedTrackPiPlus.size() != 2 || timeBasedTrackPiMinus.size() != 1) return NOERROR;
		
	// Now just to fill the histograms
	

	// For each combination of photons passing a timing cut, try to calculate the invariant mass
	
	for (unsigned int i = 0 ; i < neutralParticleVector.size(); i++){
		if (neutralParticleVector[i]->dNeutralShower->dEnergy < 0.5) continue;
		double time1 = (neutralParticleVector[i]->dNeutralShower->dSpacetimeVertex).T();

		for (unsigned int j = i+1 ; j < neutralParticleVector.size(); j++){
			if (neutralParticleVector[j]->dNeutralShower->dEnergy < 0.5) continue;
                	double time2 =  (neutralParticleVector[j]->dNeutralShower->dSpacetimeVertex).T();	

			if(fabs(time1 - time2) > 20) continue;

                	const DNeutralParticleHypothesis *gamma1 = neutralParticleVector[i]->Get_Hypothesis(Gamma);
                	const DNeutralParticleHypothesis *gamma2 = neutralParticleVector[j]->Get_Hypothesis(Gamma);

			const DLorentzVector g1 = gamma1->lorentzMomentum();
			const DLorentzVector g2 = gamma2->lorentzMomentum();
			const DLorentzVector P1 = timeBasedTrackProton[0]->lorentzMomentum();
			const DLorentzVector P2 = timeBasedTrackProton[1]->lorentzMomentum();
			const DLorentzVector PiPlus1 = timeBasedTrackPiPlus[0]->lorentzMomentum();
                        const DLorentzVector PiPlus2 = timeBasedTrackPiPlus[1]->lorentzMomentum();
			const DLorentzVector PiMinus1 = timeBasedTrackPiMinus[0]->lorentzMomentum();		
			japp->RootWriteLock();
			h3PiInvMass->Fill((PiPlus2 + PiMinus1 + g1 + g2).M());
			h3PiInvMass->Fill((PiPlus1 + PiMinus1 + g1 + g2).M());
			h2PiInvMass->Fill((PiPlus1 + PiMinus1).M());
			h2PiInvMass->Fill((PiPlus2 + PiMinus1).M());
			hPi0InvMass->Fill((g1 + g2).M());
			japp->RootUnLock();
		}	
	}



	return NOERROR;
}

//------------------
// erun
//------------------
jerror_t JEventProcessor_3pip::erun(void)
{
	// This is called whenever the run number changes, before it is
	// changed to give you a chance to clean up before processing
	// events from the next run number.
	return NOERROR;
}

//------------------
// fini
//------------------
jerror_t JEventProcessor_3pip::fini(void)
{
	// Called before program exit after event processing is finished.
	return NOERROR;
}