// $Id$
//
//    File: DEventProcessor_angle_test.cc
// Created: Thu Mar 23 13:54:41 EDT 2017
// Creator: pmatt (on Linux pmatt-VirtualBox 4.4.0-31-generic x86_64)
//

#include "DEventProcessor_angle_test.h"

// Routine used to create our DEventProcessor

extern "C"
{
	void InitPlugin(JApplication *locApplication)
	{
		InitJANAPlugin(locApplication);
		locApplication->AddProcessor(new DEventProcessor_angle_test()); //register this plugin
		locApplication->AddFactoryGenerator(new DFactoryGenerator_angle_test()); //register the factory generator
	}
} // "C"

//------------------
// init
//------------------
jerror_t DEventProcessor_angle_test::init(void)
{
	// This is called once at program startup.

	// This is called once at program startup.
	dDeltaMassVsMass = new TH2D("DeltaMassVsMass", ";#gamma#gamma mass (GeV/c^{2});#Deltam (GeV/c^{2})", 2000, 0.0, 2.0, 2000, -1.0, 1.0);
	dDeltaMassVsDeltaZ_AllMasses = new TH2D("DeltaMassVsDeltaZ_AllMasses", "All Masses;#Deltaz (cm);#Deltam (GeV/c^{2})", 1000, -50.0, 50.0, 2000, -1.0, 1.0);
	dDeltaMassVsDeltaZ_NearPi0 = new TH2D("DeltaMassVsDeltaZ_NearPi0", "Near #pi^{0};#Deltaz (cm);#Deltam (GeV/c^{2})", 1000, -50.0, 50.0, 2000, -1.0, 1.0);
	dDeltaMassVsDeltaZ_NearEta = new TH2D("DeltaMassVsDeltaZ_NearEta", "Near #eta;#Deltaz (cm);#Deltam (GeV/c^{2})", 1000, -50.0, 50.0, 2000, -1.0, 1.0);

	/*
	//OPTIONAL: Create an EventStore skim.  
	string locSkimFileName = "angle_test.idxa";
	dEventStoreSkimStream.open(locSkimFileName.c_str());
	dEventStoreSkimStream << "IDXA" << endl;
	*/

	return NOERROR;
}

//------------------
// brun
//------------------
jerror_t DEventProcessor_angle_test::brun(jana::JEventLoop* locEventLoop, int32_t locRunNumber)
{
	// This is called whenever the run number changes

	return NOERROR;
}

