// $Id$
//
//    File: DCustomAction_p3pi_hists.cc
// Created: Wed Jan 21 16:53:41 EST 2015
// Creator: jrsteven (on Linux ifarm1401 2.6.32-431.el6.x86_64 x86_64)
//

#include "DCustomAction_p3pi_hists.h"

void DCustomAction_p3pi_hists::Initialize(JEventLoop* locEventLoop)
{
	// get PID algos
        const DParticleID* locParticleID = NULL;
        locEventLoop->GetSingle(locParticleID);
        dParticleID = locParticleID;
	
	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("M3pi") == NULL) //check to see if already created by another thread
			dM3pi = new TH1D("M3pi", "#pi^{+} #pi^{-} #pi^{0} Mass; Mass (GeV/c^{2})", 200, 0.0, 2.0);
		else //already created by another thread
			dM3pi = static_cast<TH1D*>(gDirectory->Get("M3pi"));
		
		if(gDirectory->Get("MM_M3pi") == NULL) //check to see if already created by another thread
			dMM_M3pi = new TH2D("MM_M3pi", "MM off #pi^{+}#pi^{-}#pi^{0} vs M_{#pi^{+}#pi^{-}#pi^{0}}; M_{#pi^{+}#pi^{-}#pi^{0}}; MM", 200, 0.0, 2.0, 200, 0., 4.);
		else //already created by another thread
			dMM_M3pi = static_cast<TH2D*>(gDirectory->Get("MM_M3pi"));

		if(gDirectory->Get("DeltaT3pi") == NULL) //check to see if already created by another thread
                        dDeltaT3pi = new TH2D("DeltaT3pi", "Match SC-TAGGER #Delta t #pi_{1} vs #Delta t #pi_{1}; #Delta t #pi_{2}; #Delta t #pi_{1}", 100, -50.0, 50.0, 100, -50., 50.);
                else //already created by another thread
                        dDeltaT3pi = static_cast<TH2D*>(gDirectory->Get("DeltaT3pi"));

                if(gDirectory->Get("Egamma3pi") == NULL) //check to see if already created by another thread
                        dEgamma3pi = new TH1D("Egamma3pi", "TAGGER photon energy; E_{gamma}", 320, 2., 10.);
                else //already created by another thread
                        dEgamma3pi = static_cast<TH1D*>(gDirectory->Get("Egamma3pi"));

		if(gDirectory->Get("MissingDeltaPT_DeltaPhi") == NULL) //check to see if already created by another thread
                        dMissingDeltaPT_DeltaPhi = new TH2D("MissingDeltaPT_DeltaPhi", "Unused Track - MM: #Delta p_{T} vs #Delta #phi; #Delta #phi; #Delta p_{T}", 180, -60.0, 60.0, 200, -0.5, 0.5);
                else //already created by another thread
                        dMissingDeltaPT_DeltaPhi = static_cast<TH2D*>(gDirectory->Get("MissingDeltaPT_DeltaPhi"));	

		if(gDirectory->Get("MM2_DeltaPhi") == NULL) //check to see if already created by another thread
                        dMM2_DeltaPhi = new TH2D("MM2_DeltaPhi", "MM^{2} vs #Delta #phi; #Delta #phi; MM^{2}", 180, -180.0, 180.0, 100, -2., 2.);
                else //already created by another thread
                        dMM2_DeltaPhi = static_cast<TH2D*>(gDirectory->Get("MM2_DeltaPhi"));	

		if(gDirectory->Get("MM_M3pi_2SC") == NULL) //check to see if already created by another thread
			dMM_M3pi_2SC = new TH2D("MM_M3pi_2SC", "MM off #pi^{+}#pi^{-}#pi^{0} vs M_{#pi^{+}#pi^{-}#pi^{0}}; M_{#pi^{+}#pi^{-}#pi^{0}}; MM", 200, 0.0, 2.0, 200, 0., 4.);
		else //already created by another thread
			dMM_M3pi_2SC = static_cast<TH2D*>(gDirectory->Get("MM_M3pi_2SC"));

		if(gDirectory->Get("MM_M3pi_Emeson") == NULL) //check to see if already created by another thread
			dMM_M3pi_Emeson = new TH2D("MM_M3pi_Emeson", "MM off #pi^{+}#pi^{-}#pi^{0} vs M_{#pi^{+}#pi^{-}#pi^{0}}; M_{#pi^{+}#pi^{-}#pi^{0}}; MM", 200, 0.0, 2.0, 200, 0., 4.);
		else //already created by another thread
			dMM_M3pi_Emeson = static_cast<TH2D*>(gDirectory->Get("MM_M3pi_Emeson"));

		if(gDirectory->Get("CosThetaCM_t") == NULL) //check to see if already created by another thread
			dCosThetaCM_t = new TH2D("CosThetaCM_t", "Mandelstam -t vs cos#theta_{CM}^{#omega}; cos#theta_{CM}^{#omega}; -t", 1000, -1., 1., 1000, 0., 5.);
		else //already created by another thread
			dCosThetaCM_t = static_cast<TH2D*>(gDirectory->Get("CosThetaCM_t"));
	}
	japp->RootUnLock(); //RELEASE ROOT LOCK!!
}

