// $Id$
//
//    File: DEventProcessor_omega2pi.cc
// Created: Tue May 20 13:06:00 EDT 2014
// Creator: mstaib (on Linux max.phys.cmu.edu 2.6.18-371.3.1.el5 x86_64)
//

#include "DEventProcessor_omega2pi.h"

// Routine used to create our DEventProcessor

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

//------------------
// init
//------------------
jerror_t DEventProcessor_omega2pi::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();
	//  ... create historgrams or trees ...
	if(gDirectory->Get("hNumCorrectTopology") == NULL) 
	hNumCorrectTopology = new TH1I("hNumCorrectTopology","1 = 5Pi, 2 = 3PiP, 3 = 3PiN, 4 = 4Pi", 4, 0.5, 4.5);
	 gDirectory->cd("..");

	 japp->RootUnLock();
	//

	return NOERROR;
}

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

	/*
	//Recommeded: Initialize reaction-independent analysis actions (nothing should be done if already initialized)
		//Warning: If running JANA with multiple plugins, these actions should only be executed in ONE of the plugins.
			//Else you will double/triple/etc.-fill your histograms!!
	dHistogramAction_TrackMultiplicity.Initialize(locEventLoop);
	dHistogramAction_ThrownParticleKinematics.Initialize(locEventLoop);
	dHistogramAction_DetectedParticleKinematics.Initialize(locEventLoop);
	dHistogramAction_GenReconTrackComparison.Initialize(locEventLoop);
	*/

	
	//Recommeded: Create output ROOT TTrees (nothing is done if already created)
	const DEventWriterROOT* locEventWriterROOT = NULL;
	locEventLoop->GetSingle(locEventWriterROOT);
	locEventWriterROOT->Create_DataTrees(locEventLoop);
	

	return NOERROR;
}