//------------------
// evnt
//------------------
jerror_t DEventProcessor_angle_test::evnt(jana::JEventLoop* locEventLoop, uint64_t locEventNumber)
{
	// This is called for every event. Use of common resources like writing
	// to a file or filling a histogram should be mutex protected. Using
	// locEventLoop->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;
	// locEventLoop->Get(mydataclasses);
	//
	// japp->RootFillLock(this);
	//  ... fill historgrams or trees ...
	// japp->RootFillUnLock(this);

	vector<const DNeutralParticle*> locParticles;
	locEventLoop->Get(locParticles, "PreSelect");

	const DVertex* locVertex;
	locEventLoop->GetSingle(locVertex);

	//make pairs of particles
	vector<vector<const DNeutralParticle*> > locNeutralPairs;
	for(auto locIterator = locParticles.begin(); locIterator != locParticles.end(); ++locIterator)
	{
		for(auto locIterator2 = std::next(locIterator); locIterator2 != locParticles.end(); ++locIterator2)
			locNeutralPairs.emplace_back(vector<const DNeutralParticle*>{*locIterator, *locIterator2});
	}

	//calculate invariant masses of pairs using the given DVertex (default p4s)
	map<vector<const DNeutralParticle*>, double> locPairMasses;
	for(auto locPair : locNeutralPairs)
	{
		auto locHypo1 = locPair[0]->Get_Hypothesis(Gamma);
		auto locHypo2 = locPair[1]->Get_Hypothesis(Gamma);
		double locMass = (locHypo1->lorentzMomentum() + locHypo2->lorentzMomentum()).M();
		locPairMasses.emplace(locPair, locMass);
	}

	//now calculate invariant masses of pairs using the center of the target
	map<vector<const DNeutralParticle*>, double> locPairMasses_TargetCenter;
	TVector3 locTargetCenter(0.0, 0.0, 65.0);
	for(auto locPair : locNeutralPairs)
	{
		auto& locShower1 = locPair[0]->dNeutralShower;
		TVector3 locMomentum1 = locShower1->dSpacetimeVertex.Vect() - locTargetCenter;
		locMomentum1.SetMag(locShower1->dEnergy);
		TLorentzVector locP4_1(locMomentum1, locShower1->dEnergy);

		auto& locShower2 = locPair[1]->dNeutralShower;
		TVector3 locMomentum2 = locShower2->dSpacetimeVertex.Vect() - locTargetCenter;
		locMomentum2.SetMag(locShower2->dEnergy);
		TLorentzVector locP4_2(locMomentum2, locShower2->dEnergy);

		double locMass = (locP4_1 + locP4_2).M();
		locPairMasses_TargetCenter.emplace(locPair, locMass);
	}

	//study mass difference as a function of delta-z
	//at all masses, near pi0 mass, near eta mass
	double locDeltaZ = locVertex->dSpacetimeVertex.Z() - 65.0;
	vector<pair<double, double> > locMassDeltas; //first is mass_dvertex, second is delta mass
	for(auto& locMassPair : locPairMasses)
	{
		auto& locNeutralPair = locMassPair.first;
		double locMass_DVertex = locMassPair.second;
		double locMass_Target = locPairMasses_TargetCenter[locNeutralPair];
		double locDeltaMass = locMass_DVertex - locMass_Target;
		locMassDeltas.emplace_back(locMass_DVertex, locDeltaMass);
	}

	//histogram masses
	japp->RootFillLock(this);
	{
		for(auto& locMassPair : locMassDeltas)
		{
			double& locMass_DVertex = locMassPair.first;
			double& locDeltaMass = locMassPair.second;

			dDeltaMassVsMass->Fill(locMass_DVertex, locDeltaMass);
			dDeltaMassVsDeltaZ_AllMasses->Fill(locDeltaZ, locDeltaMass);
			if((locMass_DVertex >= 0.05) && (locMass_DVertex <= 0.22))
				dDeltaMassVsDeltaZ_NearPi0->Fill(locDeltaZ, locDeltaMass);
			if((locMass_DVertex >= 0.5) && (locMass_DVertex <= 0.6))
				dDeltaMassVsDeltaZ_NearEta->Fill(locDeltaZ, locDeltaMass);
		}
	}
	japp->RootFillUnLock(this);

	// DOCUMENTATION:
	// ANALYSIS library: https://halldweb1.jlab.org/wiki/index.php/GlueX_Analysis_Software

	/*********************************************************** REQUIRED ***********************************************************/

	//REQUIRED: To run an analysis, You MUST call one at least of the below code fragments. 
		//JANA is on-demand, so if you don't call one of these, then your analysis won't run. 

	/*
	//Recommended: Write surviving particle combinations (if any) to output ROOT TTree
		//If no cuts are performed by the analysis actions added to a DReaction, then this saves all of its particle combinations. 
		//The event writer gets the DAnalysisResults objects from JANA, performing the analysis. 
	// string is DReaction factory tag: will fill trees for all DReactions that are defined in the specified factory
	const DEventWriterROOT* locEventWriterROOT = NULL;
	locEventLoop->GetSingle(locEventWriterROOT);
	locEventWriterROOT->Fill_DataTrees(locEventLoop, "angle_test");
	*/

	/*
	//Optional: Get the analysis results for all DReactions. 
		//Getting these objects triggers the analysis, if it wasn't performed already. 
		//These objects contain the DParticleCombo objects that survived the DAnalysisAction cuts that were added to the DReactions
	vector<const DAnalysisResults*> locAnalysisResultsVector;
	locEventLoop->Get(locAnalysisResultsVector);
	*/

	/************************************************** OPTIONAL: FURTHER ANALYSIS **************************************************/
	/*
	//Optional: Get all particle combinations for all DReactions. 
		//If kinematic fits were requested, these contain both the measured and kinematic-fit track parameters
		//No cuts from DAnalysisActions are placed on these combos
	vector<const DParticleCombo*> locParticleCombos;
	locEventLoop->Get(locParticleCombos);
	for(size_t loc_i = 0; loc_i < locParticleCombos.size(); ++loc_i)
	{
		const DParticleCombo* locParticleCombo = locParticleCombos[loc_i];
		if(locParticleCombo->Get_Reaction()->Get_ReactionName() != "angle_test")
			continue; // particle combination was for a different reaction
		//perform further analysis steps here...
	}
	*/

	/*
	//Optional: Perform further cuts on the particle combos in the analysis results. 
	for(size_t loc_i = 0; loc_i < locAnalysisResultsVector.size(); ++loc_i)
	{
		const DAnalysisResults* locAnalysisResults = locAnalysisResultsVector[loc_i];
		if(locAnalysisResults->Get_Reaction()->Get_ReactionName() != "angle_test")
			continue; // analysis results were for a different reaction
		//get the DParticleCombo objects for this DReaction that survived all of the DAnalysisAction cuts
		deque<const DParticleCombo*> locPassedParticleCombos;
		locAnalysisResults->Get_PassedParticleCombos(locPassedParticleCombos);
		for(size_t loc_j = 0; loc_j < locPassedParticleCombos.size(); ++loc_j)
		{
			const DParticleCombo* locPassedParticleCombo = locPassedParticleCombos[loc_j];
			//perform further analysis steps here...
		}
	}
	*/

	/******************************************************** OPTIONAL: SKIMS *******************************************************/

	/*
	//Optional: Save event to output REST file. Use this to create physical skims.
	const DEventWriterREST* locEventWriterREST = NULL;
	locEventLoop->GetSingle(locEventWriterREST);
	for(size_t loc_i = 0; loc_i < locAnalysisResultsVector.size(); ++loc_i)
	{
		const DAnalysisResults* locAnalysisResults = locAnalysisResultsVector[loc_i];
		if(locAnalysisResults->Get_Reaction()->Get_ReactionName() != "angle_test")
			continue; // analysis results were for a different reaction

		//get the DParticleCombo objects for this DReaction that survived all of the DAnalysisAction cuts
		deque<const DParticleCombo*> locPassedParticleCombos;
		locAnalysisResults->Get_PassedParticleCombos(locPassedParticleCombos);

		if(!locPassedParticleCombos.empty())
			locEventWriterREST->Write_RESTEvent(locEventLoop, "angle_test"); //string is part of output file name
	}
	*/

	/*
	//Optional: Create an EventStore skim. 

	// See whether this is MC data or real data
	vector<const DMCThrown*> locMCThrowns;
	locEventLoop->Get(locMCThrowns);

	unsigned int locRunNumber = locEventLoop->GetJEvent().GetRunNumber();
	unsigned int locUniqueID = locMCThrowns.empty() ? 1 : Get_FileNumber(locEventLoop);

	// If a particle combo passed the cuts, save the event info in the output file
	for(size_t loc_i = 0; loc_i < locAnalysisResultsVector.size(); ++loc_i)
	{
		const DAnalysisResults* locAnalysisResults = locAnalysisResultsVector[loc_i];
		if(locAnalysisResults->Get_Reaction()->Get_ReactionName() != "angle_test")
			continue; // analysis results were for a different reaction
		if(locAnalysisResults->Get_NumPassedParticleCombos() == 0)
			continue; // no combos passed

		//MUST LOCK AROUND MODIFICATION OF MEMBER VARIABLES IN brun() or evnt().
		japp->WriteLock("angle_test.idxa"); //Lock is unique to this output file
		{
			dEventStoreSkimStream << locRunNumber << " " << locEventNumber << " " << locUniqueID << endl;
		}
		japp->Unlock("angle_test.idxa");
	}
	*/

	return NOERROR;
}