bool DCustomAction_p3pi_hists::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
{
	const DDetectorMatches* locDetectorMatches = NULL;
        locEventLoop->GetSingle(locDetectorMatches);
	
	// for unused tracks analysis
	vector<const DChargedTrack*> locChargedTracks;
        locEventLoop->Get(locChargedTracks, "PreSelect");
	vector<const DChargedTrackHypothesis*> locUnusedChargedTrackHypotheses;
        for(size_t loc_i = 0; loc_i < locChargedTracks.size(); ++loc_i)
                locUnusedChargedTrackHypotheses.insert(locUnusedChargedTrackHypotheses.end(), locChargedTracks[loc_i]->dChargedTrackHypotheses.begin(), locChargedTracks[loc_i]->dChargedTrackHypotheses.end());

	const DParticleComboStep* locParticleComboStep = locParticleCombo->Get_ParticleComboStep(1);
        if(locParticleComboStep->Get_InitialParticleID() != omega)
                return true;
	
        // get beam photon time
        const DKinematicData* locBeamPhoton = locParticleCombo->Get_ParticleComboStep(0)->Get_InitialParticle();
        double locBeamPhotonTime = locBeamPhoton->time();
        double locBeamPhotonEnergy = locBeamPhoton->energy();
	//if(locBeamPhotonEnergy < 6.8) return true; // cut low energy photons

	// get final state particles
        deque<const DKinematicData*> locParticles;
        if(!Get_UseKinFitResultsFlag()) //measured
                locParticleComboStep->Get_FinalParticles_Measured(locParticles);
        else
                locParticleComboStep->Get_FinalParticles(locParticles);
	
        // counter for deltaT
        int nFoundSC = 0;
        double locDeltaT[2] = {-999.,-999};
	
        // calculate time difference
        for(size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i) {
                if(locParticleComboStep->Get_FinalParticleID(loc_i) != PiPlus && locParticleComboStep->Get_FinalParticleID(loc_i) != PiMinus) continue;
		
		// get time based track
                const DChargedTrack* locChargedTrack = static_cast<const DChargedTrack*>(locParticleComboStep->Get_FinalParticle_SourceObject(loc_i));
                if(locChargedTrack == NULL) continue; // should never happen
		
                const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTrack->Get_BestFOM();
                const DTrackTimeBased* locTrackTimeBased = NULL;
                locChargedTrackHypothesis->GetSingleT(locTrackTimeBased);
		
                // get match to SC
                DSCHitMatchParams locSCHitMatchParams;
                bool foundSC = dParticleID->Get_BestSCMatchParams(locTrackTimeBased, locDetectorMatches, locSCHitMatchParams);
                if(foundSC){
                        double locSCTime = locSCHitMatchParams.dHitTime - locSCHitMatchParams.dFlightTime;
                        locDeltaT[nFoundSC] = locSCTime - locBeamPhotonTime;
                        nFoundSC++;
                }
		
		// elminate pi+/pi- from unused tracks
		for(size_t loc_j = 0; loc_j < locUnusedChargedTrackHypotheses.size(); ++loc_j) {
			if(locUnusedChargedTrackHypotheses[loc_j]->candidateid != locChargedTrackHypothesis->candidateid)
				continue;
			locUnusedChargedTrackHypotheses.erase(locUnusedChargedTrackHypotheses.begin() + loc_j);
		}
	}

	vector<const DMCThrown*> locMCThrown;
	locEventLoop->Get(locMCThrown);

	// calculate missing mass and 3pi mass
	//DLorentzVector locMissingP4 = dAnalysisUtilities->Calc_MissingP4(locParticleCombo, Get_UseKinFitResultsFlag());
	DLorentzVector locMissingP4(0,0,locBeamPhotonEnergy,locBeamPhotonEnergy+0.938);
	//DLorentzVector locOmegaP4 = dAnalysisUtilities->Calc_FinalStateP4(locParticleCombo, 1, Get_UseKinFitResultsFlag());
	//locMissingP4 -= locOmegaP4;
	DLorentzVector locOmegaP4;

	for(size_t loc_i = 0; loc_i < 2; ++loc_i) {
		//if(!Get_UseKinFitResultsFlag())
		//cout<<loc_i<<" "<<locParticles[loc_i]->PID()<<" p="<<locParticles[loc_i]->momentum().Mag()<<endl;
		locMissingP4 -= locParticles[loc_i]->lorentzMomentum();
		locOmegaP4 += locParticles[loc_i]->lorentzMomentum();
	}

	const DParticleComboStep* locParticleComboStepPi0 = locParticleCombo->Get_ParticleComboStep(2);
	deque<const DKinematicData*> locParticlesPi0;
        if(!Get_UseKinFitResultsFlag()) //measured
                locParticleComboStepPi0->Get_FinalParticles_Measured(locParticlesPi0);
        else
                locParticleComboStepPi0->Get_FinalParticles(locParticlesPi0);
	for(size_t loc_i = 0; loc_i < 2; ++loc_i) {
		//if(!Get_UseKinFitResultsFlag()) 
		//cout<<loc_i<<" "<<locParticlesPi0[loc_i]->PID()<<" p="<<locParticlesPi0[loc_i]->momentum().Mag()<<endl;
		DLorentzVector locPhotonP4 = locParticlesPi0[loc_i]->lorentzMomentum();
		if(!locMCThrown.empty()) locPhotonP4 *= 1.5;
		locMissingP4 -= locPhotonP4;
		locOmegaP4 += locPhotonP4;
	}

	// check MC tracks
	/*
	if(!locMCThrown.empty()){
		if(!Get_UseKinFitResultsFlag()) {
			cout<<"Found good event "<<locEventLoop->GetJEvent().GetEventNumber()<<": M_3pi="<<locOmegaP4.M()<<" MM="<<locMissingP4.M()<<endl;
			for(uint i=0; i<locMCThrown.size(); i++){
				cout<<locMCThrown[i]->type<<" "<<locMCThrown[i]->momentum().Mag()<<endl;
			}
		}
	}
	*/

	japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
	{
		// Fill any histograms here
		dM3pi->Fill(locOmegaP4.M());
		dMM_M3pi->Fill(locOmegaP4.M(), locMissingP4.M());
		dEgamma3pi->Fill(locBeamPhotonEnergy);

		// 3pi energy cut (Dave M. suggestion)
		if(locOmegaP4.E() > 5) {
			dMM_M3pi_Emeson->Fill(locOmegaP4.M(), locMissingP4.M());
		}

		// both tracks have SC hit
                if(nFoundSC == 2){
                        dDeltaT3pi->Fill(locDeltaT[0],locDeltaT[1]);
			dMM_M3pi_2SC->Fill(locOmegaP4.M(), locMissingP4.M());
		}

		// select omega mass
		if(locOmegaP4.M() > 0.75 && locOmegaP4.M() < 0.875){
			
			// correlate unused tracks with 3pi system (already has pre-select)
			int nUnusedChargedTrackHypotheses = 0;
			for(size_t loc_j = 0; loc_j < locUnusedChargedTrackHypotheses.size(); ++loc_j) {
				if(locUnusedChargedTrackHypotheses[loc_j]->PID() != Proton) continue;
				nUnusedChargedTrackHypotheses++;
				DLorentzVector locUnusedP4 = locUnusedChargedTrackHypotheses[loc_j]->lorentzMomentum();
				double deltaPhi = (locUnusedP4.Vect().Phi() - locMissingP4.Vect().Phi())*180./TMath::Pi();
				if(deltaPhi > 180.) deltaPhi -= 360.;
				else if(deltaPhi < -180.) deltaPhi += 360.;
				dMissingDeltaPT_DeltaPhi->Fill(deltaPhi, locUnusedP4.Vect().Perp()-locMissingP4.Vect().Perp());
				
				// threshold on recoil momentum
				if(locMissingP4.Vect().Mag() > 0.5) {
					DLorentzVector locMissingTotalP4(0,0,locBeamPhotonEnergy,locBeamPhotonEnergy+0.938);
					locMissingTotalP4 -= locUnusedP4;
					locMissingTotalP4 -= locOmegaP4;
					dMM2_DeltaPhi->Fill(deltaPhi, locMissingTotalP4.M2());
				}
			}

			// mandelstam t
			TLorentzVector locGammaP4Init(0,0,locBeamPhotonEnergy,locBeamPhotonEnergy);
			TLorentzVector locProtonP4Init(0,0,0,0.938);
			TLorentzVector locSumP4Init = locGammaP4Init; locSumP4Init += locProtonP4Init;
			TLorentzVector locDelta  = (locMissingP4 - locProtonP4Init);
			double t = locDelta.M2();

			// theta CM
			TVector3 locBoostVector = -1.*locSumP4Init.BoostVector();
			TLorentzVector locGammaP4CMFrame = locGammaP4Init;
			TLorentzVector locOmegaP4CMFrame = locOmegaP4;
			locGammaP4CMFrame.Boost(locBoostVector);
			locOmegaP4CMFrame.Boost(locBoostVector);
			
			double locThetaCM = locGammaP4CMFrame.Vect().Angle(locOmegaP4CMFrame.Vect());
			dCosThetaCM_t->Fill(TMath::Cos(locThetaCM), -1.*t);
		}
	}
	japp->RootUnLock(); //RELEASE ROOT LOCK!!

	return true; //return false if you want to use this action to apply a cut (and it fails the cut!)
}