// $Id$
//
//    File: DCustomAction_BackgroundSelector.cc
// Created: Mon Jun 16 09:56:25 EDT 2014
// Creator: mstaib (on Linux max.phys.cmu.edu 2.6.18-371.3.1.el5 x86_64)
//

#include "DCustomAction_BackgroundSelector.h"

void DCustomAction_BackgroundSelector::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();

		// Optional: Create a ROOT subfolder.
			//If another thread has already created the folder, it just changes to it. 
		// CreateAndChangeTo_Directory("MyDirName", "MyDirTitle");
			//make sub-directory content here
		// gDirectory->cd(".."); //return to the action directory

		//	(Optional) Example: Create a histogram. 
		//	if(gDirectory->Get("MyHistName") == NULL) //check to see if already created by another thread
		//		dMyHist = new TH1D("MyHistName", "MyHistTitle", 100, 0.0, 1.0);
		//	else //already created by another thread
		//		dMyHist = static_cast<TH1D*>(gDirectory->Get("MyHistName"));
	}
	japp->RootUnLock(); //RELEASE ROOT LOCK!!
	*/
}

bool DCustomAction_BackgroundSelector::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
{
	//Write custom code to perform an action on the INPUT DParticleCombo (DParticleCombo)
	//NEVER: Grab DParticleCombo or DAnalysisResults objects (of any tag!) from the JEventLoop within this function
	//NEVER: Grab objects that are created post-kinfit (e.g. DKinFitResults, etc.) from the JEventLoop if Get_UseKinFitResultsFlag() == false: CAN CAUSE INFINITE DEPENDENCY LOOP


	//Optional: check whether the user wanted to use the kinematic fit results when performing this action
	// bool locUseKinFitResultsFlag = Get_UseKinFitResultsFlag();


	/*
	//Optional: Fill histograms
	japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
	{
		// Fill any histograms here
	}
	japp->RootUnLock(); //RELEASE ROOT LOCK!!
	*/
	int isCorrectTopology = CheckIsCorrectTopology(Get_Reaction(), locEventLoop);
        int isCorrectPID = CheckIsCorrectPID(Get_Reaction(), locEventLoop, locParticleCombo);
	if ( isCorrectTopology == 0 || isCorrectPID == 0 ) return true; 
	return false; //return false if you want to use this action to apply a cut (and it fails the cut!)
}

int DCustomAction_BackgroundSelector::CheckIsCorrectTopology(const DReaction* locReaction, JEventLoop* locEventLoop)
{
  vector<const DMCThrown*> locFinalStateParticles;
  vector<const DMCThrown*> locDecayingParticles;

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

  Int_t numPiPlus = 0, numPiMinus = 0, numProton = 0, numGamma = 0, numPi0 = 0;

  if (locFinalStateParticles.size() != 7) return 0; //Needs to have the exact number of particles in the final state
  // 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++;
  }
  // 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){
    return 1;
  }
  else return 0;
}

int DCustomAction_BackgroundSelector::CheckIsCorrectPID(const DReaction* locReaction, JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
{
    float locMinThrownMatchFOM = 0.0;

    vector<const DMCThrownMatching*> locMCThrownMatchingVector;
    locEventLoop->Get(locMCThrownMatchingVector);
    const DMCThrownMatching* locMCThrownMatching = locMCThrownMatchingVector[0];
    const DMCThrown* locMCThrown;
 
    deque<const DKinematicData*> locParticles;
    locParticleCombo->Get_DetectedFinalParticles_Measured(locParticles);
 
    for(size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i)
    {
       if(ParticleCharge(locParticles[loc_i]->PID()) == 0)
       {
          double locMatchFOM = 0.0;
          const DNeutralParticleHypothesis* locNeutralParticleHypothesis = static_cast<const DNeutralParticleHypothesis*>(locParticles[loc_i]);
          locMCThrown = locMCThrownMatching->Get_MatchingMCThrown(locNeutralParticleHypothesis, locMatchFOM);
          if((locMCThrown == NULL) || (locMatchFOM < locMinThrownMatchFOM))
             return 0;
          if(((Particle_t)locMCThrown->type) != locParticles[loc_i]->PID())
             return 0;
       }
       else
       {
          double locMatchFOM = 0.0;
          const DChargedTrackHypothesis* locChargedTrackHypothesis = static_cast<const DChargedTrackHypothesis*>(locParticles[loc_i]);
          locMCThrown = locMCThrownMatching->Get_MatchingMCThrown(locChargedTrackHypothesis, locMatchFOM);
          if((locMCThrown == NULL) || (locMatchFOM < locMinThrownMatchFOM))
             return 0;
          if(((Particle_t)locMCThrown->type) != locParticles[loc_i]->PID())
             return 0;
       }
    }
    return 1;
}