int DEventProcessor_angle_test::Get_FileNumber(JEventLoop* locEventLoop) const
{
	//Assume that the file name is in the format: *_X.ext, where:
		//X is the file number (a string of numbers of any length)
		//ext is the file extension (probably .evio or .hddm)

	//get the event source
	JEventSource* locEventSource = locEventLoop->GetJEvent().GetJEventSource();
	if(locEventSource == NULL)
		return -1;

	//get the source file name (strip the path)
	string locSourceFileName = locEventSource->GetSourceName();

	//find the last "_" & "." indices
	size_t locUnderscoreIndex = locSourceFileName.rfind("_");
	size_t locDotIndex = locSourceFileName.rfind(".");
	if((locUnderscoreIndex == string::npos) || (locDotIndex == string::npos))
		return -1;

	size_t locNumberLength = locDotIndex - locUnderscoreIndex - 1;
	string locFileNumberString = locSourceFileName.substr(locUnderscoreIndex + 1, locNumberLength);

	int locFileNumber = -1;
	istringstream locFileNumberStream(locFileNumberString);
	locFileNumberStream >> locFileNumber;

	return locFileNumber;
}

//------------------
// erun
//------------------
jerror_t DEventProcessor_angle_test::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 DEventProcessor_angle_test::fini(void)
{
	// Called before program exit after event processing is finished.
	if(dEventStoreSkimStream.is_open())
		dEventStoreSkimStream.close();
	return NOERROR;
}