// $Id$
//
//    File: DCustomAction_Pi0Cuts.cc
// Created: Thu Jan 22 11:19:47 EST 2015
// Creator: jrsteven (on Linux ifarm1401 2.6.32-431.el6.x86_64 x86_64)
//

#include "DCustomAction_Pi0Cuts.h"

void DCustomAction_Pi0Cuts::Initialize(JEventLoop* locEventLoop)
{

	japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
        {
                //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();
		
                if(gDirectory->Get("Pair_E_DeltaT") == NULL) //check to see if already created by another thread
                        dPair_E_DeltaT = new TH2D("Pair_E_DeltaT", "Photon pair: Minimum Energy vs #Delta t; #Delta t; Min. Energy", 200, -10, 10, 200, 0., 10.);
                else //already created by another thread
                        dPair_E_DeltaT = static_cast<TH2D*>(gDirectory->Get("Pair_E_DeltaT"));
		
		if(gDirectory->Get("Pair_M_DeltaT") == NULL) //check to see if already created by another thread
                        dPair_M_DeltaT = new TH2D("Pair_M_DeltaT", "Photon pair: Mass vs #Delta t; #Delta t; Mass", 200, -10, 10, 200, 0., 1.);
                else //already created by another thread
                        dPair_M_DeltaT = static_cast<TH2D*>(gDirectory->Get("Pair_M_DeltaT"));
        }
        japp->RootUnLock(); //RELEASE ROOT LOCK!!
	
}

bool DCustomAction_Pi0Cuts::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
{
	const DDetectorMatches* locDetectorMatches = NULL;
        locEventLoop->GetSingle(locDetectorMatches);

	const DParticleComboStep* locParticleComboStep = locParticleCombo->Get_ParticleComboStep(2);
        if(locParticleComboStep->Get_InitialParticleID() != Pi0)
                return false;
	
	// get final state particles
        deque<const DKinematicData*> locParticles;
        if(!Get_UseKinFitResultsFlag()) //measured
                locParticleComboStep->Get_FinalParticles_Measured(locParticles);
        else
                locParticleComboStep->Get_FinalParticles(locParticles);

	double locDeltaT = locParticles[0]->time() - locParticles[1]->time();
	double energy1 = locParticles[0]->energy();
	double energy2 = locParticles[1]->energy();
	double minEnergy = TMath::Min(energy1, energy2);
	double radius[2] = {0.,0.};
	int nFCAL = 0;

	DLorentzVector pairP4 = locParticles[0]->lorentzMomentum();
	pairP4 += locParticles[1]->lorentzMomentum();
	
        // loop final state particles
        for(size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i) {
                if(locParticleComboStep->Is_FinalParticleMissing(loc_i)) continue;
		
		// get shower object
		const DNeutralShower* locNeutralShower = static_cast<const DNeutralShower*>(locParticleComboStep->Get_FinalParticle_SourceObject(loc_i));
		if(locNeutralShower == NULL) continue;
		
		// reject showers with matching tracks and count # of FCAL photons
		if(locNeutralShower->dDetectorSystem == SYS_FCAL) {
			nFCAL++;
			const DFCALShower* locFCALShower; 
			locNeutralShower->GetSingle(locFCALShower);
			if(locDetectorMatches->Get_IsMatchedToTrack(locFCALShower)) return false;
		}
		else {
			const DBCALShower* locBCALShower;
			locNeutralShower->GetSingle(locBCALShower);
			if(locDetectorMatches->Get_IsMatchedToTrack(locBCALShower)) return false;
		}
		
		// use to reject FCAL showers near the beamline
		radius[loc_i] = locNeutralShower->dSpacetimeVertex.Perp();
	}

	// keep pi0 if 1 or 2 photons in FCAL (specify parameter in DReaction setup)
	if(nFCAL != dMinFCAL) 
		return false;

	// reject inner FCAL blocks
	if(radius[0] < 20. || radius[1] < 20.)
		return false;

	dPair_E_DeltaT->Fill(locDeltaT, minEnergy);

	// reject pair with low energy shower
	if(minEnergy < 0.5) // was 0.5 for both photons in FCAL
		return false;

	dPair_M_DeltaT->Fill(locDeltaT, pairP4.M());

	// photon cuts
	if(fabs(locDeltaT) < 10) // was 10 for both photons in FCAL
		return true;	

	// reject photon pair if don't pass deltaT and energy cuts
	return false; 
}