//------------------
// evnt
//------------------
jerror_t DEventProcessor_omega2pi::evnt(jana::JEventLoop* locEventLoop, int 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);
	
	vector<const DMCThrown*> locFinalStateParticles;
        vector<const DMCThrown*> locDecayingParticles;

        locEventLoop->Get(locFinalStateParticles, "FinalState");
        locEventLoop->Get(locDecayingParticles, "Decaying");

        Int_t numPiPlus = 0, numPiMinus = 0, numPi0 = 0, numKPlus = 0 , numKMinus = 0, numK0 = 0, numGamma = 0, numProton = 0, numNeutron = 0;
	
        unsigned int numFinalStateParticles = locFinalStateParticles.size();

        // Now that we have the right number, count them.
        for (unsigned int i = 0; i < locFinalStateParticles.size(); i++){
	  const DMCThrown* locMCThrown = locFinalStateParticles[i];
	  int pdgType = locMCThrown->pdgtype;
	  if (pdgType == 211) numPiPlus++;
	  else if (pdgType == -211) numPiMinus++;
	  else if (pdgType == 2212) numProton++;
	  else if (pdgType == 22) numGamma++;
	  else if (pdgType == 321) numKPlus++;
	  else if (pdgType == -321) numKMinus++;
	  else if (pdgType == 130) numK0++;
	  else if (pdgType == 2112) numNeutron++;
        }
        // This still leaves an ambiguity about whether the photons actually came from a pi0 so lets count those

        for (unsigned int i = 0; i < locDecayingParticles.size(); i++){
	  const DMCThrown* locMCThrown = locDecayingParticles[i];
	  int pdgType = locMCThrown->pdgtype;
	  if (pdgType == 111) numPi0++;
        }

        //Now to check if we have what we want and retun
        if (numPiPlus == 2 && numPiMinus == 2 && numProton == 1 && numGamma == 2 && numPi0 == 1 && numFinalStateParticles == 7){
	  japp->RootWriteLock(); 
	  {
	  hNumCorrectTopology->Fill(1);
          }
	  japp->RootUnLock();
	}

	else if (numPiPlus == 1 && numPiMinus == 1 && numProton == 1 && numGamma == 2 && numPi0 == 1 && numFinalStateParticles == 5){
          japp->RootWriteLock();
          {
	  hNumCorrectTopology->Fill(2);
          }
	  japp->RootUnLock();
        }

	else if (numPiPlus == 2 && numPiMinus == 1 && numNeutron == 1 && numFinalStateParticles == 4){
          japp->RootWriteLock();
          {
	  hNumCorrectTopology->Fill(3);
          }
	  japp->RootUnLock();
        }
	else if (numPiPlus == 2 && numPiMinus == 2 && numProton == 1 && numFinalStateParticles == 5){
          japp->RootWriteLock();
          {
	  hNumCorrectTopology->Fill(4);
          }
	  japp->RootUnLock();
        }

	//japp->RootWriteLock();
	  //... fill historgrams or trees ...
	//japp->RootUnLock();

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

	/*
	//Recommended: Execute reaction-independent actions (fill histograms)
		//Warning: If running JANA with multiple plugins, these actions should only be executed in ONE of the plugins.
			//Else you will double/triple/etc.-fill your histograms!!
	dHistogramAction_TrackMultiplicity(locEventLoop);
	dHistogramAction_ThrownParticleKinematics(locEventLoop);
	dHistogramAction_DetectedParticleKinematics(locEventLoop);
	dHistogramAction_GenReconTrackComparison(locEventLoop);
	*/

	
	//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. 
	const DEventWriterROOT* locEventWriterROOT = NULL;
	locEventLoop->GetSingle(locEventWriterROOT);
	// string is DReaction factory tag: will fill trees for all DReactions that are defined in the specified factory
	//locEventWriterROOT->Fill_DataTrees(locEventLoop, "omega2pi");
	

	/*
	//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() != "omega2pi")
			continue; // particle combination was for a different reaction
		//perform further analysis steps here...
	}
	*/

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

	Bool_t isBackground = false;
	
	for(size_t loc_i = 0; loc_i < locAnalysisResultsVector.size(); ++loc_i){
	
		const DAnalysisResults* locAnalysisResults = locAnalysisResultsVector[loc_i];
		if(locAnalysisResults->Get_Reaction()->Get_ReactionName() == "omega2pi" && !isBackground){
			deque<const DParticleCombo*> locPassedParticleCombos;
                	locAnalysisResults->Get_PassedParticleCombos(locPassedParticleCombos);
   			if ( locPassedParticleCombos.size() == 0 ) isBackground = true;
			deque< const DParticleCombo * > bestPassedParticleCombo;
			Float_t bestKinFitFOM = -1;
                	for(size_t loc_j = 0; loc_j < locPassedParticleCombos.size(); ++loc_j)
                	{
                        	const DParticleCombo* locPassedParticleCombo = locPassedParticleCombos[loc_j];
                        	Float_t kinFitFOM = locPassedParticleCombo->Get_KinFitResults()->Get_ConfidenceLevel();
				if (kinFitFOM > bestKinFitFOM){
					bestKinFitFOM = kinFitFOM;
					bestPassedParticleCombo.clear(); 
					bestPassedParticleCombo.push_back(locPassedParticleCombo);	
				}
                       }
		       if (bestKinFitFOM > 0.) locEventWriterROOT->Fill_DataTree(locEventLoop, locAnalysisResults->Get_Reaction(), bestPassedParticleCombo);
                }
		else if (locAnalysisResults->Get_Reaction()->Get_ReactionName() == "omega2piBackground" && isBackground){
			deque<const DParticleCombo*> locPassedParticleCombos;
                        locAnalysisResults->Get_PassedParticleCombos(locPassedParticleCombos);
			deque< const DParticleCombo * > bestPassedParticleCombo;
                        Float_t bestKinFitFOM = -1;
                        for(size_t loc_j = 0; loc_j < locPassedParticleCombos.size(); ++loc_j)
                        {
                                const DParticleCombo* locPassedParticleCombo = locPassedParticleCombos[loc_j];
                                Float_t kinFitFOM = locPassedParticleCombo->Get_KinFitResults()->Get_ConfidenceLevel();
                                if (kinFitFOM > bestKinFitFOM){
                                        bestKinFitFOM = kinFitFOM;
                                        bestPassedParticleCombo.clear();
                                        bestPassedParticleCombo.push_back(locPassedParticleCombo);
                                }
                       }
                       if (bestKinFitFOM > 0.) locEventWriterROOT->Fill_DataTree(locEventLoop, locAnalysisResults->Get_Reaction(), bestPassedParticleCombo);

                 }
	}
	

	/*
	//Optional: Save event to output REST file. Use this to create skims.
	const DEventWriterREST* locEventWriterREST = NULL;
	locEventLoop->GetSingle(locEventWriterREST);
	locEventWriterREST->Write_RESTEvent(locEventLoop, "omega2pi"); //string is part of output file name
	*/

	return NOERROR;
}

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